Enhancements for EM Counterpart Detection and Bandpass Interpolation (sncosmo)

nmma/em/create_lightcurves.py

@@ -203,9 +203,9 @@ def get_parser():
     parser.add_argument(
         "--ylim",
         type=str,
-        default="22,16",
+        default="-12,-18",
         nargs="*",
-        help="Upper and lower magnitude limit for light curve plot (default: 22-16)",
+        help="Upper and lower magnitude limit for light curve plot (default: -12,-18)",
     )
     parser.add_argument(
         "--photometric-error-budget",
@@ -400,7 +400,7 @@ def main(args=None):
             "ytick.labelsize": 42,
             "text.usetex": True,
             "font.family": "Times New Roman",
-            "figure.figsize": [16, 20],
+            "figure.figsize": [18, 25],
         }
         matplotlib.rcParams.update(params)
 
@@ -410,8 +410,9 @@ def main(args=None):
 
         fig = plt.figure()
 
-        filts = list(set(mag_ds[index].keys()).difference({"t"}))
-
+        #filts = list(set(mag_ds[index].keys()).difference({"t"}))
+        filts = filters
+
         ncols = 1
         nrows = int(np.ceil(len(filts) / ncols))
         gs = fig.add_gridspec(nrows=nrows, ncols=ncols, wspace=0.6, hspace=0.5)
@@ -466,29 +467,30 @@ def return_hist(x):
             ylim = [float(x) for x in ylim]
             plt.ylim(ylim)
 
-            ax.set_ylabel(filt, fontsize=30, rotation=0, labelpad=14)
+            #ax.set_ylabel(filt, fontsize=30, rotation=0, labelpad=14)
+            ax.set_ylabel(filt, fontsize=30, rotation=90, labelpad=8)
 
             if ii == len(filts) - 1:
                 ax.set_xticks([np.ceil(x) for x in np.linspace(xlim[0], xlim[1], 8)])
             else:
                 plt.setp(ax.get_xticklabels(), visible=False)
 
             ax.set_yticks([np.ceil(x) for x in np.linspace(ylim[0], ylim[1], 8)])
-            ax.tick_params(axis="x", labelsize=42)
-            ax.tick_params(axis="y", labelsize=42)
+            ax.tick_params(axis="x", labelsize=35)
+            ax.tick_params(axis="y", labelsize=35)
             ax.grid(which="both", alpha=0.5)
 
             axs.append(ax)
 
         fig.colorbar(c, ax=axs, location="right", shrink=0.6)
         fig.text(0.4, 0.05, r"Time [days]", fontsize=42)
         fig.text(
-            0.01,
+            0.001,
             0.5,
             r"Absolute Magnitude",
             va="center",
             rotation="vertical",
-            fontsize=42,
+            fontsize=35,
         )
 
         # plt.tight_layout()

nmma/em/detect_lightcurves.py

@@ -61,12 +61,6 @@ def main():
         required=True, 
         help="Either BNS or NSBH"
     )
-    parser.add_argument(
-        "-c", 
-        "--configDirectory",
-        help="gwemopt config file directory.", 
-        required=True
-    )
     parser.add_argument(
         "--outdir", 
         type=str, 
@@ -128,7 +122,12 @@ def main():
         default=1, 
         help="Number of cores"
     )
-
+    parser.add_argument(
+        "--label", 
+        type=str, 
+        required=True, 
+        help="Label for the run"
+    )
     args = parser.parse_args()
 
     # load the injection json file
@@ -159,7 +158,7 @@ def main():
             os.makedirs(outdir)
 
         skymap_file = os.path.join(args.skymap_dir, "%d.fits" % indices[index])
-        lc_file = os.path.join(args.lightcurve_dir, "%d.dat" %  indices[index])
+        lc_file = os.path.join(args.lightcurve_dir, "%s.dat" %  f"{args.label}_{indices[index]}")
 
         # fixed scheduling time as observation plan
         #number_shot = int(1 + (args.tmax - args.tmin)/args.dt)
@@ -168,7 +167,7 @@ def main():
 
         lcs[index] = np.loadtxt(lc_file)
 
