Merge pull request #1336 from DavidT3/feat/addNTempApecModel

Feat/add n temp apec model

Author: David Turner

xga/xspec/__init__.py

@@ -1,6 +1,6 @@
 #  This code is a part of X-ray: Generate and Analyse (XGA), a module designed for the XMM Cluster Survey (XCS).
-#  Last modified by David J Turner (turne540@msu.edu) 27/04/2023, 12:16. Copyright (c) The Contributors
+#  Last modified by David J Turner (turne540@msu.edu) 28/03/2025, 11:34. Copyright (c) The Contributors
 
 from .fit import single_temp_apec, power_law, single_temp_apec_profile, blackbody, multi_temp_dem_apec, \
-    single_temp_mekal
+    single_temp_mekal, double_temp_apec
 from .run import execute_cmd, xspec_call

xga/xspec/fit/general.py

@@ -1,5 +1,5 @@
 #  This code is a part of X-ray: Generate and Analyse (XGA), a module designed for the XMM Cluster Survey (XCS).
-#  Last modified by David J Turner (turne540@msu.edu) 27/03/2025, 11:43. Copyright (c) The Contributors
+#  Last modified by David J Turner (turne540@msu.edu) 29/03/2025, 05:01. Copyright (c) The Contributors
 
 import warnings
 from inspect import signature, Parameter
