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.
DavidT3 · Sep 30, 2020 · c1e2da8 · c1e2da8
1 parent fdee4e6
commit c1e2da8
Show file tree

Hide file tree

Showing 6 changed files with 78 additions and 30 deletions.
diff --git a/experiments/deprojection_experiments.ipynb b/experiments/deprojection_experiments.ipynb
@@ -209,15 +209,15 @@
     "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",
     "# Plot the background level\n",
     "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
@@ -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
@@ -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 ([email protected]) 28/09/2020, 14:55. Copyright (c) David J Turner
+#  Last modified by David J Turner ([email protected]) 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,15 +262,18 @@ 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])
             inn_rads = np.delete(inn_rads, below[0])
 
         # 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
@@ -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 ([email protected]) 28/09/2020, 15:23. Copyright (c) David J Turner
+#  Last modified by David J Turner ([email protected]) 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
@@ -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 ([email protected]) 28/09/2020, 12:44. Copyright (c) David J Turner
+#  Last modified by David J Turner ([email protected]) 30/09/2020, 12:28. Copyright (c) David J Turner
 
 import warnings
 from typing import Tuple, Union
@@ -254,24 +254,33 @@ 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:
                 p_rt = psf_comb_rt[-1]
                 # 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
@@ -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 ([email protected]) 24/09/2020, 17:35. Copyright (c) David J Turner
+#  Last modified by David J Turner ([email protected]) 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