# Source code for xga.sources.general

#  This code is a part of X-ray: Generate and Analyse (XGA), a module designed for the XMM Cluster Survey (XCS).

from typing import Tuple, List, Union
from warnings import warn, simplefilter

import numpy as np
from astropy import wcs
from astropy.coordinates import SkyCoord
from astropy.cosmology import Planck15
from astropy.units import Quantity, UnitBase, deg, UnitConversionError
from numpy import ndarray

from .base import BaseSource
from ..exceptions import NotAssociatedError, PeakConvergenceFailedError, NoRegionsError, NoValidObservationsError, \
NoProductAvailableError
from ..products import RateMap

# This disables an annoying astropy warning that pops up all the time with XMM images
# Don't know if I should do this really
simplefilter('ignore', wcs.FITSFixedWarning)

[docs]class ExtendedSource(BaseSource):
"""
The general extended source XGA class, for extended X-ray sources that you do not yet want to analyse as
a specific type of astrophysical object. This class is subclassed by GalaxyCluster, which then adds more
specific analyses, for instance. This general class is useful for when you're doing exploratory analyses.

:param float ra: The right-ascension of the source, in degrees.
:param float dec: The declination of the source, in degrees.
:param float redshift: The redshift of the source, optional. Default is None.
:param str name: The name of the source, optional. Name will be constructed from position if None.
:param Quantity custom_region_radius: A custom analysis region radius for this source, optional.
:param bool use_peak: Whether peak position should be found and used.
:param Quantity peak_lo_en: The lower energy bound for the RateMap to calculate peak position
from. Default is 0.5keV
:param Quantity peak_hi_en: The upper energy bound for the RateMap to calculate peak position
from. Default is 2.0keV.
:param float back_inn_rad_factor: This factor is multiplied by an analysis region radius, and gives the inner
radius for the background region. Default is 1.05.
:param float back_out_rad_factor: This factor is multiplied by an analysis region radius, and gives the outer
radius for the background region. Default is 1.5.
:param cosmology: An astropy cosmology object for use throughout analysis of the source.
:param str peak_find_method: Which peak finding method should be used (if use_peak is True). Default
is hierarchical, simple may also be passed.
:param bool in_sample: A boolean argument that tells the source whether it is part of a sample or not, setting
to True suppresses some warnings so that they can be displayed at the end of the sample progress bar. Default
is False. User should only set to True to remove warnings.
"""
def __init__(self, ra: float, dec: float, redshift: float = None, name: str = None,
custom_region_radius: Quantity = None, use_peak: bool = True,
peak_lo_en: Quantity = Quantity(0.5, "keV"), peak_hi_en: Quantity = Quantity(2.0, "keV"),
load_products: bool = True, load_fits: bool = False, peak_find_method: str = "hierarchical",
in_sample: bool = False):
"""
The init for the general extended source XGA class, takes information on the position (and optionally
redshift) of source of interest, matches to extended regions, and optionally performs peak finding.

:param float ra: The right-ascension of the source, in degrees.
:param float dec: The declination of the source, in degrees.
:param float redshift: The redshift of the source, optional. Default is None.
:param str name: The name of the source, optional. Name will be constructed from position if None.
:param Quantity custom_region_radius: A custom analysis region radius for this source, optional.
:param bool use_peak: Whether peak position should be found and used.
:param Quantity peak_lo_en: The lower energy bound for the RateMap to calculate peak position
from. Default is 0.5keV
:param Quantity peak_hi_en: The upper energy bound for the RateMap to calculate peak position
from. Default is 2.0keV.
:param float back_inn_rad_factor: This factor is multiplied by an analysis region radius, and gives the inner
radius for the background region. Default is 1.05.
:param float back_out_rad_factor: This factor is multiplied by an analysis region radius, and gives the outer
radius for the background region. Default is 1.5.
:param cosmology: An astropy cosmology object for use throughout analysis of the source.
:param str peak_find_method: Which peak finding method should be used (if use_peak is True). Default
is hierarchical, simple may also be passed.
:param bool in_sample: A boolean argument that tells the source whether it is part of a sample or not, setting
to True suppresses some warnings so that they can be displayed at the end of the sample progress bar. Default
is False. User should only set to True to remove warnings.
"""
# Calling the BaseSource init method

