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

environment.yml

@@ -21,10 +21,10 @@ dependencies:
   - click
   - pymongo
   - scikit-learn=1.2
-  - ctapipe_io_lst=0.24
+  - ctapipe_io_lst >=0.25.1,<0.26a0
   - ctaplot~=0.6.4
   - pyirf~=0.10
-  - lstchain>=0.10.13
+  - lstchain>=0.10.15
   - tenacity
   - pip:
       - lstcam_calib

pyproject.toml

@@ -33,7 +33,7 @@ classifiers = [
 requires-python = ">=3.10"
 dependencies = [
     "astropy~=5.0",
-    "lstchain>=0.10.13",
+    "lstchain>=0.10.15",
     "matplotlib",
     "numpy",
     "pandas",
@@ -76,6 +76,7 @@ reprocess_longterm = "osa.scripts.reprocess_longterm:main"
 gain_selection = "osa.scripts.gain_selection:main"
 update_source_catalog = "osa.scripts.update_source_catalog:main"
 gainsel_webmaker = "osa.scripts.gainsel_webmaker:main"
+sequencer_catB_tailcuts = "osa.scripts.sequencer_catB_tailcuts:main"
 
 [tool.setuptools.packages.find]
 where = ["src"]

src/osa/configs/options.py

@@ -15,6 +15,8 @@
 warning = None
 nocheck = None
 no_dl2 = None
+no_dl1ab = None
+no_gainsel = None
 prod_id = None
 dl1_prod_id = None
 dl2_prod_id = None

src/osa/configs/sequencer.cfg

@@ -15,13 +15,16 @@ RUN_SUMMARY_DIR: %(MONITORING)s/RunSummary
 RUN_CATALOG: %(MONITORING)s/RunCatalog
 PEDESTAL_FINDER_DIR: %(BASE)s/auxiliary/PedestalFinder
 ANALYSIS_DIR: %(BASE)s/running_analysis
-CALIB_BASE_DIR: %(MONITORING)s/PixelCalibration/Cat-A
-CALIB_DIR: %(CALIB_BASE_DIR)s/calibration
-PEDESTAL_DIR: %(CALIB_BASE_DIR)s/drs4_baseline
+CALIB_BASE_DIR: %(MONITORING)s/PixelCalibration
+CAT_A_CALIB_BASE: %(CALIB_BASE_DIR)s/Cat-A
+CAT_A_CALIB_DIR: %(CAT_A_CALIB_BASE)s/calibration
+CAT_A_PEDESTAL_DIR: %(CAT_A_CALIB_BASE)s/drs4_baseline
+CAT_B_CALIB_BASE: %(CALIB_BASE_DIR)s/Cat-B
 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
@@ -45,22 +48,26 @@ DL2_PROD_ID: model2
 #charge_calibration: onsite_create_calibration_file
 drs4_baseline: lstcam_calib_onsite_create_drs4_pedestal_file
 charge_calibration: lstcam_calib_onsite_create_calibration_file
+catB_calibration: onsite_create_cat_B_calibration_file
 use_ff_heuristic_id: False
 
 # Data processing steps in datasequence script
 r0_to_dl1: lstchain_data_r0_to_dl1
 dl1ab: lstchain_dl1ab
 check_dl1: lstchain_check_dl1
 dl1_to_dl2: lstchain_dl1_to_dl2
+tailcuts_finder: lstchain_find_tailcuts
 
 # To be set by the user
 dl1a_config: /software/lstchain/data/lstchain_standard_config.json
 store_image_dl1ab: True
 merge_dl1_datacheck: True
+apply_catB_calibration: True
+apply_standard_dl1b_config: False
 use_ff_heuristic_gain_selection: False
 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: 20240918_v0.10.12
 dl3_config: /software/lstchain/data/dl3_std_config.json
 max_tries: 3
 

src/osa/conftest.py

@@ -574,10 +574,78 @@ 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"
     file.touch()
     return file
+
+
+@pytest.fixture(scope="session")
+def catB_closed_file(running_analysis_dir):
+
+    catB_closed_file = running_analysis_dir / "catB_00003.closed"
+    catB_closed_file.touch()
+    return catB_closed_file
+
+
+@pytest.fixture(scope="session")
+def catB_calib_base_dir(monitoring_dir):
+
+    catB_calib_base_dir = monitoring_dir / "PixelCalibration" / "Cat-B"
+    catB_calib_base_dir.mkdir(parents=True, exist_ok=True)
+    return catB_calib_base_dir
+
+
+@pytest.fixture(scope="session")
+def catB_calibration_file(catB_calib_dir):
+
+    catB_calib_dir = catB_calib_base_dir / "calibration" / nightdir / prod_id
+    catB_calib_file = catB_calib_dir / "cat_B_calibration_filters_52.Run00003.h5"
+    catB_calib_file.touch()
+    return catB_calib_file
+
+
+@pytest.fixture(scope="session")
+def dl1b_config_file(running_analysis_dir):
+    config_information = dedent(
+        """\
+            {
+            "tailcuts_clean_with_pedestal_threshold": {
+                "picture_thresh": 8,
+                "boundary_thresh": 4,
+                "sigma": 2.5,
+                "keep_isolated_pixels": false,
+                "min_number_picture_neighbors": 2,
+                "use_only_main_island": false,
+                "delta_time": 2
+            },
+            "dynamic_cleaning": {
+                "apply": true,
+                "threshold": 267,
+                "fraction_cleaning_intensity": 0.03
+            }
+        }"""
+    )
+    config_file = running_analysis_dir / "dl1ab_Run00003.json"
+    config_file.touch()
+    config_file.write_text(config_information)
+    return config_file
+

src/osa/job.py

@@ -20,10 +20,18 @@
     get_drive_file,
     get_summary_file,
     get_pedestal_ids_file,
+    get_dl1_prod_id_and_config,
+    get_dl2_nsb_prod_id,
 )
 from osa.utils.iofile import write_to_file
 from osa.utils.logging import myLogger
-from osa.utils.utils import date_to_dir, time_to_seconds, stringify, date_to_iso
+from osa.utils.utils import (
+    date_to_dir,
+    time_to_seconds,
+    stringify,
+    date_to_iso,
+    get_RF_model,
+)
 
 log = myLogger(logging.getLogger(__name__))
 
@@ -93,6 +101,14 @@ def are_all_jobs_correctly_finished(sequence_list):
     analysis_directory = Path(options.directory)
     for sequence in sequence_list:
         history_files_list = analysis_directory.rglob(f"*{sequence.seq}*.history")
+
+        if not options.test:
+            try:
+                next(history_files_list)
+            except StopIteration:
+                log.info("no history files found")
+                flag = False
+
         for history_file in history_files_list:
             # TODO: s.history should be SubRunObj attribute not RunObj
             # s.history only working for CALIBRATION sequence (run-wise), since it is
@@ -109,6 +125,12 @@ def are_all_jobs_correctly_finished(sequence_list):
                     f"finished up to DL1ab, but --no-dl2 option selected"
                 )
                 continue
