Noise fix

arrakis/columns_possum.py

@@ -41,7 +41,12 @@
     ("dc_min_axis_err", float, "cat", "E_DC_Min", u.arcsec),
     ("dc_pa", float, "cat", "DC_PA", u.deg),
     ("dc_pa_err", float, "cat", "E_DC_PA", u.deg),
-    ("stokesI_err", float, "cat", "Noise", u.mJy / u.beam),
+    ("stokesI_err", float, "synth", "dIFullBand", u.Jy / u.beam),
+    ("stokesQ_err", float, "synth", "dQFullBand", u.Jy / u.beam),
+    ("stokesU_err", float, "synth", "dUFullBand", u.Jy / u.beam),
+    ("stokesI_bkg", float, "synth", "bIFullBand", u.Jy / u.beam),
+    ("stokesQ_bkg", float, "synth", "bQFullBand", u.Jy / u.beam),
+    ("stokesU_bkg", float, "synth", "bUFullBand", u.Jy / u.beam),
     ("beamdist", float, "cat", "Separation_Tile_Centre", u.deg),
     ("N_Gaus", int, "cat", "N_Gaus", u.dimensionless_unscaled),
     ("cat_id", str, "cat", "Gaussian_ID", None),
@@ -400,4 +405,16 @@
         "description": "Error in Stokes I model coefficients",
         "ucd": "stat.error;stat.fit.param;phys.polarization.stokes.I",
     },
+    "stokesI_bkg": {
+        "description": "Background level of Stokes I",
+        "ucd": "phot.flux.density;phys.polarization.stokes.I",
+    },
+    "stokesQ_bkg": {
+        "description": "Background level of Stokes Q",
+        "ucd": "phot.flux.density;phys.polarization.stokes.Q",
+    },
+    "stokesU_bkg": {
+        "description": "Background level of Stokes U",
+        "ucd": "phot.flux.density;phys.polarization.stokes.U",
+    },
 }

arrakis/rmsynth_oncuts.py

@@ -33,7 +33,7 @@
 from arrakis.logger import logger
 from arrakis.utils.database import get_db, test_db
 from arrakis.utils.fitsutils import getfreq
-from arrakis.utils.fitting import fit_pl
+from arrakis.utils.fitting import fit_pl, fitted_mean, fitted_std
 from arrakis.utils.io import try_mkdir
 from arrakis.utils.pipeline import logo_str
 
@@ -248,7 +248,10 @@ def cubelet_bane(cubelet: np.ndarray, header: fits.Header) -> Tuple[np.ndarray]:
         plane = plane[np.isfinite(plane)]
         if len(plane) == 0:
             continue
-        background[chan], noise[chan] = norm.fit(plane)
+        clipped_plane = sigma_clip(
+            plane, sigma=3, cenfunc=fitted_mean, stdfunc=fitted_std, maxiters=None
+        )
+        background[chan], noise[chan] = norm.fit(clipped_plane.compressed())
 
     return background, noise
 
@@ -559,14 +562,10 @@ def rmsynthoncut1d(
 
     data = [np.array(freq)]
     bkg_data = [np.array(freq)]
-    for stokes in "iqu":
-        if noStokesI and stokes == "i":
-            continue
+    for stokes in "qu" if noStokesI else "iqu":
         data.append(filtered_stokes_spectra.__getattribute__(stokes).data)
         bkg_data.append(filtered_stokes_spectra.__getattribute__(stokes).bkg)
-    for stokes in "iqu":
-        if noStokesI and stokes == "i":
-            continue
+    for stokes in "qu" if noStokesI else "iqu":
         data.append(filtered_stokes_spectra.__getattribute__(stokes).rms)
 
     # Run 1D RM-synthesis on the spectra
@@ -606,6 +605,16 @@ def rmsynthoncut1d(
             dst = os.path.join(plotdir, base)
             copyfile(src, dst)
 
+    # Update I, Q, U noise from data
+    for stokes in "qu" if noStokesI else "iqu":
+        rms = filtered_stokes_spectra.__getattribute__(stokes).rms
+        bkg = filtered_stokes_spectra.__getattribute__(stokes).bkg
+        mean_rms = np.nanmean(rms)
+        mean_bkg = np.nanmean(bkg)
+        full_rms = mean_rms / np.sqrt(len(rms[np.isfinite(rms)]))
+        mDict[f"d{stokes.capitalize()}"] = mean_rms
+        mDict[f"d{stokes.capitalize()}FullBand"] = full_rms
+        mDict[f"b{stokes.capitalize()}FullBand"] = mean_bkg
     # Update model values if own fit was used
     if do_own_fit:
         mDict = update_rmtools_dict(mDict, stokes_i_fit_result.fit_dict)
@@ -623,6 +632,9 @@ def rmsynthoncut1d(
         mDict["fit_flag_is_not_finite"] = mDict["IfitStat"] >= 16
         mDict["fit_flag_is_not_normal"] = mDict["IfitStat"] >= 5
 
+    logger.info(f"RM-Synthesis for {cname} complete")
+    logger.info(f"RM-Synthesis results for {cname}: {pformat(mDict)}")
+
     # Ensure JSON serializable
     for k, v in mDict.items():
         if isinstance(v, np.float_):
@@ -647,6 +659,7 @@ def rmsynthoncut1d(
     head_dict.pop("", None)
     if "COMMENT" in head_dict.keys():
         head_dict["COMMENT"] = str(head_dict["COMMENT"])
+    logger.debug(f"Heading for {cname} is {pformat(head_dict)}")
 
     outer_dir = os.path.basename(os.path.dirname(filtered_stokes_spectra.i.filename))
 

arrakis/utils/fitting.py

@@ -3,13 +3,13 @@
 
 import warnings
 from functools import partial
-from typing import Tuple
+from typing import Optional, Tuple
 
 import numpy as np
 from astropy.stats import akaike_info_criterion_lsq
 from astropy.utils.exceptions import AstropyWarning
 from scipy.optimize import curve_fit
-from scipy.stats import normaltest
+from scipy.stats import norm, normaltest
 from spectral_cube.utils import SpectralCubeWarning
 
 from arrakis.logger import logger
@@ -18,6 +18,36 @@
 warnings.simplefilter("ignore", category=AstropyWarning)
 
 
+def fitted_mean(data: np.ndarray, axis: Optional[int] = None) -> float:
+    """Calculate the mean of a distribution.
+
+    Args:
+        data (np.ndarray): Data array.
+
+    Returns:
+        float: Mean.
+    """
+    if axis is not None:
+        raise NotImplementedError("Axis not implemented")
+    mean, _ = norm.fit(data)
+    return mean
+
+
+def fitted_std(data: np.ndarray, axis: Optional[int] = None) -> float:
+    """Calculate the standard deviation of a distribution.
+
+    Args:
+        data (np.ndarray): Data array.
+
+    Returns:
+        float: Standard deviation.
+    """
+    if axis is not None:
+        raise NotImplementedError("Axis not implemented")
+    _, std = norm.fit(data)
+    return std
+
+
 def chi_squared(model: np.ndarray, data: np.ndarray, error: np.ndarray) -> float:
     """Calculate chi squared.