# Setting up the custom region radius attributes

# Adding a custom radius to act as a search aperture for peak finding
# 500kpc in degrees, for the current redshift and cosmology
#  Or 5 arcminutes if no redshift information is present (that is allowed for the ExtendedSource class)
if self._redshift is not None:
search_aperture = rad_to_ang(Quantity(500, "kpc"), self._redshift, cosmo=self._cosmo)
else:
search_aperture = Quantity(5, 'arcmin').to('deg')

self._use_peak = use_peak
# Make sure the peak energy boundaries are in keV
self._peak_lo_en = peak_lo_en.to('keV')
self._peak_hi_en = peak_hi_en.to('keV')
self._peaks = {o: {} for o in self.obs_ids}
self._peaks.update({"combined": None})
self._peaks_near_edge = {o: {} for o in self.obs_ids}
self._peaks_near_edge.update({"combined": None})
self._chosen_peak_cluster = None
self._other_peak_clusters = None
self._snr = {}

# This uses the added context of the type of source to find (or not find) matches in region files
self._regions, self._alt_match_regions, self._other_regions = self._source_type_match("ext")

# Making a combined list of interloper regions
self._interloper_regions = []
for o in self._other_regions:
self._interloper_regions += self._other_regions[o]

# Run through any alternative matches and raise warnings
for o in self._alt_match_regions:
if len(self._alt_match_regions[o]) > 0:
warn_text = "There are {0} alternative matches for observation {1}, associated with " \
"source {2}".format(len(self._alt_match_regions[o]), o, self.name)
if not self._samp_member:
warn(warn_text, stacklevel=2)
else:
self._supp_warn.append(warn_text)

for obs_id in self.obs_ids:
# Generating and storing these because they should only
cur_im = self.get_products("image", obs_id)[0]

# Constructs the detected dictionary, detailing whether the source has been detected IN REGION FILES
#  in each observation.
self._detected = {o: self._regions[o] is not None for o in self._regions}

# If in some of the observations the source has not been detected, a warning will be raised
if True in self._detected.values() and False in self._detected.values():
warn_text = "{n} has not been detected in all region files, so generating and fitting products" \
" with the 'region' reg_type will not use all available data".format(n=self.name)
if not self._samp_member:
warn(warn_text, stacklevel=2)
else:
self._supp_warn.append(warn_text)

# If the source wasn't detected in ALL of the observations, then we have to rely on a custom region,
#  and if no custom region options are passed by the user then an error is raised.
elif all([det is False for det in self._detected.values()]) and self._custom_region_radius is not None:
warn_text = "{n} has not been detected in ANY region files, so generating and fitting products" \
" with the 'region' reg_type will not work".format(n=self.name)
if not self._samp_member:
warn(warn_text, stacklevel=2)
else:
self._supp_warn.append(warn_text)
elif all([det is False for det in self._detected.values()]) and self._custom_region_radius is None \
and "GalaxyCluster" not in repr(self):
raise NoRegionsError("{n} has not been detected in ANY region files, and no custom region or "
" overdensity radius has been passed. No analysis is possible.".format(n=self.name))

# Call to a method that goes through all the observations and finds the X-ray peak. Also at the same
#  time finds the X-ray peak of the combined ratemap (an essential piece of information).
self._all_peaks(peak_find_method)
if self._use_peak:
self._default_coord = self.peak