+            if out == 3 and options.no_dl1ab:
+                log.debug(
+                    f"Job {sequence.seq} ({sequence.type}) correctly "
+                    f"finished up to DL1A, but --no-dl1ab option selected"
+                )
+                continue
 
             log.warning(
                 f"Job {sequence.seq} (run {sequence.run}) not correctly finished [level {out}]"
@@ -215,6 +237,8 @@ def historylevel(history_file: Path, data_type: str):
                 # Data sequence
                 elif program == cfg.get("lstchain", "r0_to_dl1"):
                     level = 3 if exit_status == 0 else 4
+                #elif program == cfg.get("lstchain", "catB_calibration"):
+                #    level = 3 if exit_status == 0 else 4
                 elif program == cfg.get("lstchain", "dl1ab"):
                     if (exit_status == 0) and (prod_id == options.dl1_prod_id):
                         log.debug(f"DL1ab prod ID: {options.dl1_prod_id} already produced")
@@ -425,6 +449,8 @@ def data_sequence_job_template(sequence):
         commandargs.extend(("--config", f"{Path(options.configfile).resolve()}"))
     if sequence.type == "DATA" and options.no_dl2:
         commandargs.append("--no-dl2")
+    if sequence.type == "DATA" and options.no_dl1ab:
+        commandargs.append("--no-dl1ab")
 
     commandargs.extend(
         (
@@ -439,6 +465,21 @@ def data_sequence_job_template(sequence):
         )
     )
 
+    if not options.no_dl1ab:
+        dl1_prod_id, dl1b_config = get_dl1_prod_id_and_config(sequence.run)
+        sequence.dl1_prod_id = dl1_prod_id
+        sequence.dl1b_config = dl1b_config
+
+        commandargs.append(f"--dl1b-config={sequence.dl1b_config}")
+        commandargs.append(f"--dl1-prod-id={sequence.dl1_prod_id}")
+
+    if not options.no_dl2 and not options.no_dl1ab:
+        sequence.rf_model = get_RF_model(sequence.run)
+        nsb_prod_id = get_dl2_nsb_prod_id(sequence.rf_model)
+        sequence.dl2_prod_id = f"{sequence.dl1_prod_id}/{nsb_prod_id}"
+        commandargs.append(f"--rf-model-path={sequence.rf_model}")
+        commandargs.append(f"--dl2-prod-id={sequence.dl2_prod_id}")
+
     content = job_header + "\n" + PYTHON_IMPORTS
 
     if not options.test:

src/osa/nightsummary/extract.py

@@ -23,7 +23,7 @@
 from osa.job import sequence_filenames
 from osa.nightsummary import database
 from osa.nightsummary.nightsummary import run_summary_table
-from osa.paths import sequence_calibration_files, get_run_date
+from osa.paths import sequence_calibration_files, get_run_date, get_dl1_prod_id_and_config, get_dl2_prod_id
 from osa.utils.logging import myLogger
 from osa.utils.utils import date_to_iso, date_to_dir
 
@@ -263,6 +263,14 @@ def extract_sequences(date: datetime, run_obj_list: List[RunObj]) -> List[Sequen
                 f"Data sequence {sequence.seq} from run {run.run} whose parent is "
                 f"{sequence.parent} (DRS4 {required_drs4_run} & Ped-Cal {required_pedcal_run})"
             )
+            if not options.no_dl1ab and sequence.type=="DATA":
+                dl1_prod_id, dl1b_config = get_dl1_prod_id_and_config(sequence.run)
+                sequence.dl1_prod_id = dl1_prod_id
+                sequence.dl1b_config = dl1b_config
+
+            if not options.no_dl2 and not options.no_dl1ab and sequence.type=="DATA":
+                sequence.dl2_prod_id = get_dl2_prod_id(sequence.run)
+
             sequence_list.append(sequence)
 
     # Add the calibration file names