Merge pull request #108 from sagasurvey/sfr_error

Proper SFR error propagation

SAGA/objects/build3.py

@@ -601,7 +601,7 @@ def add_sfr(base):
         base["Halpha_sfr_err"] = np.float32(np.nan)
 
         mask = Query(C.valid_Halpha, "dist_estimate > 0").mask(base)
-        t = base[["EW_Halpha", "EW_Halpha_err", "SPEC_Z", "Mr"]][mask]
+        t = base[["EW_Halpha", "EW_Halpha_err", "SPEC_Z", "SPEC_Z_ERR", "Mr", "r_err"]][mask]
         with np.errstate(all="ignore"):
             log_sfr, log_sfr_err = calc_sfr.calc_SFR_Halpha(*(c.astype(np.float64) for c in t.itercols()))
 

SAGA/objects/calc_sfr.py

@@ -9,7 +9,7 @@
 _IMF_FACTOR = 0.66
 
 
-def calc_SFR_NUV(NUV_mag, NUV_mag_err, dist_mpc, internal_ext=0.9):
+def calc_SFR_NUV(NUV_mag, NUV_mag_err, dist_mpc, internal_ext=0.9, internal_ext_err=0.2):
     """
     Convert NUV magnitudes into a SFR
     Based on Iglesias-Paramo (2006), Eq 3
@@ -26,7 +26,6 @@ def calc_SFR_NUV(NUV_mag, NUV_mag_err, dist_mpc, internal_ext=0.9):
     # CONVERT GALEX m_AB TO FLUX:  erg sec-1 cm-2 Angstrom-1)
     # https://asd.gsfc.nasa.gov/archive/galex/FAQ/counts_background.html
     log_flux_nuv = -0.4 * (m_nuv_ab - 20.08 - 2.5 * np.log10(2.06e-16))
-    log_flux_nuv_err = 0.4 * (NUV_mag_err)  # ADD EXTRA ERROR FOR REDDENING OR DISTANCE?
 
     # LUMINOSITY (erg/s/A-1)
     # 796A is NUV filter width
@@ -37,20 +36,22 @@ def calc_SFR_NUV(NUV_mag, NUV_mag_err, dist_mpc, internal_ext=0.9):
 
     # CONVVERT TO SFR: EQ 3, inglesias- paramo 2006, also account for Salpeter -> Koupa IMF
     log_SFR_NUV = l_nuv_msun - 9.33 + np.log10(_IMF_FACTOR)
-    log_SFR_NUV_err = log_flux_nuv_err
+
+    # PROPAGATE ERRORS: assume ANUV_ERR and measurement errors
+    # ANUV_ERR is determined to be consistent with BD_ERR
+    log_SFR_NUV_err = 0.4 * np.hypot(NUV_mag_err, internal_ext_err)
 
     return log_SFR_NUV, log_SFR_NUV_err
 
 
-def calc_SFR_Halpha(EW_Halpha, EW_Halpha_err, spec_z, Mr, EWc=2.5, BD=3.25):
+def calc_SFR_Halpha(EW_Halpha, EW_Halpha_err, spec_z, spec_z_err, Mr, r_err, EWc=2.5, BD=3.25, BD_err=0.1):
     """
     Calculate Halpha-based EW SFR
     Bauer+ (2013) https://ui.adsabs.harvard.edu/abs/2013MNRAS.434..209B/abstract
     """
 
     # Bauer, EQ 2, term1
     term1 = (EW_Halpha + EWc) * 10 ** (-0.4 * (Mr - 34.1))
-    term1_err = EW_Halpha_err * 10 ** (-0.4 * (Mr - 34.1))
 
     # Bauer Eq 2, term2
     term2 = 3e18 / (6564.6 * (1.0 + spec_z)) ** 2
@@ -59,16 +60,18 @@ def calc_SFR_Halpha(EW_Halpha, EW_Halpha_err, spec_z, Mr, EWc=2.5, BD=3.25):
     term3 = (BD / 2.86) ** 2.36
 
     L_Halpha = term1 * term2 * term3
-    L_Halpha_err = term1_err * term2 * term3
 
     # EQ 3, Bauer et al above, also account for Salpeter -> Koupa IMF
     SFR = (L_Halpha * _IMF_FACTOR) / 1.27e34
-    SFR_err = (L_Halpha_err * _IMF_FACTOR) / 1.27e34
-
     log_Ha_SFR = np.log10(SFR)
 
-    # PROPOGATE ERRORS
-    log_SFR_err2 = SFR_err**2 * (1.0 / (SFR * np.log(10.0))) ** 2
-    log_Ha_SFR_err = np.sqrt(log_SFR_err2)
+    # PROPAGATE ERRORS: EW_err, Mr_err and AV_err
+    term1_EW_frac_err = EW_Halpha_err / (EW_Halpha + EWc)
+    term1_Mr_frac_err = 0.4 * np.log(10) * r_err
+    term1_frac_err = np.hypot(term1_EW_frac_err, term1_Mr_frac_err)
+    term2_frac_err = 2.0 * spec_z_err / (1.0 + spec_z)
+    term3_frac_err = 2.36 * (BD_err / BD)
+    L_Halpha_frac_err = np.sqrt(term1_frac_err ** 2 + term2_frac_err ** 2 + term3_frac_err ** 2)
+    log_Ha_SFR_err  = L_Halpha_frac_err / np.log(10)
 
     return log_Ha_SFR, log_Ha_SFR_err

SAGA/version.py

@@ -1,4 +1,5 @@
 """
 SAGA package version
 """
-__version__ = "0.71.1"
+
+__version__ = "0.72.0"