[docs]    def find_peak(self, rt: RateMap, method: str = "hierarchical", num_iter: int = 20, peak_unit: UnitBase = deg) \
-> Tuple[Quantity, bool, bool, ndarray, List]:
"""
A method that will find the X-ray peak for the RateMap that has been passed in. It takes
the user supplied coordinates from source initialisation as a starting point, finds the peak within a 500kpc
radius, re-centres the region, and iterates until the peak converges to within 15kpc, or until 20
20 iterations has been reached.

:param RateMap rt: The ratemap which we want to find the peak (local to our user supplied coordinates) of.
:param str method: Which peak finding method to use. Currently either hierarchical or simple can be chosen.
:param int num_iter: How many iterations should be allowed before the peak is declared as not converged.
:param UnitBase peak_unit: The unit the peak coordinate is returned in.
:return: The peak coordinate, a boolean flag as to whether the returned coordinates are near
a chip gap/edge, and a boolean flag as to whether the peak converged. It also returns the coordinates
of the points within the chosen point cluster, and a list of all point clusters that were not chosen.
:rtype: Tuple[Quantity, bool, bool, ndarray, List]
"""
all_meth = ["hierarchical", "simple"]
if method not in all_meth:
raise ValueError("{0} is not a recognised, use one of the following methods: "
"{1}".format(method, ", ".join(all_meth)))
central_coords = SkyCoord(*self.ra_dec.to("deg"))

# Iteration counter just to kill it if it doesn't converge
count = 0
# Allow 20 iterations by default before we kill this - alternatively loop will exit when centre converges
#  to within 15kpc (or 0.15arcmin).
while count < num_iter:

# Find the peak using the experimental clustering_peak method
if method == "hierarchical":
try:
peak, near_edge, chosen_coords, other_coords = rt.clustering_peak(aperture_mask, peak_unit)
except ValueError:
raise PeakConvergenceFailedError("The hierarchical clustering peak finder does not "
"have enough points to work with.")
elif method == "simple":
chosen_coords = []
other_coords = []

peak_deg = rt.coord_conv(peak, deg)
# Calculate the distance between new peak and old central coordinates
separation = Quantity(np.sqrt(abs(peak_deg[0].value - central_coords.ra.value) ** 2 +
abs(peak_deg[1].value - central_coords.dec.value) ** 2), deg)

central_coords = SkyCoord(*peak_deg.copy())
if self._redshift is not None:

if count != 0 and self._redshift is not None and separation <= Quantity(15, "kpc"):
break
elif count != 0 and self._redshift is None and separation <= Quantity(0.15, 'arcmin'):
break

count += 1

if count == num_iter:
converged = False
# To do the least amount of damage, if the peak doesn't converge then we just return the
#  user supplied coordinates
peak = self.ra_dec
near_edge = rt.near_edge(peak)
else:
converged = True

return peak, near_edge, converged, chosen_coords, other_coords

def _all_peaks(self, method: str):
"""
An internal method that finds the X-ray peaks for all of the available observations and instruments,
as well as the combined ratemap. Peak positions for individual ratemap products are allowed to not
converge, and will just write None to the peak dictionary, but if the peak of the combined ratemap fails
to converge an error will be thrown. The combined ratemap peak will also be stored by itself in an
attribute, to allow a property getter easy access.

:param str method: The method to be used for peak finding.
"""
en_key = "bound_{l}-{u}".format(l=self._peak_lo_en.value, u=self._peak_hi_en.value)
comb_rt = [rt[-1] for rt in self.get_products("combined_ratemap", just_obj=False) if en_key in rt]

if len(comb_rt) != 0:
comb_rt = comb_rt[0]
else:
# I didn't want to import this here, but otherwise circular imports become a problem
from xga.sas import emosaic
emosaic(self, "image", self._peak_lo_en, self._peak_hi_en, disable_progress=True)
emosaic(self, "expmap", self._peak_lo_en, self._peak_hi_en, disable_progress=True)
comb_rt = [rt[-1] for rt in self.get_products("combined_ratemap", just_obj=False) if en_key in rt][0]