-        efffile = os.path.join(outdir, f"efficiency_true_{indices[index]}.txt")
+        efffile = os.path.join(outdir, f"efficiency_true_{args.label}_{indices[index]}.txt")
         if os.path.isfile(efffile):
             continue
 
@@ -186,9 +185,14 @@ def main():
             [str(args.exposuretime) for i in args.filters.split(",")]
         )
 
-        command = f"gwemopt_run --telescopes {args.telescope} --do3D --doTiles --doSchedule --doSkymap --doTrueLocation --true_ra {ra} --true_dec {dec} --true_distance {dist} --doObservability --doObservabilityExit --timeallocationType powerlaw --scheduleType greedy -o {outdir} --gpstime {gpstime} --skymap {skymap_file} --filters {args.filters} --exposuretimes {exposuretime} --doSingleExposure --doAlternatingFilters --doEfficiency --lightcurveFiles {lc_file} --modelType file --configDirectory {args.configDirectory} --doPlots"
+        if args.telescope.lower() == 'ultrasat': 
+            command = f"gwemopt-run --telescopes {args.telescope} --geometry 3d --doTiles --doSchedule  --scheduleType greedy_slew --doTrueLocation --true_location  --true_ra {ra} --true_dec {dec} --true_distance {dist} --powerlaw_cl 0.9 --doObservability --doObservabilityExit --timeallocationType powerlaw   --gpstime {gpstime} --event {skymap_file} --filters {args.filters} --exposuretimes {exposuretime} --doSingleExposure --doAlternatingFilters --doEfficiency --lightcurveFiles {lc_file}  -o {outdir} --modelType file --doPlots --doUsePrimary  --ignore_observability"
+
+        else:
+            command = f"gwemopt-run --telescopes {args.telescope} --geometry 3d --doTiles --doSchedule  --scheduleType greedy_slew --doTrueLocation --true_location  --true_ra {ra} --true_dec {dec} --true_distance {dist} --powerlaw_cl 0.9 --doObservability --doObservabilityExit --timeallocationType powerlaw   --gpstime {gpstime} --event {skymap_file} --filters {args.filters} --exposuretimes {exposuretime} --doSingleExposure --doAlternatingFilters --doEfficiency --lightcurveFiles {lc_file}  -o {outdir} --modelType file --doPlots --doUsePrimary"
+
         commands.append(command)
-
+        
     print("Number of jobs remaining... %d." % len(commands))
 
     if args.parallel:
@@ -207,16 +211,22 @@ def main():
     for index, row in dataframe_from_inj.iterrows():
         outdir = os.path.join(args.outdir, str(index))
 
-        efffile = os.path.join(outdir, f"efficiency_true_{indices[index]}.txt")
+        efffile = os.path.join(outdir, f"efficiency_true_{args.label}_{indices[index]}.txt")
 
         ## Choose band for the best proxy
         # for ZTF r-band give the best proxy 
         if args.telescope == 'ZTF':
-            absmag.append(np.min(lcs[index][:, 3]))
+            absmag.append(np.min(lcs[index][:, 2]))
+
         # For Rubin i-band or z-band are best 
-        else:
-            absmag.append(np.min(lcs[index][:, 4]))
+        elif args.telescope == 'LSST':
+            absmag.append(np.min(lcs[index][:, 3]))
 
+        elif args.telescope.lower() == 'ultrasat':
+            absmag.append(np.min(lcs[index][:, 1]))
+
+        else:
+            absmag.append(np.min(lcs[index][:, 3]))
 
         if not os.path.isfile(efffile):
             fid.write("0\n")
@@ -233,8 +243,8 @@ def main():
             with open(efffile, "r") as file:
                 data_out = file.read()
             effs.append(float(data_out.split("\n")[1].split("\t")[4]))
-
-        efffile = os.path.join(outdir, f"efficiency_{indices[index]}.txt")
+            
+        efffile = os.path.join(outdir, f"efficiency_{args.label}_{indices[index]}.txt")
         if not os.path.isfile(efffile):
             probs.append(0.0)
         else:
