diff --git a/src/nectarchain/acceptance_verification_package/deadtime.py b/src/nectarchain/acceptance_verification_package/deadtime.py new file mode 100644 index 00000000..565fe04c --- /dev/null +++ b/src/nectarchain/acceptance_verification_package/deadtime.py @@ -0,0 +1,308 @@ +import argparse +import copy +import json +import logging +import os +import sys +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +from astropy import units as u +from ctapipe.containers import EventType + +from nectarchain.trr_test_suite.deadtime import ( + fit_rate_per_run, + plot_deadtime_vs_collected_trigger_rate, + plot_fitted_rate_vs_collected_trigger_rate, + run_deadtime_test_tool_process, +) +from nectarchain.trr_test_suite.utils import plot_deadtime_and_expo_fit +from nectarchain.utils.constants import ALLOWED_CAMERAS + +logging.basicConfig( + format="%(asctime)s %(name)s %(levelname)s %(message)s", + filename=f"{os.environ.get('NECTARCHAIN_LOG', '/tmp')}/{os.getpid()}/" + f"{Path(__file__).stem}_{os.getpid()}.log", + handlers=[logging.getLogger("__main__").handlers], +) +log = logging.getLogger(__name__) + +plt.style.use( + os.path.join( + os.path.abspath(os.path.dirname(__file__)), "../utils/plot_style.mpltstyle" + ) +) + + +def get_labels(): + with open("resources/source_type_labels.json", "r") as f: + source_labels = json.load(f) + + return source_labels + + +def get_args(): + """Parses command-line arguments for the deadtime test script. + + Returns + ------- + parser : argparse.ArgumentParser + The parsed command-line arguments. + """ + parser = argparse.ArgumentParser( + description="Deadtime tests B-TEL-1260 & B-TEL-1270. \n" + + "According to the nectarchain component interface, you have to set a \ + NECTARCAMDATA environment variable in the folder where you have the data \ + from your runs or where you want them to be downloaded.\n" + + "You have to provide a run number (or list of numbers), the event source, \ + a corresponding camera tag and, optionally, the number of events \ + to consider for the test and an output directory \ + to save the final plots.\n" + + "If the data is not in NECTARCAMDATA, the files will be downloaded through \ + DIRAC.\n" + + "You can optionally specify the number of events to be processed \ + (default 8000).\n", + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) + parser.add_argument( + "-r", + "--run_numbers", + type=int, + nargs="+", + help="Run number (or list of numbers)", + required=True, + default=[i for i in range(3332, 3351)] + [i for i in range(3552, 3563)], + ) + parser.add_argument( + "-s", + "--source", + type=int, + nargs="+", + choices=[member.value for member in EventType], + # [0, 1, 2, 3, 4, 15, 16, 17, 24, 32, 255] + help="Source number (or list of numbers)", + required=True, + default=[ + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 32, + 32, + 32, + 32, + 32, + 32, + 32, + 32, + 32, + 32, + 32, + 32, + 32, + 32, + 32, + 32, + 32, + 32, + 32, + 32, + ], + ) + parser.add_argument( + "-e", + "--evts", + type=int, + help="Number of events to process from each run.", + required=False, + default=8000, + ) + parser.add_argument( + "-c", + "--camera", + choices=ALLOWED_CAMERAS, + default=[camera for camera in ALLOWED_CAMERAS if "QM" in camera][0], + help="Process data for a specific NectarCAM camera.", + type=str, + ) + parser.add_argument( + "-o", + "--output", + type=str, + help="Output directory", + default=f"{os.environ.get('NECTARCHAIN_FIGURES', f'/tmp/{os.getpid()}')}", + ) + parser.add_argument( + "--temp_output", help="Temporary output directory for GUI", default=None + ) + parser.add_argument( + "-l", + "--log", + help="log level", + default="info", + type=str, + ) + + return parser + + +def main(): + """Runs the deadtime test script, which performs deadtime tests B-TEL-1260 and + B-TEL-1270, and event rate test B-MST-1280. + + The script takes command-line arguments to specify a run number, \ + corresponding event source, the camera tag, and, optionally, \ + the number of events to consider for the test and an output directory. \ + It then processes the data for each run, performs an exponential \ + fit to the deadtime distribution, and generates three plots: + + 1. A plot of the exponential function fit on the deadtime\ + distribution for each run. + 2. A plot of deadtime percentage vs. collected trigger rate, with the CTAO\ + requirement indicated. + 3. A plot of the rate from the fit vs. the collected trigger rate, with the\ + relative difference shown in the bottom panel. + + The script also saves the generated plots to the specified output directory, and\ + optionally saves the last two to a temporary output directory for use\ + in a GUI. + """ + + parser = get_args() + args = parser.parse_args() + log.setLevel(args.log.upper()) + + runlist = args.run_numbers + ids = args.source + + assert len(runlist) == len(ids), "'runlist' and 'ids' must have the same length" + + source_labels = get_labels() + + nevents = args.evts + + kwargs = copy.deepcopy(vars(args)) + kwargs.pop("camera") + camera = args.camera + + output_dir = os.path.join( + os.path.abspath(args.output), + f"av_camera_{camera}/{Path(__file__).stem}", + ) + os.makedirs(output_dir, exist_ok=True) + temp_output = os.path.abspath(args.temp_output) if args.temp_output else None + + log.info(f"Running the script with arguments: {kwargs}") + + # Drop arguments from the script after they are parsed, for the GUI to work properly + sys.argv = sys.argv[:1] + + ( + _, + _, + event_counter, + busy_counter, + collected_trigger_rates, + time_tot, + deadtime_us, + deadtime_pc, + ) = run_deadtime_test_tool_process( + runlist=runlist, nevents=nevents, ids=ids, test_type="av" + ) + + results = fit_rate_per_run(runlist=runlist, deadtime_us=deadtime_us)[-1] + + log.info(f"Output directory: {output_dir}") + log.info(f"Temporary output file: {temp_output}") + log.info(f"N max events to be considered: {nevents}") + log.info("-" * 40) + for ii, (key, values) in enumerate(results.items()): + log.info(f"For run {key}, source: {ids[ii]},") + log.info( + "Dead-Time extracted from the tool process: " + f"{np.min(deadtime_us[ii]):.3f}" + ) + log.info(f"Dead-Time from the fit: {values[0]:.3f} +- " f"{values[1]:.3f} µs") + log.info(f"Rate from the fit: {values[2]:.2f} +- " f"{values[3]:.2f} Hz") + log.info("Expected run duration from the fit: " f"{values[4]:.2f} s") + log.info("-" * 40) + + ids = np.array(ids) + runlist = np.array(runlist) + + error_deadtime_pc = [] + for run_id in range(np.array(busy_counter).shape[0]): + error_deadtime_pc.append( + np.sqrt( + (busy_counter[run_id][-1] * event_counter[run_id][-1]) + / ((busy_counter[run_id][-1] + event_counter[run_id][-1]) ** 3) + ) + ) + error_deadtime_pc = np.array(error_deadtime_pc) + + deadtime, fitted_trigger_rates, fitted_trigger_rates_err = [], [], [] + + for ii, run_num in enumerate(runlist): + results = plot_deadtime_and_expo_fit( + total_delta_t_for_busy_time=time_tot[ii], + deadtime_us=np.array(deadtime_us[ii].value), + run=run_num, + output_plot=output_dir, + run_type=EventType(ids[ii]).name, + temp_output=temp_output, + ) + deadtime.append(results[0]) + fitted_trigger_rates.append(((-1 * results[6]) * (1 / u.us)).to(u.kHz).value) + fitted_trigger_rates_err.append(((results[8]) * (1 / u.us)).to(u.kHz).value) + plt.close() + + deadtime = np.array(deadtime) + fitted_trigger_rates = np.array(fitted_trigger_rates) + fitted_trigger_rates_err = np.array(fitted_trigger_rates_err) + + deadtime_pc_fit = np.array( + [ + # the parameter_lambda is a rate value in kHz, + # so one needs to compare the deadtime in mus with the rate in kHz + # and finally make it a percentage value + deadtime[ii] * rate * 1e2 * 1e-3 + for ii, rate in enumerate(fitted_trigger_rates) + ] + ) + + if len(runlist) > 1: + plot_deadtime_vs_collected_trigger_rate( + np.array(ids), + np.array(collected_trigger_rates), + np.array(deadtime_pc), + np.array(error_deadtime_pc), + np.array(deadtime_pc_fit), + False, + source_labels, + output_dir, + temp_output, + test_type="av", + ) + plot_fitted_rate_vs_collected_trigger_rate( + np.array(ids), + np.array(collected_trigger_rates), + np.array(fitted_trigger_rates), + np.array(fitted_trigger_rates_err), + source_labels, + output_dir, + temp_output, + test_type="av", + ) + + plt.close("all") + + +if __name__ == "__main__": + main() diff --git a/src/nectarchain/acceptance_verification_package/gui.py b/src/nectarchain/acceptance_verification_package/gui.py index 2378997f..696248e0 100644 --- a/src/nectarchain/acceptance_verification_package/gui.py +++ b/src/nectarchain/acceptance_verification_package/gui.py @@ -2,7 +2,7 @@ from PyQt5.QtWidgets import QApplication -from nectarchain.acceptance_verification_package import hillas_validation +from nectarchain.acceptance_verification_package import deadtime, hillas_validation from nectarchain.trr_test_suite.gui import TestRunner @@ -10,6 +10,7 @@ class AcceptanceTestRunner(TestRunner): # redefine list of test modules test_modules = { "Hillas parameter validation": hillas_validation, + "Deadtime Test": deadtime, } diff --git a/src/nectarchain/acceptance_verification_package/resources/source_type_labels.json b/src/nectarchain/acceptance_verification_package/resources/source_type_labels.json new file mode 100644 index 00000000..ea039a61 --- /dev/null +++ b/src/nectarchain/acceptance_verification_package/resources/source_type_labels.json @@ -0,0 +1,13 @@ +{ + "0": {"source": "Flatfield", "color": "red"}, + "1": {"source": "Single Photo Electron", "color": "blue"}, + "2": {"source": "Sky Pedestal", "color": "purple"}, + "3": {"source": "Dark Pedestal", "color": "green"}, + "4": {"source": "Electronic Pedestal", "color": "orange"}, + "15": {"source": "Other Calibration", "color": "cyan"}, + "16": {"source": "Muon", "color": "magenta"}, + "17": {"source": "Hardware Stereo", "color": "yellow"}, + "24": {"source": "DAQ software trigger", "color": "brown"}, + "32": {"source": "Standard Physics stereo trigger", "color": "pink"}, + "255": {"source": "Unknown", "color": "black"} +} \ No newline at end of file diff --git a/src/nectarchain/trr_test_suite/deadtime.py b/src/nectarchain/trr_test_suite/deadtime.py index 928837f9..941993ad 100644 --- a/src/nectarchain/trr_test_suite/deadtime.py +++ b/src/nectarchain/trr_test_suite/deadtime.py @@ -43,6 +43,7 @@ def plot_deadtime_vs_collected_trigger_rate( labels, output_dir, temp_output, + test_type, ): """Plot the deadtime percentage vs the trigger rates @@ -62,19 +63,31 @@ def plot_deadtime_vs_collected_trigger_rate( show_camera_client : bool Whether to show the shadead areas for the percentage values extracted from the camera client counters - labels : list + labels : dict Labels with the source names for the plot output_dir : str Path to the output directory to save the plot temp_output : str Path to temporary output directory for the GUI + test_type : str + Test type to specify the source ids. + Accepted options are 'trr' and 'av', + for 'Test-Readiness Review' and 'Acceptance Verification'. """ + if test_type not in ["trr", "av"]: + log.warning("Invalid chosen 'test_type', falling back to 'trr'.") + test_type = "trr" + available_ids = np.unique(ids).tolist() + fig, ax = plt.subplots() - for source in range(0, 3): + for source in available_ids: mask_source = np.where(ids == source)[0] + if test_type == "av": + source = str(source) + if show_camera_client: deadtime_pc_source = deadtime_pc[mask_source] error_deadtime_pc_source = error_deadtime_pc[mask_source] @@ -121,7 +134,7 @@ def plot_deadtime_vs_collected_trigger_rate( ax.text( 3.5, 6.75, - "CTA requirement", + "CTAO requirement", color="gray", fontsize=20, alpha=0.7, @@ -129,7 +142,7 @@ def plot_deadtime_vs_collected_trigger_rate( verticalalignment="center", ) - plt.legend() + plt.legend(fontsize=12) plt.xlim(0.0, 15) plt.yscale("log") @@ -152,6 +165,7 @@ def plot_fitted_rate_vs_collected_trigger_rate( labels, output_dir, temp_output, + test_type, ): """Plot the fitted vs collected trigger rates and the relative difference @@ -166,14 +180,23 @@ def plot_fitted_rate_vs_collected_trigger_rate( lambda_from_fit_err : np.ndarray Error on lambda parameter values from the exponential fit with `fit_rate_per_run` - labels : list + labels : dict Labels with the source names for the plot output_dir : str Path to the output directory to save the plot temp_output : str Path to temporary output directory for the GUI + test_type : str + Test type to specify the source ids. + Accepted options are 'trr' and 'av', + for 'Test-Readiness Review' and 'Acceptance Verification'. """ + if test_type not in ["trr", "av"]: + log.warning("Invalid chosen 'test_type', falling back to 'trr'.") + test_type = "trr" + available_ids = np.unique(ids).tolist() + fig, ((ax1, ax2)) = plt.subplots( 2, 1, @@ -190,8 +213,12 @@ def plot_fitted_rate_vs_collected_trigger_rate( x_err = 0 err_ratio = relative * (((rate_err + x_err) / (yy - xx)) + x_err / xx) - for source in range(0, 3): + for source in available_ids: runl = np.where(ids == source)[0] + + if test_type == "av": + source = str(source) + ax2.errorbar( xx[runl], relative[runl], @@ -203,7 +230,7 @@ def plot_fitted_rate_vs_collected_trigger_rate( color=labels[source]["color"], ) - ax2.set_ylim(-25, 25) + ax2.set_ylim(-15, 15) xx = range(0, 60) @@ -224,9 +251,12 @@ def plot_fitted_rate_vs_collected_trigger_rate( ax1.set_ylim(1e0, 60) ax2.set_xlim(1e0, 60) - for source in range(0, 3): + for source in available_ids: runl = np.where(ids == source)[0] + if test_type == "av": + source = str(source) + ax1.errorbar( collected_trigger_rates[runl] / 1000, rate[runl], @@ -238,7 +268,7 @@ def plot_fitted_rate_vs_collected_trigger_rate( label=labels[source]["source"], ) - ax1.legend() + ax1.legend(fontsize=12) plot_path = os.path.join(output_dir, "fitted_vs_collected_trigger_rates.png") plt.savefig(plot_path) @@ -389,8 +419,10 @@ def fit_rate_per_run(runlist: list, deadtime_us: np.ndarray): ) -def run_deadtime_test_tool_process(runlist: list, nevents: int, ids: np.ndarray): - """Run `DeadtimeTestTool` from `utils.py` over the provided run list +def run_deadtime_test_tool_process( + runlist: list, nevents: int, ids: np.ndarray, test_type: str = "trr" +): + """Run `DeadtimeTestTool` from `tools_components.py` over the provided run list Parameters ---------- @@ -399,7 +431,11 @@ def run_deadtime_test_tool_process(runlist: list, nevents: int, ids: np.ndarray) nevents : int max number of events ids : np.ndarray - Source ids for all the runs. + Source ids for all the runs + test_type : str + Test type to specify the source ids. + Accepted options are 'trr' and 'av', + for 'Test-Readiness Review' and 'Acceptance Verification'. Returns ------- @@ -421,12 +457,18 @@ def run_deadtime_test_tool_process(runlist: list, nevents: int, ids: np.ndarray) The deadtime percentage value computed with the counters for each run """ + if test_type not in ["trr", "av"]: + log.warning("Invalid chosen 'test_type', falling back to 'trr'.") + test_type = "trr" + ucts_timestamps, ucts_deltat = [], [] event_counter, busy_counter = [], [] collected_trigger_rates = [] time_tot = [] deadtime_us, deadtime_pc = [], [] + log.info(f"Starting `DeadtimeTestTool` for test {test_type}") + for run, id in zip(runlist, ids): log.info("Processing `DeadtimeTestTool` on run {}".format(run)) tool = DeadtimeTestTool( @@ -442,7 +484,7 @@ def run_deadtime_test_tool_process(runlist: list, nevents: int, ids: np.ndarray) tool.initialize() tool.setup() tool.start() - output = tool.finish(id=id) + output = tool.finish(id=id, test_type=test_type) ucts_timestamps.append(output[0]) ucts_deltat.append(output[1]) diff --git a/src/nectarchain/trr_test_suite/tools_components.py b/src/nectarchain/trr_test_suite/tools_components.py index cd35fd05..2f1b8739 100644 --- a/src/nectarchain/trr_test_suite/tools_components.py +++ b/src/nectarchain/trr_test_suite/tools_components.py @@ -466,16 +466,20 @@ class DeadtimeTestTool(EventsLoopNectarCAMCalibrationTool): def finish(self, *args, **kwargs): id = kwargs.pop("id") + test_type = kwargs.pop("test_type") output = super().finish(return_output_component=True, *args, **kwargs) # Specify event type - # NOTE: will probably need to be revisited - if id == 0: # FFCLS - event_type = EventType.FLATFIELD - elif id == 1: # NSB - event_type = EventType.SUBARRAY - elif id == 2: # Laser - event_type = EventType.SUBARRAY + # NOTE: may probably need to be revisited + if test_type == "trr": + if id == 0: # FFCLS + event_type = EventType.FLATFIELD + elif id == 1: # NSB + event_type = EventType.SUBARRAY + elif id == 2: # Laser + event_type = EventType.SUBARRAY + elif test_type == "av": + event_type = EventType(id) charge_container = output[0].containers[event_type] diff --git a/src/nectarchain/trr_test_suite/utils.py b/src/nectarchain/trr_test_suite/utils.py index 8ccef1be..6bdbe3e1 100644 --- a/src/nectarchain/trr_test_suite/utils.py +++ b/src/nectarchain/trr_test_suite/utils.py @@ -1,5 +1,6 @@ import logging import os +import pickle import matplotlib.pyplot as plt import numpy as np @@ -559,7 +560,13 @@ def pois(x, A, R): # function by Federica def plot_deadtime_and_expo_fit( - total_delta_t_for_busy_time, deadtime_us, run, verbose=False, output_plot=None + total_delta_t_for_busy_time, + deadtime_us, + run, + verbose=False, + output_plot=None, + run_type=None, + temp_output=None, ): """Compute the deadtime and exponential fit parameters for a given dataset. @@ -575,6 +582,11 @@ def plot_deadtime_and_expo_fit( Whether to print the fit results or not. output_plot : str, optional The path to save the output plot. + run_type : str, optional + The run type, as extracted from ctapipe.containters.EventType + temp_output : str, optional + The path for the temporary output directory + for the GUI for the deadtime tests Returns ------- @@ -656,7 +668,7 @@ def plot_deadtime_and_expo_fit( parameter_R_err_new = result.params["R"].stderr if output_plot: - _, ax = plt.subplots(1, 1, figsize=(10, 10 / 1.61), layout="constrained") + fig, ax = plt.subplots(1, 1, figsize=(10, 10 / 1.61), layout="constrained") # plot deltaT distribution plt.hist( deadtime_us, @@ -689,6 +701,7 @@ def plot_deadtime_and_expo_fit( max_x_value_for_plot - 50, np.max(entries) - 10, f"Run {run}" + + (f", {run_type}" if run_type else "") + "\n" + r"$f(\Delta t) = A \cdot$" + r"$\exp({-R \cdot (\delta t - \Delta_{\mathrm{min}})})$" @@ -726,6 +739,12 @@ def plot_deadtime_and_expo_fit( plt.savefig(plot_path) log.info(f"Plot saved at {plot_path}") + if temp_output: + with open( + os.path.join(temp_output, f"plot_deadtime_expo_fit_run{run}.pkl"), "wb" + ) as f: + pickle.dump(fig, f) + return ( deadtime, deadtime_bin,