if self._use_peak:
coord, near_edge, converged, cluster_coords, other_coords = self.find_peak(comb_rt, method=method)
# Updating nH for new coord, probably won't make a difference most of the time
self._nH = nh_lookup(coord)[0]
else:
# If we don't care about peak finding then this is the one to go for
coord = self.ra_dec
near_edge = comb_rt.near_edge(coord)
converged = True
cluster_coords = np.ndarray([])
other_coords = []

# Unfortunately if the peak convergence fails for the combined ratemap I have to raise an error
if converged:
self._peaks["combined"] = coord
self._peaks_near_edge["combined"] = near_edge
# I'm only going to save the point cluster positions for the combined ratemap
self._chosen_peak_cluster = cluster_coords
self._other_peak_clusters = other_coords
else:
raise PeakConvergenceFailedError("Peak finding on the combined ratemap failed to converge within "
"15kpc for {n} in the {l}-{u} energy "
"band.".format(n=self.name, l=self._peak_lo_en, u=self._peak_hi_en))

# for obs in self.obs_ids:
#     for rt in self.get_products("ratemap", obs_id=obs, extra_key=en_key, just_obj=True):
#         if self._use_peak:
#             coord, near_edge, converged, cluster_coords, other_coords = self.find_peak(rt)
#             if converged:
#                 self._peaks[obs][rt.instrument] = coord
#                 self._peaks_near_edge[obs][rt.instrument] = near_edge
#             else:
#                 self._peaks[obs][rt.instrument] = None
#                 self._peaks_near_edge[obs][rt.instrument] = None
#         else:
#             self._peaks[obs][rt.instrument] = self.ra_dec
#             self._peaks_near_edge[obs][rt.instrument] = rt.near_edge(self.ra_dec)

[docs]    def get_peaks(self, obs_id: str = None, inst: str = None) -> Quantity:
"""
A get method to return the peak of the X-ray emission of this GalaxyCluster.

:param str obs_id: The ObsID to return the X-ray peak coordinates for.
:param str inst: The instrument to return the X-ray peak coordinates for.
:return: The X-ray peak coordinates for the input parameters.
:rtype: Quantity
"""
# Common sense checks, are the obsids/instruments associated with this source etc.
if obs_id is not None and obs_id not in self.obs_ids:
raise NotAssociatedError("The ObsID {o} is not associated with {s}.".format(o=obs_id, s=self.name))
elif obs_id is None and inst is not None:
raise ValueError("If obs_id is None, inst cannot be None as well.")
elif obs_id is not None and inst is not None and inst not in self._peaks[obs_id]:
raise NotAssociatedError("The instrument {i} is not associated with observation {o} of this "
"source.".format(i=inst, o=obs_id))
elif obs_id is None and inst is None:
chosen = self._peaks
elif obs_id is not None and inst is None:
chosen = self._peaks[obs_id]
else:
chosen = self._peaks[obs_id][inst]

return chosen