@@ -372,7 +382,7 @@ def main():
 
     plt.axes(ax4)
     plt.axis("off")
-    cbar = plt.colorbar(sm, shrink=0.5, orientation="horizontal")
+    cbar = fig.colorbar(sm, ax=ax4, orientation='horizontal')
     cbar.set_label(r"Detection Efficiency")
 
     plt.axes(ax3)

nmma/em/utils.py

@@ -8,7 +8,7 @@
 import scipy.signal
 import scipy.constants
 import scipy.stats
-from sncosmo.bandpasses import _BANDPASSES
+from sncosmo.bandpasses import _BANDPASSES, _BANDPASS_INTERPOLATORS
 
 import sncosmo
 import dust_extinction.shapes as dustShp
@@ -124,6 +124,20 @@ def getRedShift(parameters):
             z = 0.0
     return z
 
+def get_all_bandpass_metadata():
+    """
+    Retrieves and combines the metadata for all registered bandpasses and interpolators.
+
+    Returns:
+        list: Combined list of metadata dictionaries from bandpasses and interpolators for sncosmo.
+    """
+
+    bandpass_metadata = _BANDPASSES.get_loaders_metadata()
+    interpolator_metadata = _BANDPASS_INTERPOLATORS.get_loaders_metadata()
+
+    combined_metadata = bandpass_metadata + interpolator_metadata
+
+    return combined_metadata
 
 def getFilteredMag(mag, filt):
     unprocessed_filt = [
@@ -151,7 +165,7 @@ def getFilteredMag(mag, filt):
         "swope2__J",
         "swope2__H",
     ]
-    sncosmo_filts = [val["name"] for val in _BANDPASSES.get_loaders_metadata()]
+    sncosmo_filts = [val["name"] for val in get_all_bandpass_metadata()]
     sncosmo_maps = {
         name.replace(":", "_"): name.replace(":", "_") for name in sncosmo_filts
     }
@@ -275,10 +289,17 @@ def get_default_filts_lambdas(filters=None):
         [lambdas_sloan, lambdas_bessel, lambdas_radio, lambdas_Xray]
     )
 
-    bandpasses = [
-        sncosmo.get_bandpass(val["name"]) for val in _BANDPASSES.get_loaders_metadata()
-    ]
-
+    bandpasses = []
+    for val in get_all_bandpass_metadata():
+        try:
+            bandpass = sncosmo.get_bandpass(val["name"])
+        except Exception :
+            if val["name"] == "ultrasat":
+                bandpass = sncosmo.get_bandpass(val["name"], 3)
+                bandpass.name = bandpass.name.split()[0]
+
+        bandpasses.append(bandpass)
+
     filts = filts + [band.name for band in bandpasses]
     lambdas = np.concatenate([lambdas, [1e-10 * band.wave_eff for band in bandpasses]])
 

nmma/eos/create_injection.py

@@ -58,20 +58,23 @@ def file_to_dataframe(
         injection_values["mass_2"].append(min(float(row["mass1"]), float(row["mass2"])))
         injection_values["luminosity_distance"].append(float(row["distance"]))
 
-        if "polarization" in row:
+        if "polarization" in row.colnames:
             injection_values["psi"].append(float(row["polarization"]))
         else:
             injection_values["psi"].append(0.0)
 
-        if "coa_phase" in row:
+        if "coa_phase" in row.colnames:
             coa_phase = float(row["coa_phase"])
             injection_values["phase"].append(float(row["coa_phase"]))
         else:
             coa_phase = 0.0
             injection_values["phase"].append(0.0)
 
-        if "geocent_end_time" in row:
-            injection_values["geocent_time"].append(float(row["geocent_end_time"]))
+        if "geocent_end_time" in row.colnames:
+            if "geocent_end_time_ns" in row.colnames:
+                injection_values["geocent_time"].append(float(row["geocent_end_time"]) + float(row["geocent_end_time_ns"]) * (10 ** -9))
+            else:
+                injection_values["geocent_time"].append(float(row["geocent_end_time"]))
         else:
             injection_values["geocent_time"].append(trigger_time)