add telescope lim_mag, update test

frb/scripts/limiting_mag.py

@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 """
-Estimate the limiting luminosity given the magnitude limit and DM and coord
+Estimate the limiting luminosity given the magnitude limit and DMs and coord and telescope
 """
 from IPython import embed
 
@@ -12,8 +12,13 @@ def parser(options=None):
     parser.add_argument("DM_FRB", type=float, help="FRB DM")
     parser.add_argument("mag_limit", type=float, help="Magnitude limit in filter *without* extinction correction")
     parser.add_argument("--filter", type=str, default='DECaL_r', help="Filter -- only used for extinction correction.  Must be a Repo approved choice")
-    parser.add_argument("--dm_hostmw", type=float, default=100., help="Assumed DM contribution from MW and Host")
+    parser.add_argument("--dm_host", type=float, default=50., help="Assumed DM contribution from the Host. Default = 50")
+    parser.add_argument("--dm_mwhalo", type=float, default=50., help="Assumed DM contribution from the MW halo. Default = 50")
     #parser.add_argument("-v", "--verbose", default=False, action="store_true", help="Overwhelm the screen?")
+    parser.add_argument("--telescope", type=str, default='perfect', help="telescope model for the DM-z grid: CHIME, DSA, Parkes, FAST, CRAFT, \
+                        CRAFT_ICS_892/1300/1632, perfect. Default = perfect")
+    parser.add_argument("--cl", type=tuple, default=(2.5,97.5), 
+                        help="Confidence limits for the z estimate [default is a 95 percent c.l., (2.5,97.5)]")
 
     if options is None:
         pargs = parser.parse_args()
@@ -50,23 +55,58 @@ def main(pargs):
     DM_ISM = mw.ismDM(icoord)
     print(f"NE2001 = {DM_ISM}")
 
-    # DM Cosmic
-    DM_cosmic = pargs.DM_FRB - DM_ISM.value - pargs.dm_hostmw
+    # DM cosmic and EG
+    DM_extragalactic = pargs.DM_FRB - DM_ISM.value - pargs.dm_mwhalo
+    DM_cosmic = DM_extragalactic - pargs.dm_host
 
     # Redshift estimates
 
-    # Load
-    sdict = prob_dmz.grab_repo_grid()
+    # Load the telescope specific grid
+    telescope_dict = {
+        'CHIME': 'CHIME_pzdm.npz',
+        'DSA': 'DSA_pzdm.npy',
+        'Parkes': 'parkes_mb_class_I_and_II_pzdm.npy',
+        'CRAFT': 'CRAFT_class_I_and_II_pzdm.npy',
+        'CRAFT_ICS_1300': 'CRAFT_ICS_1300_pzdm.npy',
+        'CRAFT_ICS_892': 'CRAFT_ICS_892_pzdm.npy',
+        'CRAFT_ICS_1632': 'CRAFT_ICS_1632_pzdm.npy',
+        'FAST': 'FAST_pzdm.npy',
+        'perfect': 'PDM_z.npz'
+    }
+
+    # Get the perfect telescope grid (default)
+    sdict = prob_dmz.grab_repo_grid(telescope_dict['perfect'])
     PDM_z = sdict['PDM_z']
     z = sdict['z']
     DM = sdict['DM']
 
-    # Do it
+    # Grab the right entry
     iDM = np.argmin(np.abs(DM - DM_cosmic))
     PzDM = PDM_z[iDM, :] / np.sum(PDM_z[iDM, :])
 
+
+    # Get the telescope specific PZDM grid
+    if pargs.telescope and pargs.telescope != 'CHIME' and pargs.telescope != 'perfect':
+        if pargs.telescope not in telescope_dict:
+            raise ValueError(f"Unknown telescope: {pargs.telescope}")
+        zdict = prob_dmz.grab_repo_grid(telescope_dict['CHIME'])
+        z = zdict['z']
+        DM = zdict['DM']
+        PDM_z = prob_dmz.grab_repo_grid(telescope_dict[pargs.telescope])
+        iDM = np.argmin(np.abs(DM - DM_extragalactic))
+        PzDM = PDM_z[:,iDM] / np.sum(PDM_z[:,iDM])
+
+
+    if pargs.telescope and pargs.telescope == 'CHIME':
+        sdict = prob_dmz.grab_repo_grid(telescope_dict['CHIME'])
+        PDM_z = sdict['pzdm']
+        z = sdict['z']
+        DM = sdict['DM']
+        iDM = np.argmin(np.abs(DM - DM_extragalactic))
+        PzDM = PDM_z[:,iDM] / np.sum(PDM_z[:,iDM])
+
     cum_sum = np.cumsum(PzDM)
-    limits = [10, 90]
+    limits = pargs.cl
 
     z_min = z[np.argmin(np.abs(cum_sum-limits[0]/100.))]
     z_max = z[np.argmin(np.abs(cum_sum-limits[1]/100.))]
@@ -92,8 +132,18 @@ def main(pargs):
 
     # Finish
     print("-----------------------------------------------------")
-    print(f"For z_{limits[0]}={z_min:.2f}, the limiting magnitude corresponds to L={frac_Lstar_min:.5f}L*")
-    print(f"For z_{limits[1]}={z_max:.2f}, the limiting magnitude corresponds to L={frac_Lstar_max:.5f}L*")
+    print("")
+    print(f"Allowing for the MW halo, DM_MW_halo = {int(pargs.dm_mwhalo)} pc/cm^3")
+    print(f"Allowing for the Host, DM_host = {int(pargs.dm_host)} pc/cm^3")
+    print("")
+    if not pargs.telescope or pargs.telescope == 'perfect':
+        print("WARNING: This all assumes a perfect telescope and a model of the scatter in DM_cosmic (Macquart+2020)")
+    else:
+        print("This assumes the "+(str(pargs.telescope))+" telescope and a model of the scatter in DM_cosmic (Macquart+2020)")
+    print("")
+
+    print(f"For z_({limits[0]} %)={z_min:.2f}, the limiting magnitude corresponds to L={frac_Lstar_min:.5f}L*")
+    print(f"For z_({limits[1]} %)={z_max:.2f}, the limiting magnitude corresponds to L={frac_Lstar_max:.5f}L*")
 
     return frac_Lstar_min, frac_Lstar_max
 

frb/scripts/pz_dm.py

@@ -124,9 +124,9 @@ def main(pargs):
     print(f"z = [{z_min:.3f}, {z_max:.3f}]")
     print("")
     if not pargs.telescope or pargs.telescope == 'perfect':
-        print("WARNING: This all assumes a perfect telescope and a model of the scatter in DM_cosmic (Macqurt+2020)")
+        print("WARNING: This all assumes a perfect telescope and a model of the scatter in DM_cosmic (Macquart+2020)")
     else:
-        print("This assumes the "+(str(pargs.telescope))+" telescope and a model of the scatter in DM_cosmic (Macqurt+2020)")
+        print("This assumes the "+(str(pargs.telescope))+" telescope and a model of the scatter in DM_cosmic (Macquart+2020)")
     print("-----------------------------------------------------")
 
 

frb/tests/test_scripts.py

@@ -22,7 +22,7 @@ def test_frb_mag_limit():
     pargs = limiting_mag.parser(['J151849.52+122235.8', '200.', '23.'])
     Lmin, Lmax = limiting_mag.main(pargs)
 
-    assert np.isclose(Lmax, 0.018052542432481264)
+    assert np.isclose(Lmax, 0.023780306538345886)
 
 def test_frb_pz_dm():
     # Requires a file on disk that is too slow to generate in CI