[docs]    def get_1d_brightness_profile(self, outer_rad: Union[Quantity, str], obs_id: str = 'combined',
inst: str = 'combined', central_coord: Quantity = None, radii: Quantity = None,
lo_en: Quantity = None, hi_en: Quantity = None, pix_step: int = 1,
min_snr: Union[float, int] = 0.0, psf_corr: bool = False, psf_model: str = "ELLBETA",
psf_bins: int = 4, psf_algo: str = "rl", psf_iter: int = 15):
"""
A specific get method for 1D brightness profiles. Should provide a relatively simple way of retrieving
specific brightness profiles from XGA's storage system. Please note that there is not a separate get method
for brightness profiles made from combined data, instead this method will search for combined profiles if
either obs_id or inst is set to 'combined'. - Retrieving combined profiles is the default behaviour.

:param Quantity/str outer_rad: The outermost radius of the profile, either as a Quantity or a name (e.g. r500).
:param str obs_id: The ObsID used to generate the profile in question, default is None. If this is set to
combined then this method will search for profiles based on combined data.
:param str inst: The instrument used to generate the profile in question, default is None. If this is set to
combined then this method will search for profiles based on combined data.
:param Quantity central_coord: The central coordinate from which the profile was generated. Default is
None, which means we shall use the default coordinate of this source.
:param Quantity lo_en: The lower energy bound of the RateMap used to generate the profile.
:param Quantity hi_en: The upper energy bound of the RateMap used to generate the profile.
:param int pix_step: The width of each annulus in pixels used to generate the profile.
:param float min_snr: The minimum signal to noise imposed upon the profile.
:param bool psf_corr: Is the brightness profile corrected for PSF effects?
:param str psf_model: If PSF corrected, the PSF model used.
:param int psf_bins: If PSF corrected, the number of bins per side.
:param str psf_algo: If PSF corrected, the algorithm used.
:param int psf_iter: If PSF corrected, the number of algorithm iterations.
:return:
"""
# Makes sure its in our standard unit
else:
raise ValueError("Outer radius may only be a string or an astropy quantity")

if obs_id == "combined" or inst == "combined":
interim_prods = self.get_combined_profiles("brightness", central_coord, radii, lo_en, hi_en)
else:
interim_prods = self.get_profiles("brightness", obs_id, inst, central_coord, radii, lo_en, hi_en)

# The methods I used to get this far will already have gotten upset if there are no matches, so I don't need
#  to check they exist, but I do need to check if I have a list or a single object
if not isinstance(interim_prods, list):
interim_prods = [interim_prods]

matched_prods = []
for p in interim_prods:
if not psf_corr and p.outer_radius == outer_rad and p.pix_step == pix_step and p.min_snr == min_snr \
and p.psf_corrected == psf_corr:
matched_prods.append(p)
elif psf_corr and p.outer_radius == outer_rad and p.pix_step == pix_step and p.min_snr == min_snr \
and p.psf_corrected == psf_corr and p.psf_model == psf_model and p.psf_bins == psf_bins \
and p.psf_algorithm == psf_algo and p.psf_iterations == psf_iter:
matched_prods.append(p)

if len(matched_prods) == 1:
matched_prods = matched_prods[0]
elif len(matched_prods) == 0:
raise NoProductAvailableError("Cannot find any brightness profiles matching your input.")

return matched_prods

# Property SPECIFICALLY FOR THE COMBINED PEAK - as this is the peak we should be using mostly.
@property
def peak(self) -> Quantity:
"""
A property getter for the combined X-ray peak coordinates. Most analysis will be centered
on these coordinates.

:return: The X-ray peak coordinates for the combined ratemap.
:rtype: Quantity
"""
return self._peaks["combined"]

# I'm allowing this a setter, as some users may want to update the peak from outside (as is
#  done in the ClusterSample init).
@peak.setter
def peak(self, new_peak: Quantity):
"""
Allows the user to update the peak value used during analyses manually.

:param Quantity new_peak: A new RA-DEC peak coordinate, in degrees.
"""
if not new_peak.unit.is_equivalent("deg"):
raise UnitConversionError("The new peak value must be in RA and DEC coordinates")
elif len(new_peak) != 2:
raise ValueError("Please pass an astropy Quantity, in units of degrees, with two entries - "
"one for RA and one for DEC.")
self._peaks["combined"] = new_peak.to("deg")

@property
"""
A getter for the custom region that can be defined on initialisation.

:return: The radius (in kpc) of the user defined custom region.
:rtype: Quantity
"""

@property
def point_clusters(self) -> Tuple[ndarray, List[ndarray]]:
"""
This allows you to retrieve the point cluster positions from the hierarchical clustering
peak finding method run on the combined ratemap. This includes both the chosen cluster and
all others that were found.

:return: A numpy array of the positions of points of the chosen cluster (not galaxy cluster,
a cluster of points). A list of numpy arrays with the same information for all the other clusters
that were found
:rtype: Tuple[ndarray, List[ndarray]]
"""
return self._chosen_peak_cluster, self._other_peak_clusters

