added doublet ratio criteria for rejecting false positives

qsoabsfind/absfinder.py

@@ -11,7 +11,7 @@
     median_selection_after_combining,
     remove_Mg_falsely_come_from_Fe_absorber, z_abs_from_same_metal_absorber,
     contiguous_pixel_remover, estimate_snr_for_lines, absorber_search_window,
-    find_valid_indices
+    find_valid_indices, calculate_doublet_ratio
 )
 from .ew import measure_absorber_properties_double_gaussian
 from .config import load_constants
@@ -226,8 +226,12 @@ def convolution_method_absorber_finder_in_QSO_spectra(spec_index, absorber='MgII
                         sn1, sn2 = estimate_snr_for_lines(c0, c1, lam_rest, residual, error, logwave)
                         # resolution corrected velocity dispersion (should be greater than 0)
                         vel1, vel2 = vel_dispersion(c0, c1, gaussian_parameters[2], gaussian_parameters[5], resolution)
+                        # calculate best -fit doublet ratio and errors and check if they are within the range.
+                        # usually 1 < DR < f1/f2 (doublet ratio =2, for MgII, CIV)
+                        dr, dr_error = calculate_doublet_ratio(EW_first_temp_mean[0], EW_second_temp_mean[0], EW_first_error_temp[0], EW_second_error_temp[0])
+                        min_dr, max_dr = 1 - 2 * dr_error, f1/f2 + 2 * dr_error
 
-                        if (gaussian_parameters > bound[0] + 0.01).all() and (gaussian_parameters < bound[1] - 0.01).all() and lower_del_lam <= c1 - c0 <= upper_del_lam and sn1 > sn_line1 and sn2 > sn_line2 and vel1 > 0 and vel2 > 0:
+                        if (gaussian_parameters > bound[0] + 0.01).all() and (gaussian_parameters < bound[1] - 0.01).all() and lower_del_lam <= c1 - c0 <= upper_del_lam and sn1 > sn_line1 and sn2 > sn_line2 and vel1 > 0 and vel2 > 0 and min_dr <= dr <=max_dr:
                             pure_z_abs[m] = z_new
                             pure_gauss_fit[m] = fit_param_temp[0]
                             pure_gauss_fit_std[m] = fit_param_std_temp[0]

qsoabsfind/absorberutils.py

@@ -88,6 +88,33 @@ def estimate_local_sigma_conv_array(conv_array, pm_pixel):
 
     return sigma_cr
 
+@jit(nopython=True)
+def calculate_doublet_ratio(ew1, ew2, ew1_error, ew2_error):
+    """
+    Calculate the doublet ratio and its associated error.
+
+    Args:
+        ew1 (float): Equivalent width of the first line.
+        ew2 (float): Equivalent width of the second line.
+        ew1_error (float): Error associated with the first equivalent width.
+        ew2_error (float): Error associated with the second equivalent width.
+
+    Returns:
+        tuple: A tuple containing:
+            - doublet_ratio (float): The ratio of ew1 to ew2.
+            - doublet_ratio_error (float): The propagated error for the doublet ratio.
+    """
+    if ew2 == 0:
+        raise ValueError("ew2 cannot be zero, as it would result in division by zero.")
+
+    # Calculate the doublet ratio
+    doublet_ratio = ew1 / ew2
+
+    # Propagate the error using standard error propagation formula for division
+    doublet_ratio_error = doublet_ratio * np.sqrt((ew1_error / ew1) ** 2 + (ew2_error / ew2) ** 2)
+
+    return doublet_ratio, doublet_ratio_error
+
 @jit(nopython=True)
 def estimate_snr_for_lines(l1, l2, lam_rest, residual, error, log):
     """