@@ -926,34 +926,26 @@ def blackbody(sources: Union[BaseSource, BaseSample], outer_radius: Union[str, Q
     return script_paths, outfile_paths, num_cores, run_type, src_inds, None, timeout, model, fit_confs, inv_ents
 
 
-# If the spectrum checking step of the XSPEC fit is enabled (using the boolean flag spectrum_checking), then
-#     each individual spectrum available for a given source will be fitted, and if the measured temperature is less
-#     than or equal to 0.01keV, or greater than 20keV, or the temperature uncertainty is greater than 15keV, then
-#     that spectrum will be rejected and not included in the final fit. Spectrum checking also involves rejecting any
-#     spectra with fewer than 10 noticed channels.
-#
-#     Freezing the temperature value of the fit is also possible, in cases where the data may not be sufficient to
-#     constrain it, and an external temperature constrain is used (by passing to the 'start_temp' argument).
-
-
 @xspec_call
-def n_temp_apec(sources: Union[BaseSource, BaseSample], outer_radius: Union[str, Quantity],
-                num_temp_comp: Union[int, List[int]], inner_radius: Union[str, Quantity] = Quantity(0, 'arcsec'),
-                     start_temp: Quantity = Quantity(3.0, "keV"), start_met: float = 0.3,
+def double_temp_apec(sources: Union[BaseSource, BaseSample], outer_radius: Union[str, Quantity],
+                     inner_radius: Union[str, Quantity] = Quantity(0, 'arcsec'),
+                     start_temp_one: Quantity = Quantity(3.0, "keV"), start_temp_two: Quantity = Quantity(5.0, "keV"),
+                     start_met_one: float = 0.3, start_met_two: float = 0.3,
                      lum_en: Quantity = Quantity([[0.5, 2.0], [0.01, 100.0]], "keV"), freeze_nh: bool = True,
-                     freeze_met: bool = True, freeze_temp: bool = False, lo_en: Quantity = Quantity(0.3, "keV"),
+                     freeze_met_one: bool = True, freeze_met_two: bool = True, freeze_temp_one: bool = False,
+                     freeze_temp_two: bool = False, lo_en: Quantity = Quantity(0.3, "keV"),
                      hi_en: Quantity = Quantity(7.9, "keV"), par_fit_stat: float = 1., lum_conf: float = 68.,
                      abund_table: str = "angr", fit_method: str = "leven", group_spec: bool = True, min_counts: int = 5,
                      min_sn: float = None, over_sample: float = None, one_rmf: bool = True, num_cores: int = NUM_CORES,
-                     spectrum_checking: bool = True, timeout: Quantity = Quantity(1, 'hr')):
+                     spectrum_checking: bool = False, timeout: Quantity = Quantity(1, 'hr')):
     """
-    This is a convenience function for fitting an absorbed multi-temperature apec model(constant*tbabs*apec) to an
-    object's spectrum - the number of temperature components (i.e. APEC models) is set using the 'num_temp_comp'
-    argument. If there are no existing spectra with the passed settings, then they will be generated automatically.
+    This is a convenience function for fitting an absorbed double-temperature apec model (constant*tbabs*(apec+apec))
+    to an object's spectrum. If there are no existing spectra with the passed settings, then they will be
+    generated automatically.
 
     BE AWARE that higher quality data (i.e. higher signal-to-noise spectra, with a greater number of X-ray counts) are
-    required to fit models with more free parameters, so the more temperature components you include, the less likely
-    the model fit is to be successful.
+    required to fit models with more free parameters, so this model fit is less likely to converge than single
+    temperature apec fits for lower quality data.
 
     :param List[BaseSource] sources: A single source object, or a sample of sources.
     :param str/Quantity outer_radius: The name or value of the outer radius of the region that the
@@ -965,14 +957,24 @@ def n_temp_apec(sources: Union[BaseSource, BaseSample], outer_radius: Union[str,
         desired spectrum covers (for instance 'r200' would be acceptable for a GalaxyCluster,
         or Quantity(1000, 'kpc')). By default this is zero arcseconds, resulting in a circular spectrum. If
         you are fitting for multiple sources then you can also pass a Quantity with one entry per source.
-    :param Quantity start_temp: The initial temperature for the fit, the default is 3 keV. This value can also be
-        a non-scalar Quantity, with an entry for every source in a sample (this is most useful when used with the
-        'freeze_temp' argument, to provide some external constraint on temperature for objects with poor data).
-    :param start_met: The initial metallicity for the fit (in ZSun).
+    :param Quantity start_temp_one: The initial temperature of the first APEC model, the default is 3 keV. This value
+        can also be a non-scalar Quantity, with an entry for every source in a sample (this is most useful when
+        used with the 'freeze_temp' argument, to provide some external constraint on temperature for objects with
+        poor data).
+    :param Quantity start_temp_two: The initial temperature of the second APEC model, the default is 3 keV. This value
+        can also be a non-scalar Quantity, with an entry for every source in a sample (this is most useful when
+        used with the 'freeze_temp' argument, to provide some external constraint on temperature for objects with
+        poor data).
+    :param start_met_one: The initial metallicity of the first APEC model (in ZSun).
+    :param start_met_two: The initial metallicity of the second APEC model (in ZSun).
     :param Quantity lum_en: Energy bands in which to measure luminosity.
     :param bool freeze_nh: Whether the hydrogen column density should be frozen. Default is True.
-    :param bool freeze_met: Whether the metallicity parameter in the fit should be frozen. Default is True.
-    :param bool freeze_temp: Whether the temperature parameter in the fit should be frozen. Default is False
+    :param bool freeze_met_one: Whether the metallicity of the first APEC model should be frozen. Default is True.
+    :param bool freeze_met_two: Whether the metallicity of the second APEC model should be frozen. Default is True.
+    :param bool freeze_temp_one: Whether the temperature parameter of the first APEC model should be
+        frozen. Default is False
+    :param bool freeze_temp_two: Whether the temperature parameter of the second APEC model should be
+    frozen. Default is False
     :param Quantity lo_en: The lower energy limit for the data to be fitted.
     :param Quantity hi_en: The upper energy limit for the data to be fitted.
     :param float par_fit_stat: The delta fit statistic for the XSPEC 'error' command, default is 1.0 which should be
@@ -985,7 +987,7 @@ def n_temp_apec(sources: Union[BaseSource, BaseSample], outer_radius: Union[str,
     :param float min_counts: If generating a grouped spectrum, this is the minimum number of counts per channel.
         To disable minimum counts set this parameter to None.
     :param float min_sn: If generating a grouped spectrum, this is the minimum signal to noise in each channel.
-        To disable minimum signal to noise set this parameter to None.
+        To disable minimum signal-to-noise set this parameter to None.
     :param float over_sample: The minimum energy resolution for each group, set to None to disable. e.g. if
         over_sample=3 then the minimum width of a group is 1/3 of the resolution FWHM at that energy.
     :param bool one_rmf: This flag tells the method whether it should only generate one RMF for a particular
@@ -998,37 +1000,56 @@ def n_temp_apec(sources: Union[BaseSource, BaseSample], outer_radius: Union[str,
         Please note that this is not a timeout for the entire fitting process, but a timeout to individual source
         fits.
     """
-    raise NotImplementedError("Barely begun")
+
+    # This is a little cheesy as we haven't decided what exactly to do yet, but the spectrum checking behaviour
+    #  is super inconsistent with double_temp_apec - presumably just because the greater number of parameters
+    #  makes estimating errors much less stable for individual spectral fits
+    if spectrum_checking:
+        raise NotImplementedError("The double_temp_apec function does not currently support "
+                                  "'spectrum_checking=True', as it exhibits very unstable behaviour.")
+
     sources, inn_rad_vals, out_rad_vals = _pregen_spectra(sources, outer_radius, inner_radius, group_spec, min_counts,
                                                           min_sn, over_sample, one_rmf, num_cores)
     sources = _check_inputs(sources, lum_en, lo_en, hi_en, fit_method, abund_table, timeout)
 
-    # Have to check that every source has a start temperature entry, if the user decided to pass a set of them
-    if not start_temp.isscalar and len(start_temp) != len(sources):
-        raise ValueError("If a non-scalar Quantity is passed for 'start_temp', it must have one entry for each "
-                         "source. It currently has {n} for {s} sources.".format(n=len(start_temp), s=len(sources)))
+    # Have to check that every source has start temperature entries, if the user decided to pass a set of them - this
+    #  is identical to the checks we perform at the beginning of 'single_temp_apec'
+    if ((not start_temp_one.isscalar and len(start_temp_one) != len(sources)) or
+            (not start_temp_two.isscalar and len(start_temp_two) != len(sources))):
+        raise ValueError("If a non-scalar Quantity is passed for 'start_temp_one' or 'start_temp_two', there must be "
+                         "one entry for each source. They currently have {n1} and {n2} respectively, for {s} "
+                         "sources.".format(n1=len(start_temp_one), n2=len(start_temp_two), s=len(sources)))
+
     # Want to make sure that the start_temp variable is always a non-scalar Quantity with an entry for every source
     #  after this point, it means we normalise how we deal with it.
-    elif start_temp.isscalar:
+    if start_temp_one.isscalar:
         # Doing it like this, defining a new variable and then redeclaring the start_temp in each bit of the loop
         #  below is important, as the fit_conf generation code accesses the current value of start_temp specifically
-        all_start_temps = Quantity([start_temp.value.copy()]*len(sources), start_temp.unit)
+        all_start_temp_ones = Quantity([start_temp_one.value.copy()]*len(sources), start_temp_one.unit)
     else:
-        all_start_temps = start_temp.copy()
+        all_start_temp_ones = start_temp_one.copy()
 
-    # This function is for a set model, absorbed apec, so I can hard code all of this stuff.
-    # These will be inserted into the general XSPEC script template, so lists of parameters need to be in the form
-    #  of TCL lists.
-    model = "constant*tbabs*apec"
-    par_names = "{factor nH kT Abundanc Redshift norm}"
+    # And the same thing for start temperature the second
+    if start_temp_two.isscalar:
+        # Doing it like this, defining a new variable and then redeclaring the start_temp in each bit of the loop
+        #  below is important, as the fit_conf generation code accesses the current value of start_temp specifically
+        all_start_temp_twos = Quantity([start_temp_two.value.copy()]*len(sources), start_temp_two.unit)
+    else:
+        all_start_temp_twos = start_temp_two.copy()
+
+    # This function is for a set model, absorbed double apec, so I can hard code all of this stuff.
+    model = "constant*tbabs*(apec+apec)"
+    # These will be inserted into the general XSPEC script template, so lists of parameters need to be
+    #  in the form of TCL lists.
+    par_names = "{factor nH kT Abundanc Redshift norm kT Abundanc Redshift norm}"
     lum_low_lims = "{" + " ".join(lum_en[:, 0].to("keV").value.astype(str)) + "}"
     lum_upp_lims = "{" + " ".join(lum_en[:, 1].to("keV").value.astype(str)) + "}"
 
     # Here we generate the fit configuration storage key from those arguments to this function that control the fit
     #  and how it behaves
-    rel_args = FIT_FUNC_ARGS['single_temp_apec']
+    rel_args = FIT_FUNC_ARGS['double_temp_apec']
 
-    sig = signature(single_temp_apec)
+    sig = signature(double_temp_apec)
     cur_args = {k: v.default for k, v in sig.parameters.items() if v.default is not Parameter.empty}
 
     # This is purely for developers, as a check to make sure that the FIT_FUNC_ARGS dictionary is updated if the
@@ -1063,22 +1084,27 @@ def n_temp_apec(sources: Union[BaseSource, BaseSample], outer_radius: Union[str,
         if source.redshift is None:
             raise ValueError("You cannot supply a source without a redshift to this model.")
 
-        # Whatever start temperature is passed gets converted to keV, this will be put in the template
-        start_temp = all_start_temps[src_ind].to("keV", equivalencies=u.temperature_energy())
+        # Whatever start temperatures are passed get converted to keV - then will be put in the template
+        start_temp_one = all_start_temp_ones[src_ind].to("keV", equivalencies=u.temperature_energy())
+        start_temp_two = all_start_temp_twos[src_ind].to("keV", equivalencies=u.temperature_energy())
+
         # Another TCL list, this time of the parameter start values for this model.
-        par_values = "{{{0} {1} {2} {3} {4} {5}}}".format(1., source.nH.to("10^22 cm^-2").value, start_temp.value,
-                                                          start_met, source.redshift, 1.)
+        par_values = ("{{{0} {1} {2} {3} {4} {5} {6} {7} {8} "
+                      "{9}}}").format(1., source.nH.to("10^22 cm^-2").value, start_temp_one.value, start_met_one,
+                                      source.redshift, 1., start_temp_two.value, start_met_two, source.redshift, 1.)
 
         # Set up the TCL list that defines which parameters are frozen, dependent on user input - this can now
         #  include the temperature, if the user wants it fixed at the start value
-        freezing = "{{F {n} {t} {a} T F}}".format(n="T" if freeze_nh else "F",
-                                                  t="T" if freeze_temp else "F",
-                                                  a="T" if freeze_met else "F")
+        freezing = "{{F {n} {t1} {a2} T F {t2} {a2} T F}}".format(n="T" if freeze_nh else "F",
+                                                                  t1="T" if freeze_temp_one else "F",
+                                                                  a1="T" if freeze_met_one else "F",
+                                                                  t2="T" if freeze_temp_two else "F",
+                                                                  a2="T" if freeze_met_two else "F")
 
         # Set up the TCL list that defines which parameters are linked across different spectra, only the
         #  multiplicative constant that accounts for variation in normalisation over different observations is not
         #  linked
-        linking = "{F T T T T T}"
+        linking = "{F T T T T T T T T T}"
 
         # If the user wants the spectrum cleaning step to be run, then we have to setup some acceptable
         #  limits. For this function they will be hardcoded, for simplicities sake, and we're only going to
@@ -1135,5 +1161,5 @@ def n_temp_apec(sources: Union[BaseSource, BaseSample], outer_radius: Union[str,
                         'constant*tbabs*zpowerlw': power_law,
                         'constant*tbabs*powerlaw': power_law,
                         'constant*tbabs*zbbody': blackbody,
-                        'constant*tbabs*bbody': blackbody
-                        }
+                        'constant*tbabs*bbody': blackbody,
+                        'constant*tbabs*(apec+apec)': double_temp_apec}

xga/xspec/fitconfgen.py

@@ -1,13 +1,13 @@
 #  This code is a part of X-ray: Generate and Analyse (XGA), a module designed for the XMM Cluster Survey (XCS).
-#  Last modified by David J Turner (turne540@msu.edu) 12/03/2025, 16:02. Copyright (c) The Contributors
+#  Last modified by David J Turner (turne540@msu.edu) 27/03/2025, 22:31. Copyright (c) The Contributors
 
 from inspect import signature, Parameter
 from types import FunctionType
 from typing import Union
 
 from astropy.units import Quantity
 
-# This constant very importantly defined whether each argument to the XGA XSPEC fitting functions should be included
+# This constant very importantly defines whether each argument to the XGA XSPEC fitting functions should be included
 #  in the fit configuration storage key - we define it here so that this information can be accessed outside the
 #  actual fit function. We will also check the arguments of each fit function against these entries to see whether
 #  an argument has been added that is not accounted for here.
@@ -18,6 +18,14 @@
                          'group_spec': False, 'min_counts': False, 'min_sn': False, 'over_sample': False,
                          'one_rmf': False, 'num_cores': False, 'spectrum_checking': False, 'timeout': False},
 
+    'double_temp_apec': {'inner_radius': False, 'start_temp_one': True, 'start_temp_two': True,
+                         'start_met_one': True, 'start_met_two': True, 'lum_en': False, 'freeze_nh': True,
+                         'freeze_met_one': True, 'freeze_met_two': True, 'freeze_temp_one': True,
+                         'freeze_temp_two': True, 'lo_en': True, 'hi_en': True, 'par_fit_stat': True,
+                         'lum_conf': False, 'abund_table': True, 'fit_method': True,
+                         'group_spec': False, 'min_counts': False, 'min_sn': False, 'over_sample': False,
+                         'one_rmf': False, 'num_cores': False, 'spectrum_checking': False, 'timeout': False},
+
     'single_temp_mekal': {'inner_radius': False, 'start_temp': True, 'start_met': True, 'lum_en': False,
                           'freeze_nh': True, 'freeze_met': True, 'freeze_temp': True, 'lo_en': True, 'hi_en': True,
                           'par_fit_stat': True, 'lum_conf': False, 'abund_table': True, 'fit_method': True,