Look for the closest declination line RF model to produce the DL2 files

src/osa/configs/sequencer.cfg

@@ -22,6 +22,7 @@ DL1_DIR: %(BASE)s/DL1
 DL1AB_DIR: %(BASE)s/DL1
 DL2_DIR: %(BASE)s/DL2
 DL3_DIR: %(BASE)s/DL3
+RF_MODELS: %(BASE)s/models/AllSky
 OSA_DIR: %(BASE)s/OSA
 CLOSER_DIR: %(OSA_DIR)s/Closer
 HIGH_LEVEL_DIR: %(OSA_DIR)s/HighLevel
@@ -56,7 +57,7 @@ store_image_dl1ab: True
 merge_dl1_datacheck: True
 dl1b_config: /software/lstchain/data/lstchain_standard_config.json
 dl2_config: /software/lstchain/data/lstchain_standard_config.json
-rf_models: /data/models/prod5/zenith_20deg/20201023_v0.6.3
+mc_prod: 20240131_allsky_v0.10.5_all_dec_base
 dl3_config: /software/lstchain/data/dl3_std_config.json
 max_tries: 3
 

src/osa/conftest.py

@@ -574,8 +574,23 @@ def database(osa_dir):
 
 
 @pytest.fixture(scope="session")
-def gain_selection_flag_file(osa_dir):
+def rf_models_allsky_basedir(base_test_dir):
+    directory = base_test_dir / "models/AllSky"
+    directory.mkdir(parents=True, exist_ok=True)
+    return directory
+
 
+@pytest.fixture(scope="session")
+def rf_model_path(rf_models_allsky_basedir):
+    mc_prod = "20240131_allsky_v0.10.5_all_dec_base"
+    declination_str = "dec_2276"
+    rf_model_path = rf_models_allsky_basedir / mc_prod / declination_str
+    rf_model_path.mkdir(parents=True, exist_ok=True)
+    return rf_model_path
+
+
+@pytest.fixture(scope="session")
+def gain_selection_flag_file(osa_dir):
     GainSel_dir = osa_dir / "GainSel" / "20200117"
     GainSel_dir.mkdir(parents=True, exist_ok=True)
     file = GainSel_dir / "GainSelFinished.txt"

src/osa/paths.py

@@ -7,9 +7,13 @@
 from typing import List
 import subprocess
 import time
+import os
 
 import lstchain
 from astropy.table import Table
+from astropy import units as u
+from astropy.coordinates import SkyCoord
+from gammapy.data import observatory_locations
 from lstchain.onsite import (find_systematics_correction_file,
                              find_time_calibration_file,
                              find_filter_wheels)
@@ -397,10 +401,98 @@ 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) -> u.Quantity:
+    """Return the declination angle in degrees corresponding to a 
+    given string of the form "dec_XXXX" or "dec_min_XXXX"."""
+
+    # Check if dec_str has a valid format
+    pattern = r'^dec_(\d{3,4})$|^dec_min_(\d{3,4})$'
+    if re.match(pattern, dec_str):
+
+        # 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*u.deg
+
+
+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_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(cfg.get(options.tel_id, "RUN_CATALOG"))
+    run_catalog_file = run_catalog_dir / f"RunCatalog_{utils.date_to_dir(options.date)}.ecsv"
+    run_catalog = Table.read(run_catalog_file)
+    run = run_catalog[run_catalog["run_id"]==int(run_str)]
+    target_name = run["source_name"]
+
+    if options.test or options.simulate:
+        source_coordinates = SkyCoord.from_name(target_name[0])
+        source_dec = source_coordinates.dec
+    else:
+        tcu_server = cfg.get("database", "tcu_db")
+        source_dec = utils.get_source_dec_from_TCU(target_name[0], tcu_server)
+
+    source_culmination = utils.culmination_angle(source_dec)
+
+    rf_models_dir = Path(cfg.get("LST1", "RF_MODELS"))
+    mc_prod = cfg.get("lstchain", "mc_prod")
+    dec_list = os.listdir(rf_models_dir / mc_prod)
+
+    # Convert each string in the list to numerical values
+    dec_values = [convert_dec_string(dec) for dec in dec_list]
+    dec_values = [dec for dec in dec_values if dec is not None]
+
+    closest_declination = min(dec_values, key=lambda x: abs(x - source_dec))
+    closest_dec_culmination = utils.culmination_angle(closest_declination)
+
+    location = observatory_locations["cta_north"]
+    Lat = location.lat  # latitude of the LST1 site    
+    closest_lines = sorted(sorted(dec_values, key=lambda x: abs(x - Lat))[:2])
+
+    if source_dec < closest_lines[0] or source_dec > closest_lines[1]:
+        # If the source declination is between the two MC lines closest to the latitude of 
+        # the LST1 site, this check is not necessary.
+        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."
+        )
+        while 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)
+            declination_str = corresponding_dict[closest_declination]
+            dec_values.remove(closest_declination)
+            dec_list.remove(declination_str)
+            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 line to use for the DL2 production is: {closest_declination}")
+
+    corresponding_dict = get_corresponding_string(dec_list, dec_values)
+    declination_str = corresponding_dict[closest_declination]
+
+    rf_model_path = rf_models_dir / mc_prod / declination_str
+
+    return rf_model_path

