From c1e2da87311864b1caaaaf4b30335c1182683a46 Mon Sep 17 00:00:00 2001 From: dt237 Date: Wed, 30 Sep 2020 12:28:01 +0100 Subject: [PATCH] Added the calculation of uncertainties to radial_brightness_profile, and updated all dependants so that they support the new return. Closes issue #172 (for now - might change error calculation method as mentioned in the comments of #172). view_brightness_profile now plots the uncertainties. The RateMap class now stores the input Image and ExpMap objects in attributes, and can return them as properties. --- experiments/deprojection_experiments.ipynb | 4 +-- experiments/xga_psf_experiments.ipynb | 8 +++--- xga/imagetools/profile.py | 28 ++++++++++++++---- xga/products/phot.py | 25 ++++++++++++++-- xga/sources/extended.py | 33 ++++++++++++++-------- xga/sourcetools/stack.py | 10 ++++--- 6 files changed, 78 insertions(+), 30 deletions(-) diff --git a/experiments/deprojection_experiments.ipynb b/experiments/deprojection_experiments.ipynb index 947764cd..f34c2d11 100644 --- a/experiments/deprojection_experiments.ipynb +++ b/experiments/deprojection_experiments.ipynb @@ -209,7 +209,7 @@ "plt.figure(figsize=(8, 5))\n", "ax = plt.gca()\n", "\n", - "og_brightness, radii, rad_err, og_background = radial_brightness(a907_rt, source_mask, background_mask, pix_peak,\n", + "og_brightness, og_br_err, radii, rad_err, og_background = radial_brightness(a907_rt, source_mask, background_mask, pix_peak,\n", " rad, a907.redshift, 1, kpc, a907.cosmo)\n", "plt.plot(radii, og_brightness, label=\"Emission\")\n", "\n", @@ -217,7 +217,7 @@ "plt.axhline(og_background, color=\"black\", linestyle=\"dashed\", label=\"Background\")\n", "\n", "psf_peak = a907_psf_rt.clustering_peak(source_mask, pix)[0]\n", - "psf_brightness, radii, rad_err, psf_background = radial_brightness(a907_psf_rt, source_mask, background_mask, psf_peak,\n", + "psf_brightness, br_err, radii, rad_err, psf_background = radial_brightness(a907_psf_rt, source_mask, background_mask, psf_peak,\n", " rad, a907.redshift, 1, kpc, a907.cosmo)\n", "plt.plot(radii, psf_brightness, label=\"PSF Corr Emission\")\n", "\n", diff --git a/experiments/xga_psf_experiments.ipynb b/experiments/xga_psf_experiments.ipynb index 3bb34bc2..1ee858db 100644 --- a/experiments/xga_psf_experiments.ipynb +++ b/experiments/xga_psf_experiments.ipynb @@ -915,13 +915,13 @@ "# Centering at the original peak used to generate the PSFs\n", "pix_loc = a907_mos1_im.coord_conv(a907_mos1_psf.ra_dec, pix)\n", "# Getting the radial profile for the normal stacked ratemaps\n", - "brightness, radii, background, inn_r, out_r = radial_brightness(a907_comb_rt, flat_mask, flat_mask, pix_loc,\n", + "brightness, br_err, radii, background, inn_r, out_r = radial_brightness(a907_comb_rt, flat_mask, flat_mask, pix_loc,\n", " Quantity(1000, 'kpc'), 0.16, 1, kpc)\n", "brightness /= pix_kpc.value**2\n", "plt.plot(radii, brightness, label=\"Original Brightness Profile\", color=\"tab:blue\")\n", "\n", "# Getting the radial profile for the deconvolved stacked ratemaps\n", - "brightness, radii, background, inn_r, out_r = radial_brightness(a907_comb_deconv_rt, flat_mask, flat_mask, pix_loc,\n", + "brightness, br_err, radii, background, inn_r, out_r = radial_brightness(a907_comb_deconv_rt, flat_mask, flat_mask, pix_loc,\n", " Quantity(1000, 'kpc'), 0.16, 1, kpc)\n", "brightness /= pix_kpc.value**2\n", "\n", @@ -1339,13 +1339,13 @@ "# Centering at the original peak used to generate the PSFs\n", "pix_loc = off_mos1_im.coord_conv(off_mos1_psf.ra_dec, pix)\n", "# Getting the radial profile for the normal stacked ratemaps\n", - "brightness, radii, background, inn_r, out_r = radial_brightness(off_comb_rt, flat_mask, flat_mask, pix_loc,\n", + "brightness, br_err, radii, background, inn_r, out_r = radial_brightness(off_comb_rt, flat_mask, flat_mask, pix_loc,\n", " Quantity(1000, 'kpc'), 0.13500871, 1, kpc)\n", "brightness /= pix_kpc.value**2\n", "plt.plot(radii, brightness, label=\"Original Brightness Profile\", color=\"tab:blue\")\n", "\n", "# Getting the radial profile for the deconvolved stacked ratemaps\n", - "brightness, radii, background, inn_r, out_r = radial_brightness(off_comb_deconv_rt, flat_mask, flat_mask, pix_loc,\n", + "brightness, br_err, radii, background, inn_r, out_r = radial_brightness(off_comb_deconv_rt, flat_mask, flat_mask, pix_loc,\n", " Quantity(1000, 'kpc'), 0.13500871, 1, kpc)\n", "brightness /= pix_kpc.value**2\n", "\n", diff --git a/xga/imagetools/profile.py b/xga/imagetools/profile.py index b308ce26..4cfc4c29 100644 --- a/xga/imagetools/profile.py +++ b/xga/imagetools/profile.py @@ -1,5 +1,5 @@ # This code is a part of XMM: Generate and Analyse (XGA), a module designed for the XMM Cluster Survey (XCS). -# Last modified by David J Turner (david.turner@sussex.ac.uk) 28/09/2020, 14:55. Copyright (c) David J Turner +# Last modified by David J Turner (david.turner@sussex.ac.uk) 30/09/2020, 12:28. Copyright (c) David J Turner from typing import Tuple @@ -152,7 +152,8 @@ def ann_radii(im_prod: Image, centre: Quantity, rad: Quantity, z: float = None, def radial_brightness(rt: RateMap, src_mask: np.ndarray, back_mask: np.ndarray, centre: Quantity, rad: Quantity, z: float = None, pix_step: int = 1, cen_rad_units: UnitBase = arcsec, - cosmo=Planck15, min_snr: float = 0.0) -> Tuple[np.ndarray, Quantity, Quantity, np.float64, bool]: + cosmo=Planck15, min_snr: float = 0.0) \ + -> Tuple[np.ndarray, np.ndarray, Quantity, Quantity, np.float64, bool]: """ A simple method to calculate the average brightness in circular annuli upto the radius of the chosen region. The annuli are one pixel in width, and as this uses the masks that were generated @@ -172,12 +173,16 @@ def radial_brightness(rt: RateMap, src_mask: np.ndarray, back_mask: np.ndarray, :return: The brightness is returned in a flat numpy array, then the radii at the centre of the bins are returned in units of kpc, the width of the bins, and finally the average brightness in the background region is returned. - :rtype: Tuple[np.ndarray, Quantity, Quantity, np.float64, bool] + :rtype: Tuple[np.ndarray, np.ndarray, Quantity, Quantity, np.float64, bool] """ if rt.shape != src_mask.shape: raise ValueError("The shape of the src_mask array ({0}) must be the same as that of im_prod " "({1}).".format(src_mask.shape, rt.shape)) + # TODO REMOVE THIS IF SHIT - trying to make an uncertainty map + cr_err_map = np.divide(np.sqrt(rt.image.data), rt.expmap.data, out=np.zeros_like(rt.image.data), + where=rt.expmap.data != 0) + # Returns conversion factor to degrees, so multiplying by 60 goes to arcminutes # Getting this because we want to be able to convert pixel distance into arcminutes for dividing by the area to_arcmin = pix_deg_scale(centre, rt.radec_wcs) * 60 @@ -202,14 +207,19 @@ def radial_brightness(rt: RateMap, src_mask: np.ndarray, back_mask: np.ndarray, # dimension. Helps with broadcasting the annular masks with the region src_mask that gets rid of interlopers masks = annular_mask(pix_cen, inn_rads, out_rads, rt.shape) * src_mask[..., None] * rt.sensor_mask[..., None] # This calculates the area of each annulus mask - areas = np.sum(masks, axis=(0, 1)) * to_arcmin**2 + num_pix = np.sum(masks, axis=(0, 1)) + areas = num_pix * to_arcmin**2 # Creates a 3D array of the masked data masked_data = masks * rt.data[..., None] # Calculates the sum of the pixel count rates for each annular radius, masking out other known sources cr = np.sum(masked_data, axis=(0, 1)) + # TODO REMOVE IF SHIT + cr_errs = np.sqrt(np.sum(masks * cr_err_map[..., None]**2, axis=(0, 1))) + # Then calculates the actual brightness profile by dividing by the area of each annulus br = cr / areas + br_errs = cr_errs / areas # Calculates the signal to noise profile, defined as the ratio between brightness profile and background snr_prof = br / bg @@ -220,6 +230,7 @@ def radial_brightness(rt: RateMap, src_mask: np.ndarray, back_mask: np.ndarray, # Making copies in case rebinning fails init_br = br.copy() + init_br_errs = br_errs.copy() init_inn = inn_rads.copy() init_out = out_rads.copy() @@ -232,9 +243,12 @@ def radial_brightness(rt: RateMap, src_mask: np.ndarray, back_mask: np.ndarray, # We deal with and modify the count rates and areas separately as you have to combine bins in the count # rate and area regimes separately, then calculate brightness by dividing cr by area. cr[below[0]] = cr[below[0]] + cr[below[0] + 1] + # ADDING ERRORS IN QUADRATURE FOR NOW, DON'T KNOW IF I'LL KEEP IT LIKE THIS + cr_errs[below[0]] = np.sqrt(cr_errs[below[0]]**2 + cr_errs[below[0] + 1]**2) areas[below[0]] = areas[below[0]] + areas[below[0] + 1] # Then the donor bin that was added to the problem bin is deleted cr = np.delete(cr, below[0] + 1) + cr_errs = np.delete(cr_errs, below[0] + 1) areas = np.delete(areas, below[0] + 1) # The new outer radius of the problem bin is set to the outer radius of the donor bin out_rads[below[0]] = out_rads[below[0] + 1] @@ -248,8 +262,10 @@ def radial_brightness(rt: RateMap, src_mask: np.ndarray, back_mask: np.ndarray, # As this is the last bin in the profile, there is no extra bin in the outward direction to add to our # problem bin. As such the problem bin is added inwards, to the N-1th bin in the profile cr[below[0] - 1] = cr[below[0] - 1] + cr[below[0]] + cr_errs[below[0] - 1] = np.sqrt(cr_errs[below[0] - 1]**2 + cr_errs[below[0]]**2) areas[below[0] - 1] = areas[below[0] - 1] + areas[below[0]] cr = np.delete(cr, below[0]) + cr_errs = np.delete(cr_errs, below[0]) areas = np.delete(areas, below[0]) out_rads[below[0] - 1] = out_rads[below[0]] out_rads = np.delete(out_rads, below[0]) @@ -257,6 +273,7 @@ def radial_brightness(rt: RateMap, src_mask: np.ndarray, back_mask: np.ndarray, # Calculate the new brightness with our combined count rates and areas br = cr/areas + br_errs = cr_errs / areas # Recalculate the SNR profile after the re-binning in this iteration snr_prof = br / bg # Find out which bins are still below the SNR threshold (if any) @@ -266,6 +283,7 @@ def radial_brightness(rt: RateMap, src_mask: np.ndarray, back_mask: np.ndarray, inn_rads = init_inn out_rads = init_out br = init_br + br_errs = init_br_errs succeeded = False else: succeeded = True @@ -275,7 +293,7 @@ def radial_brightness(rt: RateMap, src_mask: np.ndarray, back_mask: np.ndarray, cen_rads = (inn_rads + out_rads) / 2 rad_err = (out_rads-inn_rads) / 2 - return br, cen_rads, rad_err, bg, succeeded + return br, br_errs, cen_rads, rad_err, bg, succeeded # TODO At some point implement minimum SNR for this also diff --git a/xga/products/phot.py b/xga/products/phot.py index 59fd5bc8..8ea21c75 100644 --- a/xga/products/phot.py +++ b/xga/products/phot.py @@ -1,5 +1,5 @@ # This code is a part of XMM: Generate and Analyse (XGA), a module designed for the XMM Cluster Survey (XCS). -# Last modified by David J Turner (david.turner@sussex.ac.uk) 28/09/2020, 15:23. Copyright (c) David J Turner +# Last modified by David J Turner (david.turner@sussex.ac.uk) 30/09/2020, 12:28. Copyright (c) David J Turner import warnings @@ -598,6 +598,9 @@ def __init__(self, xga_image: Image, xga_expmap: ExpMap): self._im_path = xga_image.path self._expmap_path = xga_expmap.path + self._im_obj = xga_image + self._ex_obj = xga_expmap + self._im_data = None self._ex_data = None self._data = None @@ -706,8 +709,6 @@ def data(self) -> np.ndarray: self._construct_on_demand() return self._data - - def get_rate(self, at_coord: Quantity) -> float: """ A simple method that converts the given coordinates to pixels, then finds the rate (in photons @@ -974,6 +975,24 @@ def expmap_path(self) -> str: """ return self._expmap_path + @property + def image(self) -> Image: + """ + This property allows the user to access the input Image object for this ratemap. + :return: The input XGA Image object used to create this ratemap. + :rtype: Image + """ + return self._im_obj + + @property + def expmap(self) -> ExpMap: + """ + This property allows the user to access the input ExpMap object for this ratemap. + :return: The input XGA ExpMap object used to create this ratemap. + :rtype: ExpMap + """ + return self._ex_obj + class PSF(Image): def __init__(self, path: str, psf_model: str, obs_id: str, instrument: str, stdout_str: str, stderr_str: str, diff --git a/xga/sources/extended.py b/xga/sources/extended.py index 2376869f..74d63123 100644 --- a/xga/sources/extended.py +++ b/xga/sources/extended.py @@ -1,5 +1,5 @@ # This code is a part of XMM: Generate and Analyse (XGA), a module designed for the XMM Cluster Survey (XCS). -# Last modified by David J Turner (david.turner@sussex.ac.uk) 28/09/2020, 12:44. Copyright (c) David J Turner +# Last modified by David J Turner (david.turner@sussex.ac.uk) 30/09/2020, 12:28. Copyright (c) David J Turner import warnings from typing import Tuple, Union @@ -254,11 +254,16 @@ def view_brightness_profile(self, reg_type: str, profile_type: str = "radial", n if profile_type == "radial": ax.set_title("{n} - {l}-{u}keV Radial Brightness Profile".format(n=self.name, l=self._peak_lo_en.value, u=self._peak_hi_en.value)) - brightness, radii, radii_bins, og_background, success = radial_brightness(comb_rt, source_mask, - background_mask, pix_peak, rad, - self._redshift, pix_step, kpc, - self.cosmo, min_snr=min_snr) - prof = plt.errorbar(radii.value, brightness, xerr=radii_bins.value, fmt='x', capsize=2, label="Emission") + brightness, brightness_errs, radii, radii_bins, og_background, success = radial_brightness(comb_rt, + source_mask, + background_mask, + pix_peak, rad, + self._redshift, + pix_step, kpc, + self.cosmo, + min_snr=min_snr) + prof = plt.errorbar(radii.value, brightness, xerr=radii_bins.value, yerr=brightness_errs, fmt='x', + capsize=2, label="Emission") plt.plot(radii.value, brightness, color=prof[0].get_color()) for psf_comb_rt in psf_comb_rts: @@ -266,12 +271,16 @@ def view_brightness_profile(self, reg_type: str, profile_type: str = "radial", n # If the user wants to use individual peaks, we have to find them here. if not same_peak: pix_peak = self.find_peak(p_rt)[0] - brightness, radii, radii_bins, background, success = radial_brightness(psf_comb_rt[-1], source_mask, - background_mask, pix_peak, rad, - self._redshift, pix_step, kpc, - self.cosmo, min_snr=min_snr) - prof = plt.errorbar(radii.value, brightness, xerr=radii_bins.value, capsize=2, fmt="x", - label="{m} PSF Corrected".format(m=p_rt.psf_model)) + brightness, brightness_errs, radii, radii_bins, background, success = radial_brightness(psf_comb_rt[-1], + source_mask, + background_mask, + pix_peak, rad, + self._redshift, + pix_step, kpc, + self.cosmo, + min_snr=min_snr) + prof = plt.errorbar(radii.value, brightness, xerr=radii_bins.value, yerr=brightness_errs, capsize=2, + fmt="x", label="{m} PSF Corrected".format(m=p_rt.psf_model)) plt.plot(radii.value, brightness, color=prof[0].get_color()) plt.axhline(background, color=prof[0].get_color(), linestyle="dashed", diff --git a/xga/sourcetools/stack.py b/xga/sourcetools/stack.py index 79cc7fa5..b37b0c20 100644 --- a/xga/sourcetools/stack.py +++ b/xga/sourcetools/stack.py @@ -1,5 +1,5 @@ # This code is a part of XMM: Generate and Analyse (XGA), a module designed for the XMM Cluster Survey (XCS). -# Last modified by David J Turner (david.turner@sussex.ac.uk) 24/09/2020, 17:35. Copyright (c) David J Turner +# Last modified by David J Turner (david.turner@sussex.ac.uk) 30/09/2020, 12:28. Copyright (c) David J Turner from multiprocessing.dummy import Pool from typing import List, Tuple, Union @@ -94,9 +94,11 @@ def construct_profile(src: GalaxyCluster, src_id: int, lower: Quantity, upper: Q source_mask = src.get_interloper_mask() rad = Quantity(src.source_back_regions(scale_radius)[0].to_pixel(rt.radec_wcs).radius, pix) - brightness, cen_rad, rad_bins, bck, success = radial_brightness(rt, source_mask, background_mask, pix_peak, - rad, src.redshift, pix_step, pix, src.cosmo, - min_snr=min_snr) + brightness, brightness_err, cen_rad, rad_bins, bck, success = radial_brightness(rt, source_mask, + background_mask, pix_peak, + rad, src.redshift, pix_step, + pix, src.cosmo, + min_snr=min_snr) # Subtracting the background in the simplest way possible brightness -= bck