[docs]class PointSource(BaseSource):
"""
The general point source XGA class, for point X-ray sources that you do not yet want to analyse as
a specific type of astrophysical object. This general class is useful for when you're doing exploratory analyses.

:param float ra: The right-ascension of the point source, in degrees.
:param float dec: The declination of the point source, in degrees.
:param float redshift: The redshift of the point source, optional. Default is None.
:param str name: The name of the point source, optional. If no names are supplied
then they will be constructed from the supplied coordinates.
:param Quantity point_radius: The point source analysis region radius for this sample. An astropy quantity
containing the radius should be passed; remember that units like kpc will also need redshift
information. Default is 30 arcsecond radius.
:param bool use_peak: Whether peak position should be found and used. For PointSource the 'simple' peak
finding method is the only one available, other methods are allowed for extended sources.
:param Quantity peak_lo_en: The lower energy bound for the RateMap to calculate peak
position from. Default is 0.5keV.
:param Quantity peak_hi_en: The upper energy bound for the RateMap to calculate peak
position from. Default is 2.0keV.
:param float back_inn_rad_factor: This factor is multiplied by an analysis region radius, and gives the inner
radius for the background region. Default is 1.05.
:param float back_out_rad_factor: This factor is multiplied by an analysis region radius, and gives the outer
radius for the background region. Default is 1.5.
:param cosmology: An astropy cosmology object for use throughout analysis of the source.
:param bool regen_merged: Should merged images/exposure maps be regenerated after cleaning. Default is
True. This option is here so that sample objects can regenerate all merged products at once, which is
more efficient as it can exploit parallelisation more fully - user probably doesn't need to touch this.
:param bool in_sample: A boolean argument that tells the source whether it is part of a sample or not, setting
to True suppresses some warnings so that they can be displayed at the end of the sample progress bar. Default
is False. User should only set to True to remove warnings.
"""
def __init__(self, ra, dec, redshift=None, name=None, point_radius=Quantity(30, 'arcsec'), use_peak=False,
regen_merged: bool = True, in_sample: bool = False):
"""
The init of the general XGA point source class.

:param float ra: The right-ascension of the point source, in degrees.
:param float dec: The declination of the point source, in degrees.
:param float redshift: The redshift of the point source, optional. Default is None.
:param str name: The name of the point source, optional. If no names are supplied
then they will be constructed from the supplied coordinates.
:param Quantity point_radius: The point source analysis region radius for this sample. An astropy quantity
containing the radius should be passed; remember that units like kpc will also need redshift
information. Default is 30 arcsecond radius.
:param bool use_peak: Whether peak position should be found and used. For PointSource the 'simple' peak
finding method is the only one available, other methods are allowed for extended sources.
:param Quantity peak_lo_en: The lower energy bound for the RateMap to calculate peak
position from. Default is 0.5keV.
:param Quantity peak_hi_en: The upper energy bound for the RateMap to calculate peak
position from. Default is 2.0keV.
:param float back_inn_rad_factor: This factor is multiplied by an analysis region radius, and gives the inner
radius for the background region. Default is 1.05.
:param float back_out_rad_factor: This factor is multiplied by an analysis region radius, and gives the outer
radius for the background region. Default is 1.5.
:param cosmology: An astropy cosmology object for use throughout analysis of the source.
:param bool regen_merged: Should merged images/exposure maps be regenerated after cleaning. Default is
True. This option is here so that sample objects can regenerate all merged products at once, which is
more efficient as it can exploit parallelisation more fully - user probably doesn't need to touch this.
:param bool in_sample: A boolean argument that tells the source whether it is part of a sample or not, setting
to True suppresses some warnings so that they can be displayed at the end of the sample progress bar. Default
is False. User should only set to True to remove warnings.
"""