src/osa/provenance/utils.py

@@ -40,7 +40,7 @@ def parse_variables(class_instance):
     configfile_dl1b = cfg.get("lstchain", "dl1b_config")
     configfile_dl2 = cfg.get("lstchain", "dl2_config")
     raw_dir = Path(cfg.get("LST1", "R0_DIR"))
-    rf_models_directory = Path(cfg.get("lstchain", "RF_MODELS"))
+    rf_models_directory = Path(cfg.get("LST1", "RF_MODELS"))
     dl1_dir = Path(cfg.get("LST1", "DL1_DIR"))
     dl2_dir = Path(cfg.get("LST1", "DL2_DIR"))
     calib_dir = Path(cfg.get("LST1", "CALIB_DIR"))

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[:5])
     dl1_file = dl1ab_subdirectory / f"dl1_LST-1.Run{run_str}.h5"
 
     command = cfg.get("lstchain", "dl1_to_dl2")

src/osa/scripts/tests/test_osa_scripts.py

@@ -106,10 +106,10 @@ def test_simulate_processing(
 
     with open(json_file_dl2) as file:
         dl2 = yaml.safe_load(file)
-    assert len(dl2["entity"]) == 25
-    assert len(dl2["activity"]) == 6
-    assert len(dl2["used"]) == 21
-    assert len(dl2["wasGeneratedBy"]) == 12
+    assert len(dl2["entity"]) == 19
+    assert len(dl2["activity"]) == 5
+    assert len(dl2["used"]) == 15
+    assert len(dl2["wasGeneratedBy"]) == 10
 
     rc = run_program("simulate_processing", "-p")
     assert rc.returncode == 0
@@ -240,17 +240,28 @@ def test_closer(
     assert closed_seq_file.exists()
 
 
-def test_datasequence(running_analysis_dir):
+def test_datasequence(
+    running_analysis_dir,
+    run_catalog,
+    run_catalog_dir,
+    rf_models_allsky_basedir,
+    rf_model_path
+):
     drs4_file = "drs4_pedestal.Run00001.0000.fits"
     calib_file = "calibration.Run00002.0000.hdf5"
     timecalib_file = "time_calibration.Run00002.0000.hdf5"
     systematic_correction_file = "no_sys_corrected_calibration_scan_fit_20210514.0000.h5"
     drive_file = "DrivePosition_20200117.txt"
     runsummary_file = "RunSummary_20200117.ecsv"
     prod_id = "v0.1.0"
-    run_number = "00003.0000"
+    run_number = "01807.0000"
     options.directory = running_analysis_dir
 
+    assert run_catalog_dir.exists()
+    assert run_catalog.exists()
+    assert rf_models_allsky_basedir.exists()
+    assert rf_model_path.exists()
+
     output = run_program(
         "datasequence",
         "--date=2020-01-17",

src/osa/utils/utils.py

@@ -8,13 +8,17 @@
 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
+from pymongo import MongoClient
 
 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 +289,29 @@ 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: u.Quantity) -> u.Quantity:  
+    """
+    Calculate culmination angle for a given declination.
+
+    Parameters
+    ----------
+    dec: Quantity
+        declination coordinate in degrees
+
+    Returns
+    -------
+    Culmination angle in degrees
+    """
+    location = observatory_locations["cta_north"]
+    Lat = location.lat  # latitude of the LST1 site    
+    return abs(Lat - dec)
+
+
+def get_source_dec_from_TCU(source_name: str, tcu_server: str) -> u.Quantity:
+    """Get the declination of a given source from the TCU database."""
+    client = MongoClient(tcu_server)
+    collection = client["lst1_config"]["structure_configurations"]
+    source_dec = collection.find_one({"name": source_name})["target"]["source_dec"]
+    return source_dec*u.deg