Look for the closest declination line and adequate NSB random forest model for DL2 production

src/osa/paths.py

@@ -7,9 +7,11 @@
 from typing import List
 import subprocess
 import time
+import os
 
 import lstchain
 from astropy.table import Table
+from astropy.coordinates import SkyCoord
 from lstchain.onsite import (find_systematics_correction_file,
                              find_time_calibration_file,
                              find_filter_wheels)
@@ -397,10 +399,100 @@ def create_longterm_symlink(cherenkov_job_id: str = None):
     else:
         log.warning(f"Job {cherenkov_job_id} (lstchain_cherenkov_transparency) did not finish successfully.")
 
+
 def dl1_datacheck_longterm_file_exits() -> bool:
     """Return true if the longterm DL1 datacheck file was already produced."""
     nightdir = utils.date_to_dir(options.date)
     longterm_dir = Path(cfg.get("LST1", "LONGTERM_DIR"))
     longterm_file = longterm_dir / options.prod_id / nightdir / f"DL1_datacheck_{nightdir}.h5"
     return longterm_file.exists()
 
+
+def convert_dec_string(dec_str: str) -> float:
+    """Return the declination angle in degrees corresponding to a 
+    given string of the form "dec_XXXX" or "dec_min_XXXX"."""
+    # Split the string into parts
+    parts = dec_str.split('_')
+
+    # Extract the sign, degrees, and minutes
+    sign = 1 if 'min' not in parts else -1
+    degrees = int(parts[-1])
+
+    # Calculate the numerical value
+    dec_value = sign * (degrees / 100)
+
+    return dec_value
+
+
+def get_corresponding_string(list1: list, list2: list) -> dict:
+    """Return a dictionary created from two given lists."""
+    corresponding_dict = {}
+    for index, element in enumerate(list2):
+        corresponding_dict[element] = list1[index]
+    return corresponding_dict
+
+
+def get_latest_RF_model_path(dec_str: str) -> Path:
+    """Get the path of the most recent version of RF models for a given declination,
+    excluding the ones produced for the source-dependent analysis."""
+    BASE_MODELS = Path("/fefs/aswg/data/models/AllSky")
+    # make sure the RF models correspond to the current version of lstchain
+    current_version = get_major_version(utils.get_lstchain_version())
+    list_nodes = sorted(BASE_MODELS.rglob(f"*{current_version}*/{dec_str}"), key=os.path.getmtime)
+
+    log.debug(f"Found len(list_nodes) paths with {current_version} corresponding to {dec_str}:")
+    for path in list_nodes:
+        log.debug(path)
+
+    # remove from the list the models produced for the source-dependent analysis        
+    for i in list_nodes:
+        if "srcdep" in str(i):
+            list_nodes.remove(i)
+
+    return list_nodes[-1]
+
+
+def get_RF_model(run_str: str) -> Path:
+    """Get the path of the RF model to be used in the DL2 production for a given run."""
+    run_catalog_dir = Path("/fefs/aswg/data/real/monitoring/RunCatalog")
+    run_catalog_file = run_catalog_dir / f"RunCatalog_{options.date}.ecsv"
+    run_catalog = Table.read(run_catalog_file)
+    run = run_catalog[run_catalog["run_id"]==int(run_str)]
+
+    try:
+        target_name = run["source_name"]
+        source_coordinates = SkyCoord.from_name(target_name)
+        source_dec = source_coordinates.dec.value
+
+    except TypeError:
+        source_dec = run["source_dec"][0]
+
+    source_culmination = utils.culmination_angle(source_dec)
+
+    dec_list = os.listdir("/fefs/aswg/data/mc/DL0/LSTProd2/TrainingDataset/Protons")[:-2]
+
+    # Convert each string in the list to numerical values
+    dec_values = [convert_dec_string(dec) for dec in dec_list]
+
+    closest_declination = min(dec_values, key=lambda x: abs(x - source_dec))
+    closest_dec_culmination = utils.culmination_angle(closest_declination)
+    log.debug(
+        f"The declination closest to {source_dec} is: {closest_declination}."
+        "Checking if the culmination angle is larger than the one of the target source."
+    )
+
+    if closest_dec_culmination > source_culmination:
+        # If the culmination angle of the closest declination line is larger than for the source, 
+        # remove it from the declination lines list and look for the second closest declination line.
+        corresponding_dict = get_corresponding_string(dec_list, dec_values)
+        corresponding_string = corresponding_dict[closest_declination]
+        dec_values.remove(closest_declination)
+        dec_list.remove(corresponding_string)
+        closest_declination = min(dec_values, key=lambda x: abs(x - source_dec))
+
+        log.debug(f"The declination line to use for the DL2 production is: {closest_declination}")
+
+    corresponding_dict = get_corresponding_string(dec_list, dec_values)
+    corresponding_string = corresponding_dict[closest_declination]
+
+    return get_latest_RF_model_path(corresponding_string)

src/osa/scripts/datasequence.py

@@ -12,6 +12,7 @@
 from osa.utils.cliopts import data_sequence_cli_parsing
 from osa.utils.logging import myLogger
 from osa.utils.utils import date_to_dir
+from osa.paths import get_RF_model
 
 __all__ = ["data_sequence", "r0_to_dl1", "dl1_to_dl2", "dl1ab", "dl1_datacheck"]
 
@@ -265,7 +266,7 @@ def dl1_to_dl2(run_str: str) -> int:
     dl1ab_subdirectory = Path(options.directory) / options.dl1_prod_id
     dl2_subdirectory = Path(options.directory) / options.dl2_prod_id
     dl2_config = Path(cfg.get("lstchain", "dl2_config"))
-    rf_models_directory = Path(cfg.get("lstchain", "RF_MODELS"))
+    rf_models_directory = get_RF_model(run_str)
     dl1_file = dl1ab_subdirectory / f"dl1_LST-1.Run{run_str}.h5"
 
     command = cfg.get("lstchain", "dl1_to_dl2")

src/osa/utils/utils.py

@@ -8,13 +8,16 @@
 from datetime import datetime, timedelta
 from pathlib import Path
 from socket import gethostname
+from gammapy.data import observatory_locations
+from astropy import units as u
 
 import osa.paths
 from osa.configs import options
 from osa.configs.config import cfg
 from osa.utils.iofile import write_to_file
 from osa.utils.logging import myLogger
 
+
 __all__ = [
     "get_lstchain_version",
     "date_to_dir",
@@ -285,3 +288,21 @@ def wait_for_daytime(start=8, end=18):
     while time.localtime().tm_hour <= start or time.localtime().tm_hour >= end:
         log.info("Waiting for sunrise to not interfere with the data-taking. Sleeping.")
         time.sleep(3600)
+
+
+def culmination_angle(dec: int) -> float:  
+    """
+    Calculate culmination angle for a given declination.
+
+    Parameters
+    ----------
+    dec: int
+        declination in degrees
+
+    Returns
+    -------
+    Culmination angle in degrees
+    """
+    location = observatory_locations["cta_north"]
+    Lat = location.lat  # latitude of the LST1 site    
+    return abs(Lat - dec*u.deg).value