# This uses the added context of the type of source to find (or not find) matches in region files
# This is the internal dictionary where all regions, defined by reg-files or by users, will be stored
self._regions, self._alt_match_regions, self._other_regions = self._source_type_match("pnt")
self._detected = {o: self._regions[o] is not None for o in self._regions}

# Making a combined list of interloper regions
self._interloper_regions = []
for o in self._other_regions:
self._interloper_regions += self._other_regions[o]

if self._redshift is not None and point_radius.unit.is_equivalent("kpc"):
else:
raise UnitConversionError("Can't convert {u} to a XGA supported length unit".format(u=point_radius.unit))

# This generates masks to remove interloper regions
for obs_id in self.obs_ids:
# Generating and storing these because they should only
cur_im = self.get_products("image", obs_id)[0]

# Here we automatically clean the observations, to make sure the point source does actually lie
#  on the detector and not just near it
# Use a pretty harsh acceptance fraction
reject_dict = self.obs_check("point", 0.9)
if len(reject_dict) != 0:
# Use the source method to remove data we've decided isn't worth keeping
self.disassociate_obs(reject_dict)
if len(self._obs) == 0:
raise NoValidObservationsError("Observation cleaning has been run and there are no remaining"
" observations. ")

if regen_merged:
from ..sas import emosaic
emosaic(self, "image", self._peak_lo_en, self._peak_hi_en, disable_progress=True)
emosaic(self, "expmap", self._peak_lo_en, self._peak_hi_en, disable_progress=True)

self._use_peak = use_peak
# Make sure the peak energy boundaries are in keV
self._peak_lo_en = peak_lo_en.to('keV')
self._peak_hi_en = peak_hi_en.to('keV')
self._peaks = {o: {} for o in self.obs_ids}
self._peaks.update({"combined": None})
self._peaks_near_edge = {o: {} for o in self.obs_ids}

self._all_peaks()

if self._use_peak:
self._default_coord = self.peak

@property
"""
Property getter to access the point_radius declared on initialisation of the source, the radius
of the region that is used for point source analysis.

:return: The radius of the point source analysis region.
:rtype: Quantity
"""

def _all_peaks(self):
en_key = "bound_{l}-{u}".format(l=self._peak_lo_en.value, u=self._peak_hi_en.value)
comb_rt = self.get_products("combined_ratemap", extra_key=en_key)

if len(comb_rt) != 0:
comb_rt = comb_rt[0]
else:
from xga.sas import emosaic
emosaic(self, "image", self._peak_lo_en, self._peak_hi_en, disable_progress=True)
emosaic(self, "expmap", self._peak_lo_en, self._peak_hi_en, disable_progress=True)
comb_rt = self.get_products("combined_ratemap", extra_key=en_key)[0]

if self._use_peak:
coord, near_edge = self.find_peak(comb_rt)
# Updating nH for new coord, probably won't make a difference most of the time
self._nH = nh_lookup(coord)[0]
else:
# If we don't care about peak finding then this is the boi to go for
coord = self.ra_dec
near_edge = comb_rt.near_edge(coord)

self._peaks["combined"] = coord
self._peaks_near_edge["combined"] = near_edge

[docs]    def find_peak(self, rt: RateMap, peak_unit: UnitBase = deg) -> Tuple[Quantity, bool]:
"""
Uses a simple 'brightest pixel' method to measure a peak coordinate for the point source.

:param RateMap rt: The RateMap to measure the peak from.
:param UnitBase peak_unit: The desired output unit of the peak.
:return:  The peak, and a boolean flag as to whether the peak is near an edge.
:rtype: Tuple[Quantity, bool]
"""
central_coords = SkyCoord(*self.ra_dec.to("deg"))

return peak, near_edge

@property
def peak(self) -> Quantity:
"""
A property getter for the combined X-ray peak coordinates.

:return: The X-ray peak coordinates for the combined ratemap.
:rtype: Quantity
"""
return self._peaks["combined"]