From 95da131a1bbf627b5fd390e51b1d0cd02e8c6c2e Mon Sep 17 00:00:00 2001 From: MOUSADI Despoina Date: Wed, 9 Oct 2024 11:25:55 +0200 Subject: [PATCH] Converted test_scripts submodule to a regular directory (ignore previous 2 commits) --- .../user_scripts/dmousadi/test_scripts | 1 - .../dmousadi/test_scripts/README.md | 101 ++ .../user_scripts/dmousadi/test_scripts/gui.py | 443 ++++++ .../dmousadi/test_scripts/requirements.txt | 16 + .../dmousadi/test_scripts/tests/__init__.py | 0 .../test_scripts/tests/deadtime_test.py | 559 +++++++ .../test_scripts/tests/linearity_test.py | 292 ++++ .../test_scripts/tests/pedestal_test.py | 127 ++ .../tests/pix_couple_tim_uncertainty_test.py | 131 ++ .../tests/pix_tim_uncertainty_test.py | 305 ++++ .../tests/test_tools_components.py | 1291 +++++++++++++++++ .../dmousadi/test_scripts/tests/utils.py | 500 +++++++ 12 files changed, 3765 insertions(+), 1 deletion(-) delete mode 160000 src/nectarchain/user_scripts/dmousadi/test_scripts create mode 100644 src/nectarchain/user_scripts/dmousadi/test_scripts/README.md create mode 100644 src/nectarchain/user_scripts/dmousadi/test_scripts/gui.py create mode 100644 src/nectarchain/user_scripts/dmousadi/test_scripts/requirements.txt create mode 100644 src/nectarchain/user_scripts/dmousadi/test_scripts/tests/__init__.py create mode 100644 src/nectarchain/user_scripts/dmousadi/test_scripts/tests/deadtime_test.py create mode 100644 src/nectarchain/user_scripts/dmousadi/test_scripts/tests/linearity_test.py create mode 100644 src/nectarchain/user_scripts/dmousadi/test_scripts/tests/pedestal_test.py create mode 100644 src/nectarchain/user_scripts/dmousadi/test_scripts/tests/pix_couple_tim_uncertainty_test.py create mode 100644 src/nectarchain/user_scripts/dmousadi/test_scripts/tests/pix_tim_uncertainty_test.py create mode 100644 src/nectarchain/user_scripts/dmousadi/test_scripts/tests/test_tools_components.py create mode 100644 src/nectarchain/user_scripts/dmousadi/test_scripts/tests/utils.py diff --git a/src/nectarchain/user_scripts/dmousadi/test_scripts b/src/nectarchain/user_scripts/dmousadi/test_scripts deleted file mode 160000 index defff7f4..00000000 --- a/src/nectarchain/user_scripts/dmousadi/test_scripts +++ /dev/null @@ -1 +0,0 @@ -Subproject commit defff7f4713e6a986f0f2b82e1dcf656203f23eb diff --git a/src/nectarchain/user_scripts/dmousadi/test_scripts/README.md b/src/nectarchain/user_scripts/dmousadi/test_scripts/README.md new file mode 100644 index 00000000..beec9ddb --- /dev/null +++ b/src/nectarchain/user_scripts/dmousadi/test_scripts/README.md @@ -0,0 +1,101 @@ + +# NectarCAM Test Scripts + +Automated tests from the Test Readiness Review document for the CTA NectarCAM based on [nectarchain](https://github.com/cta-observatory/nectarchain) and [ctapipe](https://github.com/cta-observatory/ctapipe). This project can test for the following requirements: +- B-TEL-1010 Intensity resolution & B-TEL-1390 Linearity (script [linearity_test.py](tests/linearity_test.py)) +- B-TEL-1030 Time resolution (script [pix_couple_tim_uncertainty_test.py](tests/pix_couple_tim_uncertainty_test.py)) +- B-TEL-1260 Deadtime, B-TEL-1270 Deadtime Measurement & B-MST-1280 Event Rate (script [deadtime_test.py](tests/deadtime_test.py)) +- B-TEL-1370 Pedestal subtraction (script [pedestal_test.py](tests/pedestal_test.py)) +- B-TEL-1380 Systematic pixel timing uncertainty (script [pix_tim_uncertainty_test.py](tests/pix_tim_uncertainty_test.py)) + + +## Table of Contents + +- [Installation](#installation) +- [Usage](#usage) +- [Features](#features) +- [Contributing](#contributing) + + +## Installation + +Instructions on how to install and set up the project: +``` +git clone https://drf-gitlab.cea.fr/dmousadi/camera-test-scripts +cd camera-test-scripts +pip install -r requirements.txt +``` +If you want to automatically download your data, one of the requirements is also DIRAC, for which you need to have a grid certificate. It is not necessary for this repo, if you have your NectarCAM runs (fits files) locally stored. You can find more information about DIRAC [here](https://gitlab.cta-observatory.org/cta-computing/dpps/workload/CTADIRAC). If you are installing these packages for the first time and getting 'error building wheel', you might need to (re)install some of these: swig, ca-certificates, openssl, boost, protobuff, cmake. + +Once you have set up your environment, if you're not already a nectarchain user you need to set the NECTARCAMDATA environment variable to the directory where you have the NectarCAM runs: +``` +conda activate your_env_name +mkdir -p $CONDA_PREFIX/etc/conda/activate.d +mkdir -p $CONDA_PREFIX/etc/conda/deactivate.d +nano $CONDA_PREFIX/etc/conda/activate.d/env_vars.sh +nano $CONDA_PREFIX/etc/conda/deactivate.d/env_vars.sh +``` +In activate.d/env_vars.sh add: +``` +export NECTARCAMDATA=your_path_to_nectarcamdata +``` +And then in deactivate.d/env_vars.sh: +``` +unset NECTARCAMDATA +``` + + +According to nectarchain, the directory should have this structure: +``` +NECTARCAMDATA/runs/{fits_files} +``` +If the structure is different in your server (testbench server, looking at you) you could do a symbolic link of the form: +``` +ln -s /data/ZFITS/*/*/* nectarcamdatafolder/runs +``` + + +## Usage + +The test scripts can be run autonomously, or through a GUI. To run the gui you can just start it with python: +``` +python gui.py +``` +In the gui you can pick the test and set the parameters needed or keep the default ones. The most important parameter in every test is runlist, where you state which runs you want to process. If you have this runs in your NECTARCAMDATA, they will be found there and processed, otherwise if you don't have them but you're using DIRAC, they will be downloaded there. All tests have an output parameter, where you state where you want to save the plots produced by the tests. The evts parameter indicates the amount of events that are going to be processed (more events = more accurate results but more runtime). The linearity test also requires the corresponding transmission for each of the runs in the list, while in the timing uncertainty between couples of pixels you have to add the path to the csv file containing the PMT transit time values (calculated by Federica, might be stored in a database later). +To have a better look at the parameters, you can also run the tests on their own with the help option: +``` +python tests/linearity_test.py --help +``` +When everything is ready, you can press on run test. + + +## Features + +List of key features of the project: + +- GUI to run tests automatically +- If runs aren't already present, they will be downloaded with DIRAC +- Tests process runs and create a plot showing the results. Intermediate results for each separate run are saved in HDF5 tables in a folder in NECTARCAMDATA. The plot is saved in a specified folder, as well as in a temporary file which is read by the GUI +- Plot display with zoom and navigation +- Output logs and parameters displayed dynamically + +The test scripts essentially loop through the runs using the ctapipe Tool and Component classes, and then process the runs all together to calculate the results that go into the final plot. + +## Contributing + +You can contribute by adding more tests. The full list of tests is [here](https://docs.google.com/spreadsheets/d/1t5Z9ZHRESB6BYzmMbyFCDDALgTqvSb1KpJ8JMEiZgnI/edit?usp=sharing). It is planned to discuss about the implementation of more of the tests. + +Guidelines for contributing to the project: + +1. Fork the repository +2. Create a new branch (`git checkout -b feature-branch`) +3. Make your changes +4. Commit your changes (`git commit -am 'Add new feature'`) +5. Push to the branch (`git push origin feature-branch`) +6. Create a new Pull Request + +## Notes + +This project uses a lot of methods developped by Federica Bradascio, and has also taken inspiration from Guillaume Grolleron's tutorials in nectarchain about using the ctapipe Tool class. + +For the in-code documentation I have used Cody (Copilot's free-ish cousin). I have checked everything but I still could have missed some mistake in explaining the code's functions. diff --git a/src/nectarchain/user_scripts/dmousadi/test_scripts/gui.py b/src/nectarchain/user_scripts/dmousadi/test_scripts/gui.py new file mode 100644 index 00000000..1cba6059 --- /dev/null +++ b/src/nectarchain/user_scripts/dmousadi/test_scripts/gui.py @@ -0,0 +1,443 @@ +""" +The `TestRunner` class is a GUI application that allows the user to run various tests and display the results. + +The class provides the following functionality: +- Allows the user to select a test from a dropdown menu. +- Dynamically generates input fields based on the selected test. +- Runs the selected test and displays the output in a text box. +- Displays the test results in a plot canvas, with navigation buttons to switch between multiple plots. +- Provides a dark-themed UI with custom styling for various UI elements. + +The class uses the PyQt5 library for the GUI implementation and the Matplotlib library for plotting the test results. +""" +import sys +import os +import argparse +import tempfile +from PyQt5.QtWidgets import (QApplication, QWidget, QVBoxLayout, QLabel, + QComboBox, QLineEdit, QPushButton, QMessageBox, + QTextEdit, QHBoxLayout, QFileDialog, QSpacerItem, QSizePolicy, + QGroupBox) +from PyQt5.QtCore import QProcess, QTimer +from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas +from matplotlib.backends.backend_qt5 import NavigationToolbar2QT as NavigationToolbar +from matplotlib.figure import Figure +import matplotlib.pyplot as plt +import pickle + +# Ensure the tests directory is in sys.path +test_dir = os.path.abspath('tests') +if test_dir not in sys.path: + sys.path.append(test_dir) + +# Import test modules +import linearity_test +import deadtime_test +import pedestal_test +import pix_couple_tim_uncertainty_test +import pix_tim_uncertainty_test + +class TestRunner(QWidget): + def __init__(self): + super().__init__() + self.params = {} + self.process = None + self.plot_files = [] # Store the list of plot files + self.current_plot_index = 0 # Index to track which plot is being displayed + self.figure = Figure(figsize=(8, 6)) + self.init_ui() + + def init_ui(self): + # Main layout: vertical, dividing into two sections (top for controls/plot, bottom for output) + main_layout = QVBoxLayout() + + self.setStyleSheet(""" + QWidget { + background-color: #2e2e2e; /* Dark background */ + color: #ffffff; /* Light text */ + } + QLabel { + font-weight: bold; + color: #ffffff; /* Light text */ + } + QComboBox { + background-color: #444444; /* Dark combo box */ + color: #ffffff; /* Light text */ + border: 1px solid #888888; /* Light border */ + min-width: 200px; /* Set a minimum width */ + } + QLineEdit { + background-color: #444444; /* Dark input field */ + color: #ffffff; /* Light text */ + border: 1px solid #888888; /* Light border */ + padding: 5px; /* Add padding */ + fixed-height: 30px; /* Set a fixed height */ + min-width: 200px; /* Fixed width */ + } + QPushButton { + background-color: #4caf50; /* Green button */ + color: white; /* White text */ + border: none; /* No border */ + padding: 10px; /* Add padding */ + border-radius: 5px; /* Rounded corners */ + } + QPushButton:hover { + background-color: #45a049; /* Darker green on hover */ + } + QTextEdit { + background-color: #1e1e1e; /* Dark output box */ + color: #ffffff; /* Light text */ + border: 1px solid #888888; /* Light border */ + padding: 5px; /* Add padding */ + fixed-height: 150px; /* Set a fixed height */ + min-width: 800px; /* Set a minimum width to match the canvas */ + } + QTextEdit:focus { + border: 1px solid #00ff00; /* Green border on focus for visibility */ + } + QPushButton { + background-color: #4caf50; /* Green button */ + color: white; /* White text */ + border: none; /* No border */ + padding: 10px; /* Add padding */ + border-radius: 5px; /* Rounded corners */ + } + QPushButton:disabled { + background-color: rgba(76, 175, 80, 0.5); /* Transparent green when disabled */ + color: rgba(255, 255, 255, 0.5); /* Light text when disabled */ + } + QPushButton:hover { + background-color: #45a049; /* Darker green on hover */ + } + """) + + + # Horizontal layout for test options (left) and plot canvas (right) + top_layout = QHBoxLayout() + + # Create a QGroupBox for controls + controls_box = QGroupBox("Test Controls", self) + controls_box.setFixedHeight(600) + controls_layout = QVBoxLayout() # Layout for the controls + + # Dropdown for selecting the test + self.label = QLabel('Select Test:', self) + self.label.setFixedSize(100, 20) + controls_layout.addWidget(self.label) + + self.test_selector = QComboBox(self) + self.test_selector.addItems([ + "Linearity Test", "Deadtime Test", "Pedestal Test", + "Pixel Time Uncertainty Test", "Time Uncertainty Between Couples of Pixels" + ]) + self.test_selector.setFixedWidth(400) # Fixed width for the dropdown + self.test_selector.currentIndexChanged.connect(self.update_parameters) + controls_layout.addWidget(self.test_selector) + + # Container for parameter fields + self.param_widgets = QWidget(self) + self.param_widgets.setFixedSize(400, 400) + self.param_layout = QVBoxLayout(self.param_widgets) + controls_layout.addWidget(self.param_widgets) + + # Button to run the test + self.run_button = QPushButton('Run Test', self) + self.run_button.clicked.connect(self.run_test) + controls_layout.addWidget(self.run_button) + + # Set the controls layout to the group box + controls_box.setLayout(controls_layout) + + # Add the controls box to the top layout (left side) + top_layout.addWidget(controls_box) + + # Add a stretchable spacer to push the canvas to the right + top_layout.addSpacerItem(QSpacerItem(40, 20, QSizePolicy.Expanding, QSizePolicy.Minimum)) + + # Create a vertical layout for the plot container + self.plot_container = QWidget(self) + self.plot_layout = QVBoxLayout(self.plot_container) + + # Fixed size for the container + self.plot_container.setFixedSize(800, 650) # Set desired fixed size for the container + + # Create a vertical layout for the plot (canvas and toolbar) + self.canvas = FigureCanvas(self.figure) + self.canvas.setFixedSize(800, 600) # Fixed size for the canvas + + # Add toolbar for zooming and panning + self.toolbar = NavigationToolbar(self.canvas, self) + self.toolbar.setFixedHeight(50) # Fixed height for the toolbar + + # Add the toolbar and canvas to the plot layout + self.plot_layout.addWidget(self.toolbar) # Toolbar stays on top + self.plot_layout.addWidget(self.canvas) # Canvas below toolbar + + # Add the plot container to the top layout (to the right) + top_layout.addWidget(self.plot_container) + + # Add the top layout (controls + canvas) to the main layout + main_layout.addLayout(top_layout) + + # Navigation buttons + nav_layout = QHBoxLayout() + self.prev_button = QPushButton('Previous Plot', self) + self.prev_button.clicked.connect(self.show_previous_plot) + self.prev_button.setEnabled(False) # Initially disabled + nav_layout.addWidget(self.prev_button) + + self.next_button = QPushButton('Next Plot', self) + self.next_button.clicked.connect(self.show_next_plot) + self.next_button.setEnabled(False) # Initially disabled + nav_layout.addWidget(self.next_button) + + main_layout.addLayout(nav_layout) + + # Output text box (bottom section of the main layout) + self.output_text_edit = QTextEdit(self) + self.output_text_edit.setReadOnly(True) # To prevent user editing + self.output_text_edit.setFixedHeight(150) # Set a fixed height for the output box + self.output_text_edit.setMinimumWidth(800) # Set a minimum width to match the canvas + main_layout.addWidget(self.output_text_edit) + + # Set the main layout to the window + self.setLayout(main_layout) + self.setWindowTitle('Test Runner GUI') + self.showFullScreen() + + + + def get_parameters_from_module(self, module): + # Fetch parameters from the module + if hasattr(module, 'get_args'): + parser = module.get_args() + return {arg.dest: arg.default for arg in parser._actions if isinstance(arg, argparse._StoreAction)} + else: + raise RuntimeError('No get_args function found in module.') + + def debug_layout(self): + for i in range(self.param_layout.count()): + item = self.param_layout.itemAt(i) + widget = item.widget() + if widget: + print(f"Widget in layout: {widget.objectName()}") + + def update_parameters(self): + # Clear existing parameter fields + for i in reversed(range(self.param_layout.count())): + widget = self.param_layout.itemAt(i).widget() + if widget: + widget.deleteLater() + + # Get the selected test and corresponding module + selected_test = self.test_selector.currentText() + test_modules = { + "Linearity Test": linearity_test, + "Deadtime Test": deadtime_test, + "Pedestal Test": pedestal_test, + "Pixel Time Uncertainty Test": pix_tim_uncertainty_test, + "Time Uncertainty Between Couples of Pixels": pix_couple_tim_uncertainty_test + } + + module = test_modules.get(selected_test) + if module: + try: + self.params = self.get_parameters_from_module(module) + for param, default in self.params.items(): + if param=="temp_output": + continue + label = QLabel(f'{param}:', self) + self.param_layout.addWidget(label) + entry = QLineEdit(self) + entry.setText(str(default).replace('[','').replace(']','').replace(',','')) + entry.setObjectName(param) + entry.setFixedWidth(400) # Set fixed width for QLineEdit + self.param_layout.addWidget(entry) + + # Force update the layout + self.param_widgets.update() + QTimer.singleShot(0, self.param_widgets.update) # Ensures the layout is updated + except Exception as e: + QMessageBox.critical(self, "Error", f"Failed to fetch parameters: {e}") + else: + QMessageBox.critical(self, "Error", "No test selected or test not found") + + def run_test(self): + # Clean up old plot files to avoid loading leftover files + self.cleanup_tempdir() + + selected_test = self.test_selector.currentText() + test_modules = { + "Linearity Test": linearity_test, + "Deadtime Test": deadtime_test, + "Pedestal Test": pedestal_test, + "Pixel Time Uncertainty Test": pix_tim_uncertainty_test, + "Time Uncertainty Between Couples of Pixels": pix_couple_tim_uncertainty_test + } + module = test_modules.get(selected_test) + + if module: + params = [] + self.update() + self.repaint() + + # Generate temporary output path + self.temp_output = tempfile.gettempdir() + print(f"Temporary output dir: {self.temp_output}") # Debug print + + for param, default in self.params.items(): + widget_list = self.param_widgets.findChildren(QLineEdit, param) + if widget_list: + widget = widget_list[0] + params.append(f"--{param}") + params.extend(widget.text().split(' ')) + if param == "output": + params.append(f"--output={widget.text()}") + + params.append(f"--temp_output={self.temp_output}") + else: + print(f"Widget with name {param} not found") + + test_script_path = os.path.abspath(module.__file__) + command = [sys.executable, test_script_path] + params + print(f"Command to run: {command}") # Debug print + + try: + self.output_text_edit.clear() + + self.process = QProcess(self) + self.process.setProcessChannelMode(QProcess.MergedChannels) + self.process.readyReadStandardOutput.connect(self.read_process_output) + self.process.finished.connect(self.process_finished) + + QTimer.singleShot(0, lambda: self.process.start(sys.executable, [test_script_path] + params)) + + except Exception as e: + QMessageBox.critical(self, "Error", f"Failed to run the test: {e}") + else: + QMessageBox.critical(self, "Error", "No parameters found for the selected test") + + self.plot_files = [os.path.join(self.temp_output, f'plot{i}.pkl') for i in range(1, 3)] + + + def read_process_output(self): + """Reads and displays the process output in real-time.""" + if self.process: + output = self.process.readAllStandardOutput().data().decode("utf-8").strip() + if output: + self.output_text_edit.append(output) + + def process_finished(self): + """Handle the process when it finishes.""" + if self.process.exitCode() == 0: + QMessageBox.information(self, "Test Output", "Test completed successfully.") + + # Delay to ensure file creation is complete + QTimer.singleShot(1000, self.check_and_display_plot) + else: + QMessageBox.critical(self, "Error", f"Test failed with exit code {self.process.exitCode()}") + + def check_and_display_plot(self): + plot_files = [os.path.join(self.temp_output, f'plot{i}.pkl') for i in range(1, 3)] + self.display_plot([f for f in plot_files if os.path.exists(f)]) + + def display_plot(self, plot_files): + """Loads the plots from the pickle files and displays them on the canvas.""" + self.plots = [] + self.current_plot_index = 0 + # Load all available plots from the pickle files + for plot_file in plot_files: + with open(plot_file, 'rb') as f: + self.plots.append(pickle.load(f)) # Load the plot data + # Display the first plot if there are any loaded plots + if self.plots: + self.update_plot_canvas() + + # Enable the "Next" button if there is more than one plot + if len(self.plots) > 1: + self.next_button.setEnabled(True) + else: + self.next_button.setEnabled(False) + + + + def update_plot_canvas(self): + """Updates the canvas with the current plot.""" + if not self.plots: + return + + try: + # Load the current figure + + with open(self.plot_files[self.current_plot_index], 'rb') as f: + loaded_figure = pickle.load(f) + #loaded_figure = self.plots[self.current_plot_index] + + # Remove the old canvas and toolbar from the plot layout + self.plot_layout.removeWidget(self.canvas) + self.canvas.deleteLater() # Properly delete the old canvas + self.plot_layout.removeWidget(self.toolbar) + self.toolbar.deleteLater() # Properly delete the old toolbar + + # Create a new canvas with the loaded figure + self.canvas = FigureCanvas(loaded_figure) + self.canvas.setFixedSize(800, 600) # Set a fixed size for the canvas + + # Adjust the plot margins to ensure the x-axis is visible + loaded_figure.subplots_adjust(bottom=0.15) # Increase the bottom margin + #loaded_figure.tight_layout(pad=2.0) # Use tight layout with padding + + self.canvas.draw() + + # Create a new toolbar with the loaded figure + self.toolbar = NavigationToolbar(self.canvas, self) + self.toolbar.setFixedHeight(50) + + # Clear the plot layout and re-add toolbar above the canvas + self.plot_layout.addWidget(self.toolbar) + self.plot_layout.addWidget(self.canvas) + + + except Exception as e: + QMessageBox.critical(self, "Error", f"Failed to load plot: {e}") + + + + + + + + + + def show_next_plot(self): + if self.current_plot_index < len(self.plots) - 1: + self.current_plot_index += 1 + self.update_plot_canvas() + self.prev_button.setEnabled(True) + if self.current_plot_index == len(self.plots) - 1: + self.next_button.setEnabled(False) + else: + self.next_button.setEnabled(False) + + def show_previous_plot(self): + if self.current_plot_index > 0: + self.current_plot_index -= 1 + self.update_plot_canvas() + self.next_button.setEnabled(True) + if self.current_plot_index == 0: + self.prev_button.setEnabled(False) + else: + self.prev_button.setEnabled(False) + + def cleanup_tempdir(self): + """Remove old plot files in temp directory.""" + for i in range(1, 3): + plot_file = os.path.join(tempfile.gettempdir(), f'plot{i}.pkl') + if os.path.exists(plot_file): + os.remove(plot_file) + + +if __name__ == '__main__': + app = QApplication(sys.argv) + ex = TestRunner() + sys.exit(app.exec_()) diff --git a/src/nectarchain/user_scripts/dmousadi/test_scripts/requirements.txt b/src/nectarchain/user_scripts/dmousadi/test_scripts/requirements.txt new file mode 100644 index 00000000..e82bfc4e --- /dev/null +++ b/src/nectarchain/user_scripts/dmousadi/test_scripts/requirements.txt @@ -0,0 +1,16 @@ +astropy==5.3.4 +ctapipe==0.19.1 +ctapipe_io_nectarcam==0.1.5 +h5py==3.10.0 +hdf5plugin==5.0.0 +iminuit==2.30.0 +lmfit==1.2.2 +matplotlib==3.8.2 +nectarchain==0.1.9 +numpy +pandas==2.2.3 +PyQt5==5.15.11 +PyQt5_sip==12.15.0 +scipy==1.14.1 +traitlets==5.14.3 + diff --git a/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/__init__.py b/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/deadtime_test.py b/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/deadtime_test.py new file mode 100644 index 00000000..6c0d1230 --- /dev/null +++ b/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/deadtime_test.py @@ -0,0 +1,559 @@ +#don't forget to set environment variable NECTARCAMDATA + +import numpy as np +import pathlib +import os +import sys + +import matplotlib.pyplot as plt + +from utils import source_ids_deadtime,deadtime_labels, pois, deadtime_and_expo_fit, err_ratio + +import argparse +from iminuit import Minuit + +import matplotlib.pyplot as plt +import numpy as np +from test_tools_components import DeadtimeTestTool +from astropy import units as u +from utils import ExponentialFitter +from iminuit import Minuit +import pickle + +def get_args(): + """ + Parses command-line arguments for the deadtime test script. + + Returns: + 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 give a list of runs (run numbers with spaces inbetween), a corresponding source list and an output directory to save the final plot.\n' + +'If the data is not in NECTARCAMDATA, the files will be downloaded through DIRAC.\n For the purposes of testing this script, default data is from the runs used for this test in the TRR document.\n' + +'You can optionally specify the number of events to be processed (default 1000).\n') + parser.add_argument('-r','--runlist', type=int,nargs='+', help='List of runs (numbers separated by space)', required=False, default = [i for i in range(3332,3350)]+[i for i in range(3552,3562)]) + parser.add_argument('-s','--source', type=int, choices = [0,1,2], nargs='+', help='List of corresponding source for each run: 0 for random generator, 1 for nsb source, 2 for laser', required=False , default = source_ids_deadtime) + parser.add_argument('-e','--evts', type = int, help='Number of events to process from each run. Default is 1000', required=False, default=1000) + parser.add_argument('-o','--output', type=str, help='Output directory. If none, plot will be saved in current directory', required=False, default='./') + parser.add_argument("--temp_output", help="Temporary output directory for GUI", default=None) + + + return parser + + + + + +def main(): + """ + Runs the deadtime test script, which performs deadtime tests B-TEL-1260 and B-TEL-1270. + + The script takes command-line arguments to specify the list of runs, corresponding sources, number of events to process, and output directory. It then processes the data for each run, performs an exponential fit to the deadtime distribution, and generates two plots: + + 1. A plot of deadtime percentage vs. collected trigger rate, with the CTA requirement indicated. + 2. 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 them to a temporary output directory for use in a GUI. + """ + + parser = get_args() + args = parser.parse_args() + + runlist = args.runlist + ids = args.source + + nevents = args.evts + + output_dir = os.path.abspath(args.output) + temp_output = os.path.abspath(args.temp_output) if args.temp_output else None + + print(f"Output directory: {output_dir}") # Debug print + print(f"Temporary output file: {temp_output}") # Debug print + + + sys.argv = sys.argv[:1] + + # ucts_timestamps = np.zeros((len(runlist),nevents)) + # delta_t = np.zeros((len(runlist),nevents-1)) + # event_counter = np.zeros((len(runlist),nevents)) + # busy_counter = np.zeros((len(runlist),nevents)) + # collected_triger_rate = np.zeros(len(runlist)) + # time_tot = np.zeros(len(runlist)) + # deadtime_us=np.zeros((len(runlist),nevents-1)) + # deadtime_pc = np.zeros(len(runlist)) + + ucts_timestamps = [] + delta_t = [] + event_counter = [] + busy_counter = [] + collected_trigger_rate = [] + time_tot = [] + deadtime_us = [] + deadtime_pc = [] + + labels = deadtime_labels + + + for i,run in enumerate(runlist): + + print("PROCESSING RUN {}".format(run)) + tool = DeadtimeTestTool( + progress_bar=True, run_number=run, max_events=nevents, events_per_slice = 10000, log_level=20, peak_height=10, window_width=16, overwrite=True + ) + tool.initialize() + tool.setup() + tool.start() + output=tool.finish() + ucts_timestamps.append(output[0]) + delta_t.append(output[1]) + event_counter.append(output[2]) + busy_counter.append(output[3]) + collected_trigger_rate.append(output[4].value) + time_tot.append(output[5].value) + deadtime_pc.append(output[6]) + + deadtime_us.append((delta_t[i]*u.ns).to(u.us)) + + + collected_trigger_rate=np.array(collected_trigger_rate) + deadtime_pc = np.array(deadtime_pc) + + + parameter_A2_new_list = [] + parameter_lambda_new_list = [] + parameter_tau_new_list = [] + parameter_A2_err_new_list = [] + parameter_lambda_err_new_list = [] + parameter_tau_err_new_list = [] + parameter_R2_new_list = [] + + + fitted_rate = [] + for i in range(len(runlist)): + print("fitting rate for run", runlist[i]) + + dt_mus = deadtime_us[i].value + dt_mus = dt_mus[dt_mus>0] + + + lim_sup_mus = 120 + lim_sup_s = lim_sup_mus*1e-6 + nr_bins=500 + + rate_initial_guess = 40000 + + data_content, bin_edges = np.histogram(dt_mus*1e-6,bins=np.linspace(0.001e-6,lim_sup_s,nr_bins)) + + init_param = [np.sum(data_content),0.6e-6,1./rate_initial_guess] + + + fitter = ExponentialFitter(data_content,bin_edges=bin_edges) + m = Minuit( fitter.compute_minus2loglike, init_param, name=('Norm','deadtime','1/Rate') ) + + # Set Parameter Limits and tolerance + + m.errors["Norm"] = 0.3*init_param[0] + m.limits["Norm"] = (0.,None) + + m.errors["deadtime"] = 0.1e-6 + m.limits["deadtime"] = ( 0.6e-6,1.1e-6) # Put some tigh constrain as the fit will be in trouble when it expect 0. and measured something instead. + + m.print_level = 2 + + + + res = m.migrad(2000000) + + + + fitted_params = np.array( [res.params[p].value for p in res.parameters] ) + #print(fitted_params) + + fitted_params_err = np.array( [res.params[p].error for p in res.parameters] ) + #print(fitted_params_err) + + print(f"Dead-Time is {1.e6*fitted_params[1]:.3f} +- {1.e6*fitted_params_err[1]:.3f} µs") + print(f"Rate is {1./fitted_params[2]:.2f} +- {fitted_params_err[2]/(fitted_params[2]**2):.2f} Hz") + print(f"Expected run duration is {fitted_params[0]*fitted_params[2]:.2f} s") + + fitted_rate.append(1./fitted_params[2]) + + # plt.savefig(figurepath + 'deadtime_exponential_fit_nsb_run{}_newfit_cutoff.png'.format(run)) + + + y = data_content + y_fit = fitter.expected_distribution(fitted_params) + # residual sum of squares + ss_res = np.sum((y - y_fit) ** 2) + + # total sum of squares + ss_tot = np.sum((y - np.mean(y)) ** 2) + + # r-squared + r2 = 1 - (ss_res / ss_tot) + #print(r2) + + parameter_A2_new_list.append(fitted_params[0]) + parameter_lambda_new_list.append(1./fitted_params[2]/1e3) #kHz + parameter_tau_new_list.append(1.e6*fitted_params[1]) #musec + parameter_A2_err_new_list.append(fitted_params_err[0]) + parameter_lambda_err_new_list.append(fitted_params_err[2]/(fitted_params[2]**2)/1e3) + parameter_tau_err_new_list.append(1.e6*fitted_params_err[1]) + + parameter_R2_new_list.append(r2) + + + + + deadtime_from_fit = (parameter_tau_new_list) + deadtime_from_fit_err = (parameter_tau_err_new_list) + lambda_from_fit = (parameter_lambda_new_list) + lambda_from_fit_err = (parameter_lambda_err_new_list) + A2_from_fit = (parameter_A2_new_list) + A2_from_fit_err = (parameter_A2_err_new_list) + R2_from_fit = (parameter_R2_new_list) + + + + ####################################### + #PLOT + # print(event_counter) + # print(busy_counter) + # print(collected_triger_rate) + # print(deadtime_pc) + # print(busy_counter[:,-1]) + # print(event_counter[:,-1]) + # print(busy_counter[:,-1]/(event_counter[:,-1]+busy_counter[:,-1])) + plt.clf() + fig, ax = plt.subplots(figsize=(10*1.1,10*1.1/1.61)) # constrained_layout=True) + ids = np.array(ids) + runlist = np.array(runlist) + + ratio_list = [] + collected_rate = [] + err = [] + + for source in range(0,3): + + # runl = np.where(ids==source)[0] + # for i in runl: + # #print(labels[ids[i]]) + + # deadtime = (deadtime_from_fit[i]*1e3)*u.ns + # delta_deadtime = deadtime_from_fit_err[i]*1e3 + # freq = (1/deadtime).to(u.kHz).value + # freq_err = delta_deadtime/deadtime.value**2 + + # rate = lambda_from_fit[i] + # rate_err=lambda_from_fit_err[i] + # ratio = rate/freq + # ratio_list.append(ratio*100) + # collected_rate.append(collected_triger_rate[i]) + # ratio_err = np.sqrt((rate_err/freq)**2 + (freq_err*rate/(freq**2))) + + + # err.append(ratio_err*100) + + Y = list(np.array(collected_trigger_rate[np.where(ids==source)[0]])/1000) + X = list(np.array(deadtime_pc[np.where(ids==source)[0]])) + #err = list(err) + X_sorted = [x for y, x in sorted(zip(Y, X))] + #err_sorted = [err for y,err in sorted(zip(Y,err))] + + plt.plot(sorted(Y), X_sorted, #yerr = err_sorted, + alpha=0.6, ls='-', marker='o',color=labels[source]['color'], label = labels[source]['source']) + + + + plt.xlabel('Collected Trigger Rate [kHz]') + plt.ylabel(r'Deadtime [%]') + + + plt.axhline(5, ls='-', color='gray', alpha=0.4) + plt.axvline(7, ls='-', color='gray', alpha=0.4,) + + ax.text(28, 6.25, 'CTA requirement', color='gray', fontsize=20, alpha=0.6, + horizontalalignment='center', + verticalalignment='center') #transform=ax.transAxes) + + plt.legend() + + plt.xlim(-0.5,16) + # plt.grid(which='both') + plt.yscale(u'log') + plt.ylim(1e-2,1e2) + plt.savefig(os.path.join(output_dir,"deadtime.png")) + + if temp_output: + with open(os.path.join(args.temp_output, 'plot1.pkl'), 'wb') as f: + pickle.dump(fig, f) + + + + + + + + ################################################## + #SECOND PLOT + plt.clf() + plt.figure(figsize=(10,10/1.61)) + fig, ((ax1, ax2)) = plt.subplots(2, 1, sharex='col', sharey='row', figsize=(10,8), + # sharex=True, sharey=True, + gridspec_kw={'height_ratios': [5,2]}) + + + + x = collected_trigger_rate/1000 + rate=np.array(lambda_from_fit) + y = lambda_from_fit + rate_err=np.array(lambda_from_fit_err) + relative = (y-x)/x * 100 + #print(np.argmin(relative)) + + x_err = 0 + err_ratio = relative * ( ((rate_err + x_err)/(y - x)) + x_err/x) + ax2.errorbar(x, relative, + xerr= x_err, yerr=err_ratio, + alpha=0.9, ls=' ', marker='o', + color='C1') + ax2.set_ylim(-25,25) + + x=range(0,60) + + ax1.set_xscale(u'log') + ax1.set_yscale(u'log') + ax1.plot(x, x, color='gray', ls='--', alpha=0.5) + + ax2.plot(x, np.zeros(len(x)), color='gray', ls='--', alpha=0.5) + ax2.fill_between(x, np.ones(len(x))*(-10), np.ones(len(x))*(10), color='gray', alpha=0.1) + + ax2.set_xlabel('Collected Trigger Rate [kHz]') + ax1.set_ylabel(r'Rate from fit [kHz]') + ax2.set_ylabel(r'Relative difference [%]') + + + ax1.set_xlim(1,60) + ax1.set_ylim(1,60) + ax2.set_xlim(1,60) + + ids = np.array(ids) + runlist = np.array(runlist) + #print(lambda_from_fit) + #print("coll",collected_triger_rate) + for source in range(0,3): + runl = np.where(ids==source)[0] + + + + # print(collected_triger_rate[runl]) + # print(rate[runl]) + ax1.errorbar(collected_trigger_rate[runl]/1000, + rate[runl], + #xerr=((df_mean_nsb[df_mean_nsb['Run']==run]['Collected_trigger_rate[Hz]_err']))/1000, + yerr=rate_err[runl], + alpha=0.9, + ls=' ', marker='o', color=labels[source]['color'], label = labels[source]['source']) + # label = 'Run {} ({} V)'.format(run, df_mean_rg[df_mean_rg['Run']==run]['Voltage[V]'].values[0])) + + ax1.legend(frameon=False, prop={'size':10}, + loc="upper left", ncol=1) + + + plt.savefig(os.path.join(output_dir,"deadtime_meas.png")) + + if temp_output: + with open(os.path.join(args.temp_output, 'plot2.pkl'), 'wb') as f: + pickle.dump(fig, f) + + plt.close('all') + + + + +if __name__ == "__main__": + main() + + + + + + +##################################PREVIOUS############################### +# collected_rate = [] + + +# parameter_A_new_list = [] +# parameter_R_new_list = [] +# parameter_A_err_new_list = [] +# parameter_R_err_new_list = [] +# first_bin_length_list = [] +# tot_nr_events_histo_list = [] +# ucts_busytime = [] +# ucts_event_counter = [] +# total_delta_t_for_busy_time_list = [] + +# deadtime = list() +# deadtime_bin = list() +# deadtime_err = list() +# deadtime_bin_length = list() + +# for i, run in enumerate(runlist): +# deadtime_run, deadtime_bin_run, deadtime_err_run, deadtime_bin_length_run, \ +# total_delta_t_for_busy_time, parameter_A_new, parameter_R_new, parameter_A_err_new, parameter_R_err_new, \ +# first_bin_length, tot_nr_events_histo = deadtime_and_expo_fit(time_tot[i],deadtime_us[i], run) +# total_delta_t_for_busy_time_list.append(total_delta_t_for_busy_time) +# parameter_A_new_list.append(parameter_A_new) +# parameter_R_new_list.append(parameter_R_new) +# parameter_A_err_new_list.append(parameter_A_err_new) +# parameter_R_err_new_list.append(parameter_R_err_new) + +# tot_nr_events_histo_list.append(first_bin_length) +# tot_nr_events_histo_list.append(tot_nr_events_histo) +# deadtime.append(deadtime_run) +# deadtime_bin.append(deadtime_bin_run) +# deadtime_err.append(deadtime_err_run) +# deadtime_bin_length.append(deadtime_bin_length_run/np.sqrt(12)) +# #print(run, parameter_A_new, parameter_R_new) + + +# deadtime_from_first_bin = np.array(deadtime_bin) +# deadtime_bin_length = np.array(deadtime_bin_length) +# rate = (np.array(parameter_R_new_list) * 1 / u.us).to(u.kHz).to_value() #in kHz +# rate_err = (np.array(parameter_R_err_new_list) * 1 / u.us).to(u.kHz).to_value() +# A_from_fit = (parameter_A_new_list) +# A_from_fit_err = (parameter_A_err_new_list) +# ucts_busy_rate = (np.array(busy_counter[:,-1]) / (np.array(time_tot) * u.s).to(u.s)).to( +# u.kHz).value +# nr_events_from_histo = (tot_nr_events_histo) +# first_bin_delta_t = first_bin_length + +# deadtime_average = np.average(deadtime_bin, +# weights=1 / (deadtime_bin_length ** 2)) +# deadtime_average_err_nsb = np.sqrt(1 / (np.sum(1 / deadtime_bin_length ** 2))) + + +# ####################################################################################### + + + + +# #B-TEL-1260 +# plt.clf() +# fig, ax = plt.subplots(figsize=(10*1.1,10*1.1/1.61)) +# ids = np.array(ids) +# runlist = np.array(runlist) + +# ratio_list = [] +# collected_rate = [] +# err = [] +# for source in range(0,3): +# runl = np.where(ids==source)[0] +# for i in runl: +# #print(labels[ids[i]]) + +# deadtime = (deadtime_bin[i]*1e3)*u.ns +# delta_deadtime = deadtime_bin_length[i]*1e3 +# freq = (1/deadtime).to(u.kHz).value +# freq_err = delta_deadtime/deadtime.value**2 + +# ratio = rate[i]/freq +# ratio_list.append(np.array(ratio)*100) +# ratio_err = np.sqrt((rate_err[i]/freq)**2 + (freq_err*rate[i]/(freq**2))) +# err.append(ratio_err*100) + + +# Y = list(np.array(collected_triger_rate[np.where(ids==source)[0]])/1000) +# X = list(ratio_list) +# err = list(err) +# X_sorted = [x for y, x in sorted(zip(Y, X))] +# err_sorted = [err for y,err in sorted(zip(Y,err))] + +# plt.errorbar(sorted(Y), X_sorted, yerr = err_sorted, alpha=0.6, ls='-', marker='o',color=labels[source]['color'], label = labels[source]['source']) + +# plt.xlabel('Collected Trigger Rate [kHz]') +# plt.ylabel(r'Deadtime [%]') + + +# plt.axhline(5, ls='-', color='gray', alpha=0.4) +# plt.axvline(7, ls='-', color='gray', alpha=0.4,) + +# ax.text(28, 6.25, 'CTA requirement', color='gray', fontsize=20, alpha=0.6, +# horizontalalignment='center', +# verticalalignment='center') +# plt.legend() + +# plt.xlim(-0.5,33) +# # plt.grid(which='both') +# plt.yscale(u'log') +# plt.ylim(1e-2,1e2) +# plt.savefig(os.path.join(output_dir,"deadtime.png")) + + +# ############################################################################# + + + + + + + + + +# #B-TEL-1670 +# plt.clf() +# f, ((ax1, ax2)) = plt.subplots(2, 1, sharex='col', sharey='row', figsize=(10,8), +# gridspec_kw={'height_ratios': [5,2]}) + + +# #plt.figure(figsize=(10,10/1.61)) +# x = collected_triger_rate/1000 +# y = rate +# relative = (y-x)/x * 100 +# # print(np.argmin(relative)) +# # print(collected_triger_rate[5]) +# # print(rate[5]) +# x_err = 0 +# err_ratio = relative * ( ((rate_err + x_err)/(y - x)) + x_err/x) +# ax2.errorbar(x, relative, +# xerr= x_err, yerr=err_ratio, +# alpha=0.9, ls=' ', marker='o', +# color='C1') +# ax2.set_ylim(-25,25) + +# x=range(0,60) + +# ax1.set_xscale(u'log') +# ax1.set_yscale(u'log') +# ax1.plot(x, x, color='gray', ls='--', alpha=0.5) + +# ax2.plot(x, np.zeros(len(x)), color='gray', ls='--', alpha=0.5) +# ax2.fill_between(x, np.ones(len(x))*(-10), np.ones(len(x))*(10), color='gray', alpha=0.1) + +# ax2.set_xlabel('Collected Trigger Rate [kHz]') +# ax1.set_ylabel(r'Rate from fit [kHz]') +# ax2.set_ylabel(r'Relative difference [%]') + +# ax1.set_xlim(1,60) +# ax1.set_ylim(1,60) +# ax2.set_xlim(1,60) + +# ids = np.array(ids) +# runlist = np.array(runlist) +# for source in range(0,3): +# runl = np.where(ids==source)[0] +# #print(collected_triger_rate[runl]) +# ax1.errorbar(collected_triger_rate[runl]/1000, +# rate[runl], +# #xerr=((df_mean_nsb[df_mean_nsb['Run']==run]['Collected_trigger_rate[Hz]_err']))/1000, +# yerr=rate_err[runl], +# alpha=0.9, +# ls=' ', marker='o', color=labels[source]['color'], label = labels[source]['source']) +# # label = 'Run {} ({} V)'.format(run, df_mean_rg[df_mean_rg['Run']==run]['Voltage[V]'].values[0])) + +# ax1.legend(frameon=False, prop={'size':10}, +# loc="upper left", ncol=1) + + +# plt.savefig(os.path.join(output_dir,"deadtime_meas.png")) diff --git a/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/linearity_test.py b/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/linearity_test.py new file mode 100644 index 00000000..ad32295e --- /dev/null +++ b/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/linearity_test.py @@ -0,0 +1,292 @@ +#don't forget to set environment variable NECTARCAMDATA + +import numpy as np +import pathlib +import os +import sys +import matplotlib.pyplot as plt + +from utils import filters, adc_to_pe, optical_density_390ns, trasmission_390ns, fit_function, err_ratio, err_sum, plot_parameters +from lmfit.models import Model +from test_tools_components import LinearityTestTool +import argparse +import pickle + + +def get_args(): + """ + Parses command-line arguments for the linearity test script. + + Returns: + argparse.ArgumentParser: The parsed command-line arguments. + """ + + parser = argparse.ArgumentParser(description='Linearity test B-TEL-1390 & Intensity resolution B-TEL-1010. \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 give a list of runs (run numbers with spaces inbetween), a corresponding transmission list and an output directory to save the final plot.\n' + +'If the data is not in NECTARCAMDATA, the files will be downloaded through DIRAC.\n For the purposes of testing this script, default data is from the runs used for this test in the TRR document.\n' + +'You can optionally specify the number of events to be processed (default 500) and the number of pixels used (default 70).\n') + parser.add_argument('-r','--runlist', type=int,nargs='+', help='List of runs (numbers separated by space)', required=False, default = [i for i in range (3404,3424)]+[i for i in range(3435,3444)]) + parser.add_argument('-t','--trans', type=float,nargs='+', help='List of corresponding transmission for each run', required=False , default = trasmission_390ns) + parser.add_argument('-e','--evts', type = int, help='Number of events to process from each run. Default is 500', required=False, default=500) + parser.add_argument('-o','--output', type=str, help='Output directory. If none, plot will be saved in current directory', required=False, default='./') + parser.add_argument("--temp_output", help="Temporary output directory for GUI", default=None) + + return parser + + + +def main(): + """ + The `main()` function is the entry point of the linearity test script. It parses the command-line arguments, processes the specified runs, and generates plots to visualize the linearity and charge resolution of the detector. The function performs the following key steps: + + 1. Parses the command-line arguments using the `get_args()` function, which sets up the argument parser and handles the input parameters. + 2. Iterates through the specified run list, processing each run using the `LinearityTestTool` class. This tool initializes, sets up, starts, and finishes the processing for each run, returning the relevant output data. + 3. Normalizes the high-gain and low-gain charge values using the charge value at 0.01 transmission. + 4. Generates three subplots: + - The first subplot shows the estimated charge vs. the true charge, with the fitted linear function for both high-gain and low-gain channels. + - The second subplot shows the residuals between the estimated and true charge, as a percentage. + - The third subplot shows the ratio of high-gain to low-gain charge, with a fitted linear function. + 5. Saves the generated plots to the specified output directory, and optionally saves temporary plot files for a GUI. + 6. Generates an additional plot to visualize the charge resolution, including the statistical limit. + 7. Saves the charge resolution plot to the specified output directory, and optionally saves a temporary plot file for a GUI. + """ + parser = get_args() + args = parser.parse_args() + + + + + runlist = args.runlist + transmission=args.trans #corresponding transmission for above data + + nevents = args.evts + + output_dir = os.path.abspath(args.output) + temp_output= os.path.abspath(args.temp_output) if args.temp_output else None + + print(f"Output directory: {output_dir}") # Debug print + print(f"Temporary output file: {temp_output}") # Debug print + + + sys.argv = sys.argv[:1] + + #runlist = [3441] + + charge = np.zeros((len(runlist),2)) + std = np.zeros((len(runlist),2)) + std_err = np.zeros((len(runlist),2)) + + index = 0 + for run in runlist: + print("PROCESSING RUN {}".format(run)) + tool = LinearityTestTool( + progress_bar=True, run_number=run, events_per_slice = 999, max_events=nevents, log_level=20, window_width=14, overwrite=True + + ) + tool.initialize() + tool.setup() + tool.start() + output = tool.finish() + + charge[index], std[index], std_err[index], npixels = output + index += 1 + + + #print("FINAL",charge) + + #we assume that they overlap at 0.01 so they should have the same value + #normalise high gain and low gain using charge value at 0.01 + transmission=np.array(transmission) + norm_factor_hg = charge[np.argwhere((transmission<1.1e-2) & (transmission>9e-3)),0][0] + #print(norm_factor_hg) + norm_factor_lg = charge[np.argwhere((transmission<1.1e-2) & (transmission>9e-3)),1][0] + #print(norm_factor_lg) + charge_norm_hg = charge[:,0]/norm_factor_hg + charge_norm_lg = charge[:,1]/norm_factor_lg + std_norm_hg = std[:,0]/norm_factor_hg + std_norm_lg = std[:,1]/norm_factor_lg + + #true charge is transmission as percentage of hg charge at 0.01 + true_charge = transmission/0.01 * norm_factor_hg + + + + fig, axs = plt.subplots(3, 1, sharex='col', sharey='row', figsize=(10,11), gridspec_kw={'height_ratios': [4,2,2]}) + axs[0].grid(True, which="both") + axs[0].set_ylabel("Estimated charge [p.e.]") + axs[0].set_yscale('log') + axs[0].set_xscale('log') + axs[0].axvspan(10,1000,alpha=0.2,color="orange") + axs[1].grid(True, which="both") + axs[1].set_ylabel("Residuals [%]") + axs[1].set_xscale('log') + axs[1].set_ylim((-100,100)) + axs[1].axvspan(10,1000,alpha=0.2,color="orange") + axs[2].grid(True, which="both") + axs[2].set_ylabel("HG/LG ratio") + axs[2].set_yscale('log') + axs[2].set_xscale('log') + axs[2].set_ylim((0.5,20)) + axs[2].axvspan(10,1000,alpha=0.2,color="orange") + axs[2].set_xlabel("Illumination charge [p.e.]") + + + + + + + for channel, (channel_charge, channel_std, name) in enumerate(zip([charge_norm_hg,charge_norm_lg],[std_norm_hg,std_norm_lg],["High Gain","Low Gain"])): + yx = zip(true_charge,channel_charge,channel_std,runlist) #sort by true charge + + ch_sorted = np.array(sorted(yx)) + #print(ch_sorted) + + #linearity + model = Model(fit_function) + params = model.make_params(a=100, b=0) + true = ch_sorted[:,0] + + ch_charge = ch_sorted[:,1] * norm_factor_hg[0] + ch_std = ch_sorted[:,2] * norm_factor_hg[0] + ch_err = ch_std/np.sqrt(npixels) + + + ch_fit = model.fit(ch_charge[plot_parameters[name]["linearity_range"][0]:plot_parameters[name]["linearity_range"][1]], params, + weights=1/ch_err[plot_parameters[name]["linearity_range"][0]:plot_parameters[name]["linearity_range"][1]], + x=true[plot_parameters[name]["linearity_range"][0]:plot_parameters[name]["linearity_range"][1]]) + + #print(ch_fit.fit_report()) + + + a = ch_fit.params['a'].value + b = ch_fit.params['b'].value + a_err = ch_fit.params['a'].stderr + b_err = ch_fit.params['b'].stderr + + + + + axs[0].errorbar(true,ch_charge,yerr = ch_err, label = name, ls='',marker='o',color = plot_parameters[name]["color"]) + axs[0].plot(true[plot_parameters[name]["linearity_range"][0]:plot_parameters[name]["linearity_range"][1]], + fit_function(true[plot_parameters[name]["linearity_range"][0]:plot_parameters[name]["linearity_range"][1]],a,b),color = plot_parameters[name]["color"]) + + + + + axs[0].text(plot_parameters[name]["label_coords"][0], plot_parameters[name]["label_coords"][0], name, + va='top', fontsize=20, color=plot_parameters[name]["color"], alpha=0.9) + axs[0].text(plot_parameters[name]["text_coords"][0], plot_parameters[name]["text_coords"][1], 'y = ax + b\n' + + r'a$_{\rm %s}=%2.2f \pm %2.2f$ '%(plot_parameters[name]["initials"],round(a,2), + round(a_err,2)) + + '\n' + r'b$_{\rm %s}=(%2.2f \pm %2.2f) $ p.e.' %(plot_parameters[name]["initials"],round(b,2), + round(b_err,2)) + + "\n" r"$\chi_{\mathrm{%s}}^2/{\mathrm{dof}} =%2.1f/%d$" %(plot_parameters[name]["initials"],round(ch_fit.chisqr,1), + int(ch_fit.chisqr/ch_fit.redchi)), + backgroundcolor='white', bbox=dict(facecolor='white', edgecolor=plot_parameters[name]["color"], lw=3, + boxstyle='round,pad=0.5'), + ha='left', va='top', fontsize=15, color='k', alpha=0.9) + + #residuals + + pred = fit_function(true,a,b) #prediction + pred_err = a_err*true + b_err #a,b uncorrelated + + resid = (ch_charge - pred)/ch_charge + + numerator_err = err_sum(ch_err,pred_err) + + resid_err = err_ratio((ch_charge-pred),ch_charge,numerator_err,ch_std) + + axs[1].errorbar(true,resid*100, yerr=abs(resid_err)*100, ls='', marker = 'o',color = plot_parameters[name]["color"]) #percentage + axs[1].axhline(8, color='k', alpha=0.4, ls='--', lw=2) + axs[1].axhline(-8, color='k', alpha=0.4, ls='--', lw=2) + + + #hg/lg ratio + ratio = charge[:,0]/charge[:,1] + ratio_std = err_ratio(charge[:,0],charge[:,1],std[:,0],std[:,1]) + + yx = zip(true_charge,ratio,ratio_std) #sort by true charge + ratio_sorted = np.array(sorted(yx)) + true = ratio_sorted[:,0] + ratio = ratio_sorted[:,1] + ratio_std = ratio_sorted[:,2] + + model = model = Model(fit_function) + params = model.make_params(a=100, b=0) + ratio_fit = model.fit(ratio[10:-4], params, weights=1/ratio_std[10:-4], x=true[10:-4]) + + + axs[2].set_ylabel("hg/lg") + axs[2].set_xlabel("charge(p.e.)") + axs[2].errorbar(true,ratio,yerr=ratio_std,ls='',color='C1',marker = 'o') + axs[2].plot(true[10:-4],fit_function(true[10:-4], ratio_fit.params['a'].value, ratio_fit.params['b'].value), ls='-',color='C1', alpha=0.9) + axs[2].text(8, 1, 'y = ax + b\n' + + r"$a_{\mathrm{HG/LG}}= (%2.1f\pm%2.1f)$ p.e.$^{-1}$ " %(round(ratio_fit.params['a'].value,1), + round(ratio_fit.params['a'].stderr,1)) + + "\n" r"$b_{\mathrm{HG/LG}}=%2.1f\pm%2.1f$" %(round(ratio_fit.params['b'].value,1), + round(ratio_fit.params['b'].stderr,1)) + + "\n" r"$\chi_{\mathrm{HG/LG}}^2/{\mathrm{dof}} =%2.1f/%d$" %(round(ratio_fit.chisqr,1), + int(ratio_fit.chisqr/ratio_fit.redchi)), + + backgroundcolor='white', bbox=dict(facecolor='white', edgecolor='C1', lw=3, + boxstyle='round,pad=0.5'), + ha='left', va='bottom', fontsize=11, color='k', alpha=0.9) + + + + + plt.savefig(os.path.join(output_dir,"linearity_test.png")) + if temp_output: + with open(os.path.join(args.temp_output, 'plot1.pkl'), 'wb') as f: + pickle.dump(fig, f) + + + + + + #charge resolution + charge_hg = charge[:,0] + std_hg = std[:,0] + std_hg_err = std_err[:,0] + charge_lg = charge[:,1] + std_lg = std[:,1] + std_lg_err = std_err[:,1] + + fig = plt.figure() + charge_plot = np.linspace(3e-2,1e4) + stat_limit = 1/np.sqrt(charge_plot) #statistical limit + + hg_res = (std_hg/charge_hg) + hg_res_err = err_ratio(std_hg,charge_hg,std_hg_err,std_hg) + + lg_res = (std_lg/charge_lg) + lg_res_err = err_ratio(std_lg,charge_lg,std_lg_err,std_lg) + plt.xscale(u'log') + plt.yscale(u'log') + plt.plot(charge_plot,stat_limit,color='gray', ls='-', lw=3,alpha=0.8,label='Statistical limit ') + mask = charge_hg<3e2 + + plt.errorbar(charge_hg[mask], hg_res[mask], + xerr=std_hg[mask]/np.sqrt(npixels),yerr=hg_res_err[mask]/np.sqrt(npixels), + ls = '', marker ='o', label=None,color='C0') + + mask = np.invert(mask) + plt.errorbar(charge_lg[mask]*ratio_fit.params['b'].value,lg_res[mask], + xerr=std_lg[mask]/np.sqrt(npixels),yerr=lg_res_err[mask]/np.sqrt(npixels), ls = '', marker ='o', label=None,color='C0') + + plt.xlabel(r'Charge $\overline{Q}$ [p.e.]') + plt.ylabel(r'Charge resolution $\frac{\sigma_{Q}}{\overline{Q}}$') + + plt.xlim(3e-2,4e3) + plt.legend(frameon=False) + plt.savefig(os.path.join(output_dir,"charge_resolution.png")) + if temp_output: + with open(os.path.join(args.temp_output, 'plot2.pkl'), 'wb') as f: + pickle.dump(fig, f) + plt.close('all') + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/pedestal_test.py b/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/pedestal_test.py new file mode 100644 index 00000000..7d6e2b04 --- /dev/null +++ b/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/pedestal_test.py @@ -0,0 +1,127 @@ +#don't forget to set environment variable NECTARCAMDATA + +import numpy as np +import os + +import matplotlib.pyplot as plt + +from utils import pe2photons,photons2pe +import argparse +import matplotlib.pyplot as plt +import numpy as np +from test_tools_components import PedestalTool +from utils import adc_to_pe, pe2photons +import pandas as pd +from nectarchain.makers.calibration import PedestalNectarCAMCalibrationTool +from ctapipe.containers import EventType +from ctapipe_io_nectarcam.containers import NectarCAMDataContainer +import sys +import pickle + + +def get_args(): + """ + Parses command-line arguments for the linearity test script. + + Returns: + argparse.ArgumentParser: The parsed command-line arguments. + """ + + parser = argparse.ArgumentParser(description='Pedestal substraction test B-TEL-1370.\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 give a list of runs (run numbers with spaces inbetween) and an output directory to save the final plot.\n' + +'If the data is not in NECTARCAMDATA, the files will be downloaded through DIRAC.\n For the purposes of testing this script, default data is from the runs used for this test in the TRR document.\n' + +'You can optionally specify the number of events to be processed (default 1200).\n') + parser.add_argument('-r','--runlist', type=int, nargs='+', help='List of runs (numbers separated by space)', required=False, default = [3647]) + parser.add_argument('-e','--evts', type = int, help='Number of events to process from each run. Default is 1200. 4000 or more gives best results but takes some time', required=False, default=10) + parser.add_argument('-o','--output', type=str, help='Output directory. If none, plot will be saved in current directory', required=False, default='./') + parser.add_argument("--temp_output", help="Temporary output directory for GUI", default=None) + + return parser + + + +def main(): + """ + The main function that runs the pedestal subtraction test. It parses command-line arguments, processes the specified runs, and generates two plots: + + 1. A 2D heatmap of the pedestal RMS for all events and pixels. + 2. A line plot of the mean pedestal RMS for each pixel, with the CTA requirement range highlighted. + + The function also saves the generated plots to the specified output directory, and optionally saves the first plot to a temporary output file. + """ + + parser = get_args() + args = parser.parse_args() + + runlist = args.runlist + nevents = args.evts + + # output_dir = args.output + # temp_output_file = args.temp_output + output_dir = os.path.abspath(args.output) + temp_output = os.path.abspath(args.temp_output) if args.temp_output else None + + print(f"Output directory: {output_dir}") # Debug print + print(f"Temporary output dir: {temp_output}") # Debug print + + + sys.argv = sys.argv[:1] + output = [] + + for run in runlist: + + print("PROCESSING RUN {}".format(run)) + tool = PedestalTool( + progress_bar=True, run_number=run, max_events=nevents, events_per_slice = 999, log_level=20, peak_height=10, window_width=16, overwrite=True + ) + tool.initialize() + print("OUTPUT_PATH", tool.output_path) + tool.setup() + tool.start() + output.append(tool.finish()) + + + rms_ped = pe2photons(np.array(output[0])/adc_to_pe) #in photons + plt.figure() + plt.title("Pedestal rms for all events and pixels") + plt.pcolormesh(rms_ped.T,clim=(0.8,1.5)) + plt.colorbar() + plt.savefig(os.path.join(output_dir,"pedestal_rms_2d_graph.png")) + + + + mean_rms_per_pix = np.mean(rms_ped,axis=0) + mean_value = np.mean(mean_rms_per_pix) + + fig, ax = plt.subplots() + ax.set_title("Mean pedestal rms for each pixel") + ax.set_xlabel("Pixel") + ax.set_ylabel("p.e.") + ax.plot(mean_rms_per_pix, marker='x', linestyle='') + ax.axhline(1.2, color='black', alpha=0.8) + ax.axhline(mean_value, color='red', linestyle='--', alpha=0.5) + + # Fill the region between 0.8 * mean_value and 1.2 * mean_value + ax.fill_between(range(len(mean_rms_per_pix)), 0.8 * mean_value, 1.2 * mean_value, color='red', alpha=0.3) + + right_x_position = len(mean_rms_per_pix) * 0.95 # Slightly left of the right edge of the plot + ax.arrow(right_x_position, 1.2 * mean_value, 0, -0.4 * mean_value, color='grey', width=5, + head_width=50, head_length=0.05, length_includes_head=True, zorder=3) + ax.arrow(right_x_position, 0.8 * mean_value, 0, 0.4 * mean_value, color='grey', width=5, + head_width=50, head_length=0.05, length_includes_head=True, zorder=3) + ax.text(right_x_position + 100, mean_value, "±20%", color='grey', fontsize=12, verticalalignment='center') + ax.text(10, 1.2, "CTA requirement", color='black', fontsize=12) + ax.arrow(300, 1.2, 0, -0.1, color='black', width=5, head_width=50, head_length=0.05, length_includes_head=True, zorder=5) + ax.arrow(1500, 1.2, 0, -0.1, color='black', width=5, head_width=50, head_length=0.05, length_includes_head=True, zorder=5) + + plt.savefig(os.path.join(output_dir,"mean_pedestal_rms.png")) + if temp_output: + with open(os.path.join(args.temp_output, 'plot1.pkl'), 'wb') as f: + pickle.dump(fig, f) + + plt.close('all') + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/pix_couple_tim_uncertainty_test.py b/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/pix_couple_tim_uncertainty_test.py new file mode 100644 index 00000000..532ad9d6 --- /dev/null +++ b/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/pix_couple_tim_uncertainty_test.py @@ -0,0 +1,131 @@ +import matplotlib.pyplot as plt +import numpy as np +from test_tools_components import ToMPairsTool +import os +import sys +import pickle + +import argparse + + +def get_args(): + """ + Parses command-line arguments for the pix_couple_tim_uncertainty_test.py script. + + Returns: + argparse.ArgumentParser: The parsed command-line arguments. + """ + parser = argparse.ArgumentParser(description='Time resolution (timing uncertainty between couples of pixels) test B-TEL-1030.\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 give a list of runs (run numbers with spaces inbetween) and an output directory to save the final plot.\n' + +'If the data is not in NECTARCAMDATA, the files will be downloaded through DIRAC.\n For the purposes of testing this script, default data is from the runs used for this test in the TRR document.\n' + +'You can optionally specify the number of events to be processed (default 1000). Takes a lot of time.\n') + parser.add_argument('-r','--runlist', type=int, nargs='+', help='List of runs (numbers separated by space). You can put just one run, default 3292', required=False, default = [3292]) + parser.add_argument('-e','--evts', type = int, help='Number of events to process from each run. Default is 100. 1000 or more gives best results but takes some time', required=False, default=100) + parser.add_argument('-t','--pmt_transit_time', type=str, help='.csv file with pmt transit time corrections', required=False, default="/Users/dm277349/nectarchain_data/transit_time/hv_pmt_tom_correction_laser_measurement_per_pixel_fit_sqrt_hv_newmethod.csv") + parser.add_argument('-o','--output', type=str, help='Output directory. If none, plot will be saved in current directory', required=False, default='./') + parser.add_argument("--temp_output", help="Temporary output directory for GUI", default=None) + + return parser + + + +def main(): + """ + Generates a plot of the RMS of the time-of-maximum (TOM) difference for pairs of pixels, with a visualization of the CTA requirement. + + The script processes a list of runs, calculates the TOM difference with and without transit time corrections, and plots the distribution of the RMS of the corrected TOM differences. The CTA requirement of 2 ns RMS is visualized on the plot. + + The script takes several command-line arguments, including the list of runs to process, the number of events to process per run, the path to a CSV file with PMT transit time corrections, and the output directory for the plot. + + If a temporary output directory is specified, the plot is also saved to a pickle file in that directory for the gui to use. + """ + + + parser = get_args() + args = parser.parse_args() + + + + tt_path = "/Users/dm277349/nectarchain_data/transit_time/hv_pmt_tom_correction_laser_measurement_per_pixel_fit_sqrt_hv_newmethod.csv" + + runlist = args.runlist + nevents = args.evts + tt_path = args.pmt_transit_time + output_dir = os.path.abspath(args.output) + temp_output = os.path.abspath(args.temp_output) if args.temp_output else None + + print(f"Output directory: {output_dir}") # Debug print + print(f"Temporary output file: {temp_output}") # Debug print + + + sys.argv = sys.argv[:1] + tom = [] + + pixel_ids = [] #pixel ids for run + + tom_corrected = [] + + dt_no_correction = [] + dt_corrected = [] + pixel_pairs = [] + + + for run in runlist: + + print("PROCESSING RUN {}".format(run)) + tool = ToMPairsTool( + progress_bar=True, run_number=run, events_per_slice = 501 ,max_events=nevents, log_level=20, peak_height=10, window_width=16, overwrite=True + ) + tool.initialize() + tool.setup() + tool.start() + output = tool.finish(tt_path) + tom.append(output[0]) + tom_corrected.append(output[1]) + pixel_ids.append(output[2]) + dt_no_correction.append(output[3]) + dt_corrected.append(output[4]) + pixel_pairs.append(output[5]) + + + dt_no_correction = np.array(dt_no_correction[0]) + dt_corrected = np.array(output[4]) + dt_corrected[abs(dt_corrected)>7]=np.nan + dt_corrected[abs(dt_no_correction)>7]=np.nan + std_corrected = np.nanstd(dt_corrected, axis = 1) + + fig, ax = plt.subplots(figsize=(10,10/1.61)) + plt.hist(std_corrected,range=(0,5),density=True, + histtype='step', lw=3, bins=200) + + plt.axvline(2, color='C4', alpha=0.8) + ax.text(2.1, 0.5, 'CTA requirement', color='C4', fontsize=20, + rotation=-90) + + + ax.annotate("", xy=(1.6, 0.25), xytext=(2, 0.25), color='C4', alpha=0.5, + transform=ax.transAxes, + arrowprops=dict(color='C4', alpha=0.7, lw=3, arrowstyle='->')) + + ax.annotate("", xy=(1.6, 0.75), xytext=(2, 0.75), color='C4', alpha=0.5, + transform=ax.transAxes, + arrowprops=dict(color='C4', alpha=0.7, lw=3, arrowstyle='->')) + + plt.xlabel('RMS of $\Delta t_{\mathrm{TOM}}$ for pairs of pixels [ns]') + plt.ylabel('Normalized entries') + + plt.gcf() + + plt.savefig(os.path.join(output_dir,"pix_couple_tim_uncertainty.png")) + + if temp_output: + with open(os.path.join(args.temp_output, 'plot1.pkl'), 'wb') as f: + pickle.dump(fig, f) + + plt.close('all') + + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/pix_tim_uncertainty_test.py b/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/pix_tim_uncertainty_test.py new file mode 100644 index 00000000..b5e32b21 --- /dev/null +++ b/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/pix_tim_uncertainty_test.py @@ -0,0 +1,305 @@ +#don't forget to set environment variable NECTARCAMDATA + +import numpy as np +import os + +import matplotlib.pyplot as plt + +from utils import pe2photons,photons2pe +from test_tools_components import TimingResolutionTestTool +import argparse +import sys +import pickle + + +def get_args(): + """ + Parses command-line arguments for the pixel timing uncertainty test script. + + Returns: + argparse.ArgumentParser: The parsed command-line arguments. + """ + parser = argparse.ArgumentParser(description='Systematic pixel timing uncertainty test B-TEL-1380.\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 give a list of runs (run numbers with spaces inbetween) and an output directory to save the final plot.\n' + +'If the data is not in NECTARCAMDATA, the files will be downloaded through DIRAC.\n For the purposes of testing this script, default data is from the runs used for this test in the TRR document.\n' + +'You can optionally specify the number of events to be processed (default 1200) and the number of pixels used (default 70).\n') + parser.add_argument('-r','--runlist', type=int, nargs='+', help='List of runs (numbers separated by space)', required=False, default = [i for i in range (3446,3457)]) + parser.add_argument('-e','--evts', type = int, help='Number of events to process from each run. Default is 100', required=False, default=100) + #parser.add_argument('-p','--pixels', type = int, help='Number of pixels used. Default is 70', required=False, default=70) + parser.add_argument('-o','--output', type=str, help='Output directory. If none, plot will be saved in current directory', required=False, default='./') + parser.add_argument("--temp_output", help="Temporary output directory for GUI", default=None) + return parser + + +def main(): + """ + Processes the pixel timing uncertainty test data and generates a plot. + + The function processes the data from the specified list of runs, calculates the weighted mean RMS and RMS error, and generates a plot of the results. The plot is saved to the specified output directory. + + If a temporary output directory is provided, the plot is also saved to a pickle file in that directory for the gui to use. + """ + + parser = get_args() + args = parser.parse_args() + + + runlist = args.runlist + nevents = args.evts + #npixels = args.pixels + output_dir = os.path.abspath(args.output) + temp_output = os.path.abspath(args.temp_output) if args.temp_output else None + + print(f"Output directory: {output_dir}") # Debug print + print(f"Temporary output file: {temp_output}") # Debug print + + + sys.argv = sys.argv[:1] + + # rms_mu = [] + # rms_mu_err = [] + rms_no_fit = [] + rms_no_fit_err = [] + mean_charge_pe = [] + + + for run in runlist: + print("PROCESSING RUN {}".format(run)) + tool = TimingResolutionTestTool( + progress_bar=True, run_number=run, max_events=nevents, events_per_slice = 999, log_level=20, window_width=16, overwrite=True + ) + tool.initialize() + tool.setup() + tool.start() + output = tool.finish() + # rms_mu.append(output[0]) + # rms_mu_err.append(output[1]) + rms_no_fit.append(output[0]) + rms_no_fit_err.append(output[1]) + mean_charge_pe.append(output[2]) + + + print(rms_no_fit) + photons_spline = np.array(mean_charge_pe)*100/25 + rms_no_fit_err=np.array(rms_no_fit_err) + print(rms_no_fit_err) + rms_no_fit_err[rms_no_fit_err==0]=1e-5 #almost zero + #rms_no_fit_err[rms_no_fit_err==np.nan]=1e-5 + print(rms_no_fit_err) + + # mean_rms_mu = np.mean(rms_mu,axis=1) + # mean_rms_no_fit = np.mean(rms_no_fit,axis=1) + + # weights_mu_pix = 1/(np.array(rms_mu_err)+1e-5)**2 + weights_no_fit_pix = 1/(rms_no_fit_err)**2 + weights_no_fit_pix[weights_no_fit_pix > 1e5] = 1e5 + print(weights_no_fit_pix) + + # rms_mu_weighted=[] + # rms_mu_weighted_err=[] + rms_no_fit_weighted=[] + rms_no_fit_weighted_err = [] + + for run in range(len(runlist)): + # rms_mu_weighted.append(np.sum(rms_mu[run]*weights_mu_pix[run])/np.sum(weights_mu_pix[run])) + # rms_mu_weighted_err.append(np.sqrt(1/np.sum(weights_mu_pix[run]))) + rms_no_fit_weighted.append(np.nansum(rms_no_fit[run]*weights_no_fit_pix[run])/np.nansum(weights_no_fit_pix[run])) + rms_no_fit_weighted_err.append(np.sqrt(1/np.nansum(weights_no_fit_pix[run]))) + + print(rms_no_fit_weighted) + print(rms_no_fit_weighted_err) + + + #FIGURE + fig, ax = plt.subplots(figsize=(10,7), constrained_layout=True) + + + plt.errorbar(x=photons_spline[:], + y=np.sqrt(np.array(rms_no_fit_weighted[:])**2), + yerr = rms_no_fit_weighted_err, + ls='', marker='o', + label = r'$\mathtt{scipy.signal.find\_peaks}$') + # plt.errorbar(x=photons_spline[:], + # y=np.sqrt(np.array(rms_mu_weighted[:])**2), + # yerr=rms_mu_weighted_err, + # ls='', marker='o', + # label='Gaussian fit') + + + plt.axhline(1, ls='--', color='C4', alpha=0.6) + plt.axhline(1/np.sqrt(12), ls='--', color='gray', alpha=0.7, label='Quantification rms noise') + + + plt.axvspan(20, 1000, alpha=0.1, color='C4') + + ax.text(51.5, 1.04, 'CTA requirement', color='C4', fontsize=20, + horizontalalignment='left', + verticalalignment='center') + ax.annotate("", xy=(40, 0.9), xytext=(40, 0.995), color='C4', alpha=0.5, + arrowprops=dict(color='C4', alpha=0.7, lw=3, arrowstyle='->')) + + ax.annotate("", xy=(200, 0.9), xytext=(200, 0.995), color='C4', alpha=0.5, + arrowprops=dict(color='C4', alpha=0.7, lw=3, arrowstyle='->')) + + plt.legend(frameon=True, prop={'size':18}, loc='upper right', handlelength=1.2) + plt.xlabel('Illumination charge [photons]' ) + plt.ylabel('Mean rms per pixel [ns]' ) + plt.xscale(u'log') + plt.ylim((0,2.7)) + secax = ax.secondary_xaxis('top', functions=(pe2photons, photons2pe)) + secax.set_xlabel('Illumination charge [p.e.]', labelpad=7) + plt.savefig(os.path.join(output_dir,"pix_tim_uncertainty.png")) + + if temp_output: + with open(os.path.join(args.temp_output, 'plot1.pkl'), 'wb') as f: + pickle.dump(fig, f) + + plt.close('all') + + +if __name__ == "__main__": + main() + + + + + + + + + + +# #don't forget to set environment variable NECTARCAMDATA + +# import numpy as np +# import os +# import sys + +# import matplotlib.pyplot as plt + +# from utils import pe2photons,photons2pe +# from test_tools_components import TimingResolutionTestTool +# import argparse + +# parser = argparse.ArgumentParser(description='Systematic pixel timing uncertainty test B-TEL-1380.\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 give a list of runs (run numbers with spaces inbetween) and an output directory to save the final plot.\n' +# +'If the data is not in NECTARCAMDATA, the files will be downloaded through DIRAC.\n For the purposes of testing this script, default data is from the runs used for this test in the TRR document.\n' +# +'You can optionally specify the number of events to be processed (default 1200) and the number of pixels used (default 70).\n') +# parser.add_argument('-r','--runlist', type = int, nargs='+', help='List of runs (numbers separated by space)', required=False) +# parser.add_argument('-e','--evts', type = int, help='Number of events to process from each run. Default is 1200. 4000 or more gives best results but takes some time', required=False, default=1200) +# #parser.add_argument('-p','--pixels', type = int, help='Number of pixels used. Default is 70', required=False, default=70) +# parser.add_argument('-o','--output', type=str, help='Output directory. If none, plot will be saved in current directory', required=False, default='./') + +# args = parser.parse_args() + +# runlist = args.runlist +# nevents = args.evts +# #npixels = args.pixels + +# output_dir = args.output + +# sys.argv = sys.argv[:1] + +# # rms_mu = [] +# # rms_mu_err = [] +# rms_no_fit = [] +# rms_no_fit_err = [] +# mean_charge_pe = [] + + +# for run in runlist: +# print("PROCESSING RUN {}".format(run)) +# tool = TimingResolutionTestTool( +# progress_bar=True, run_number=run, max_events=nevents, log_level=20, window_width=16, overwrite=True +# ) +# tool.initialize() +# tool.setup() +# tool.start() +# output = tool.finish() +# print(output) +# print(len(output[0]),len(output[1])) +# # rms_mu.append(output[0]) +# # rms_mu_err.append(output[1]) +# rms_no_fit.append(output[0]) +# rms_no_fit_err.append(output[1]) +# mean_charge_pe.append(output[2]) + +# rms_no_fit_err=np.array(rms_no_fit_err) +# photons_spline = np.array(mean_charge_pe)*100/25 +# print(photons_spline) +# # mean_rms_mu = np.mean(rms_mu,axis=1) +# # mean_rms_no_fit = np.mean(rms_no_fit,axis=1) +# rms_no_fit_err[rms_no_fit_err==0]=1e-5 #almost zero +# rms_no_fit_err[rms_no_fit_err==np.nan]=1e-5 + +# # print(np.isnan(rms_no_fit_err).any()) +# # print(np.where(rms_no_fit_err==0)) +# # weights_mu_pix = 1/(np.array(rms_mu_err)+1e-5)**2 +# weights_no_fit_pix = (1/(np.array(rms_no_fit_err)))**2 +# # print(np.isnan(weights_no_fit_pix)) +# # print(weights_no_fit_pix.dtype) +# # print(np.isnan(weights_no_fit_pix).any()) +# # print(np.max(weights_no_fit_pix)) +# # print(np.isinf(weights_no_fit_pix).any()) +# # print(np.where(np.isinf(weights_no_fit_pix))) +# # for i in range(len(weights_no_fit_pix[0])): +# # print(i) +# # print(weights_no_fit_pix[0][i]) +# # print(np.nansum(weights_no_fit_pix[0])) +# # print(np.sqrt(weights_no_fit_pix[0])) + +# # rms_mu_weighted=[] +# # rms_mu_weighted_err=[] +# rms_no_fit_weighted=[] +# rms_no_fit_weighted_err = [] + +# for run in range(len(runlist)): +# # rms_mu_weighted.append(np.sum(rms_mu[run]*weights_mu_pix[run])/np.sum(weights_mu_pix[run])) +# # rms_mu_weighted_err.append(np.sqrt(1/np.sum(weights_mu_pix[run]))) +# rms_no_fit_weighted.append(np.nansum(rms_no_fit[run]*weights_no_fit_pix[run])/np.nansum(weights_no_fit_pix[run])) +# rms_no_fit_weighted_err.append(np.sqrt(1/np.nansum(weights_no_fit_pix[run]))) + + + +# print("RESULT", rms_no_fit_weighted) +# #FIGURE +# fig, ax = plt.subplots(figsize=(10,7), constrained_layout=True) + + +# plt.errorbar(x=photons_spline[:], +# y=np.sqrt(np.array(rms_no_fit_weighted[:])**2), +# yerr = rms_no_fit_weighted_err, +# ls='', marker='o', +# label = r'$\mathtt{scipy.signal.find\_peaks}$') +# # plt.errorbar(x=photons_spline[:], +# # y=np.sqrt(np.array(rms_mu_weighted[:])**2), +# # yerr=rms_mu_weighted_err, +# # ls='', marker='o', +# # label='Gaussian fit') + + +# plt.axhline(1, ls='--', color='C4', alpha=0.6) +# plt.axhline(1/np.sqrt(12), ls='--', color='gray', alpha=0.7, label='Quantification rms noise') + + +# plt.axvspan(20, 1000, alpha=0.1, color='C4') + +# ax.text(51.5, 1.04, 'CTA requirement', color='C4', fontsize=20, +# horizontalalignment='left', +# verticalalignment='center') +# ax.annotate("", xy=(40, 0.9), xytext=(40, 0.995), color='C4', alpha=0.5, +# arrowprops=dict(color='C4', alpha=0.7, lw=3, arrowstyle='->')) + +# ax.annotate("", xy=(200, 0.9), xytext=(200, 0.995), color='C4', alpha=0.5, +# arrowprops=dict(color='C4', alpha=0.7, lw=3, arrowstyle='->')) + +# plt.legend(frameon=True, prop={'size':18}, loc='upper right', handlelength=1.2) +# plt.xlabel('Illumination charge [photons]' ) +# plt.ylabel('Mean rms per pixel [ns]' ) +# plt.xscale(u'log') +# #plt.ylim((0,2.7)) +# secax = ax.secondary_xaxis('top', functions=(pe2photons, photons2pe)) +# secax.set_xlabel('Illumination charge [p.e.]', labelpad=7) +# plt.savefig(os.path.join(output_dir,"pix_tim_uncertainty.png")) diff --git a/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/test_tools_components.py b/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/test_tools_components.py new file mode 100644 index 00000000..d827779d --- /dev/null +++ b/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/test_tools_components.py @@ -0,0 +1,1291 @@ +from astropy import units as u +from ctapipe.containers import EventType, Field +from ctapipe.core.traits import ComponentNameList, Integer +from ctapipe_io_nectarcam import constants +from ctapipe_io_nectarcam.containers import NectarCAMDataContainer +from itertools import combinations +from lmfit.models import Model +from nectarchain.data.container import ( + NectarCAMContainer, +) +from nectarchain.makers import EventsLoopNectarCAMCalibrationTool +from nectarchain.makers.component import NectarCAMComponent +from nectarchain.utils import ComponentUtils +import h5py +import hdf5plugin +import matplotlib.pyplot as plt +import numpy as np +import os +import pandas as pd +import pathlib +import random +from scipy.interpolate import InterpolatedUnivariateSpline +from scipy.signal import find_peaks +from utils import ( + adc_to_pe, + argmedian, + closest_value, +) + + +#overriding so we can have maxevents in the path +def _init_output_path(self): + """ + Initializes the output path for the NectarCAMCalibrationTool. + + If `max_events` is `None`, the output file name will be in the format `{self.name}_run{self.run_number}.h5`. Otherwise, the file name will be in the format `{self.name}_run{self.run_number}_maxevents{self.max_events}.h5`. + + The output path is constructed by joining the `NECTARCAMDATA` environment variable (or `/tmp` if not set) with the `tests` subdirectory and the generated file name. + """ + + if self.max_events is None: + filename = f"{self.name}_run{self.run_number}.h5" + else: + filename = f"{self.name}_run{self.run_number}_maxevents{self.max_events}.h5" + self.output_path = pathlib.Path( + f"{os.environ.get('NECTARCAMDATA','/tmp')}/tests/{filename}" + ) + +EventsLoopNectarCAMCalibrationTool._init_output_path = _init_output_path + + + + +class ChargeContainer(NectarCAMContainer): + + """ + This class contains fields that store various properties and data related to NectarCAM events, including: + + - `run_number`: The run number associated with the waveforms. + - `npixels`: The number of effective pixels. + - `pixels_id`: An array of pixel IDs. + - `ucts_timestamp`: An array of UCTS timestamps for the events. + - `event_type`: An array of trigger event types. + - `event_id`: An array of event IDs. + - `charge_hg`: A 2D array of high gain charge values. + - `charge_lg`: A 2D array of low gain charge values. + """ + run_number = Field( + type=np.uint16, + description="run number associated to the waveforms", + ) + npixels = Field( + type=np.uint16, + description="number of effective pixels", + ) + pixels_id = Field(type=np.ndarray, dtype=np.uint16, ndim=1, description="pixel ids") + ucts_timestamp = Field( + type=np.ndarray, dtype=np.uint64, ndim=1, description="events ucts timestamp" + ) + event_type = Field( + type=np.ndarray, dtype=np.uint8, ndim=1, description="trigger event type" + ) + event_id = Field(type=np.ndarray, dtype=np.uint32, ndim=1, description="event ids") + + + charge_hg = Field( + type=np.ndarray, dtype=np.float64, ndim=2, description="The high gain charge" + ) + charge_lg = Field( + type=np.ndarray, dtype=np.float64, ndim=2, description="The low gain charge" + ) + + +class ChargeComp(NectarCAMComponent): + + """ + This class `ChargeComp` is a NectarCAMComponent that processes NectarCAM event data. It extracts the charge information from the waveforms of each event, handling cases of saturated or noisy events. The class has the following configurable parameters: + + - `window_shift`: The time in ns before the peak to extract the charge. + - `window_width`: The duration of the charge extraction window in ns. + + The `__init__` method initializes important members of the component, such as timestamps, event type, event ids, pedestal and charge for both gain channels. + The `__call__` method is the main processing logic, which is called for each event. It extracts the charge information for both high gain and low gain channels, handling various cases such as saturated events and events with no signal. + The `finish` method collects all the processed data and returns a `ChargeContainer` object containing the run number, number of pixels, pixel IDs, UCTS timestamps, event types, event IDs, and the high and low gain charge values. + """ + window_shift = Integer( + default_value=6, + help="the time in ns before the peak to extract charge", + ).tag(config=True) + + window_width = Integer( + default_value=12, + help="the duration of the extraction window in ns", + ).tag(config=True) + + + + + def __init__(self, subarray, config=None, parent=None, *args, **kwargs): + super().__init__( + subarray=subarray, config=config, parent=parent, *args, **kwargs + ) + ## If you want you can add here members of MyComp, they will contain interesting quantity during the event loop process + + self.__ucts_timestamp = [] + self.__event_type = [] + self.__event_id = [] + + self.__pedestal_hg = [] + self.__pedestal_lg = [] + + self.__charge_hg = [] + self.__charge_lg = [] + + + + ##This method need to be defined ! + def __call__(self, event: NectarCAMDataContainer, *args, **kwargs): + self.__event_id.append(np.uint32(event.index.event_id)) + + #print(event.index.event_id) + + self.__event_type.append(event.trigger.event_type.value) + self.__ucts_timestamp.append(event.nectarcam.tel[0].evt.ucts_timestamp) + + wfs = [] + wfs.append(event.r0.tel[0].waveform[constants.HIGH_GAIN][self.pixels_id]) + wfs.append(event.r0.tel[0].waveform[constants.LOW_GAIN][self.pixels_id]) + + + + + #####THE JOB IS HERE###### + for i, (pedestal,charge) in enumerate(zip([self.__pedestal_hg, self.__pedestal_lg],[self.__charge_hg,self.__charge_lg])): + wf = np.array(wfs[i],dtype=np.float16) + + index_peak = np.argmax(wf,axis=1) + index_peak[index_peak<20]=20 + index_peak[index_peak>40]=40 + + signal_start = index_peak - self.window_shift + signal_stop = index_peak + self.window_width - self.window_shift + + + + + if event.trigger.event_type==EventType.FLATFIELD: + integral = np.zeros(len(self.pixels_id)) + + for pix in range(len(self.pixels_id)): + #search for saturated events or events with no signal + ped = np.round(np.mean(wf[pix, 0:signal_start[pix]])) + + + peaks_sat = find_peaks(wf[pix,20:45],height=1000,plateau_size=self.window_width) + if len(peaks_sat[0])==1: + #saturated + + #print("saturated") + + signal_start[pix]=argmedian(wf[pix,20:45])+20-self.window_shift + signal_stop[pix]=signal_start[pix]+self.window_width + integral[pix] = np.sum(wf[pix,signal_start[pix]:signal_stop[pix]])-ped*(signal_stop[pix]-signal_start[pix]) + + else: + peaks_signal = find_peaks(wf[pix],prominence=10) + if len(peaks_signal[0])>=12: + #print("noisy event") + integral[pix]=0 + + elif len(peaks_signal[0])<1: + #flat + integral[pix]=0 + + else: + # x = np.linspace(0,signal_stop[pix]-signal_start[pix],signal_stop[pix]-signal_start[pix]) + # spl = UnivariateSpline(x,y) + # integral[pix] = spl.integral(0,signal_stop[pix]-signal_start[pix]) + + + integral[pix] = np.sum(wf[pix,signal_start[pix]:signal_stop[pix]])-ped*(signal_stop[pix]-signal_start[pix]) + + + chg = integral + + charge.append(chg) + + + + + + + + + + + + + + ##This method need to be defined ! + def finish(self): + + output = ChargeContainer( + run_number=ChargeContainer.fields["run_number"].type(self._run_number), + npixels=ChargeContainer.fields["npixels"].type(self._npixels), + pixels_id=ChargeContainer.fields["pixels_id"].dtype.type(self._pixels_id), + ucts_timestamp=ChargeContainer.fields["ucts_timestamp"].dtype.type( + self.__ucts_timestamp + ), + event_type=ChargeContainer.fields["event_type"].dtype.type(self.__event_type), + event_id=ChargeContainer.fields["event_id"].dtype.type(self.__event_id), + + charge_hg=ChargeContainer.fields["charge_hg"].dtype.type( + np.array(self.__charge_hg) + ), + charge_lg=ChargeContainer.fields["charge_lg"].dtype.type( + np.array(self.__charge_lg) + ), + ) + + + return output + + + +class LinearityTestTool(EventsLoopNectarCAMCalibrationTool): + """ + This class, `LinearityTestTool`, is a subclass of `EventsLoopNectarCAMCalibrationTool`. It is responsible for performing a linearity test on NectarCAM data. The class has a `componentsList` attribute that specifies the list of NectarCAM components to be applied. + + The `finish` method is the main functionality of this class. It reads the charge data from the output file, calculates the mean charge, standard deviation, and standard error for both the high gain and low gain channels, and returns these values. This information can be used to assess the linearity of the NectarCAM system. + """ + + name = "LinearityTestTool" + + componentsList = ComponentNameList( + NectarCAMComponent, + default_value=["ChargeComp"], + help="List of Component names to be apply, the order will be respected", + ).tag(config=True) + + + def finish(self, *args, **kwargs): + + super().finish(return_output_component=True, *args, **kwargs) + + + mean_charge = [0,0] #per channel + std_charge = [0,0] + std_err = [0,0] + + charge_hg = [] + charge_lg = [] + + + output_file = h5py.File(self.output_path) + + + for thing in output_file: + group = output_file[thing] + dataset = group['ChargeContainer'] + data = dataset[:] + #print("data",data) + for tup in data: + try: + npixels = tup[1] + charge_hg.extend(tup[6]) + charge_lg.extend(tup[7]) + + except: + break + + output_file.close() + + + charge_pe_hg=np.array(charge_hg)/adc_to_pe + charge_pe_lg=np.array(charge_lg)/adc_to_pe + + for channel,charge in enumerate([charge_pe_hg,charge_pe_lg]): + + + pix_mean_charge = np.mean(charge, axis=0) #in pe + + pix_std_charge = np.std(charge, axis=0) + + #average of all pixels + mean_charge[channel] = np.mean(pix_mean_charge) + + std_charge[channel] = np.mean(pix_std_charge) + #for the charge resolution + std_err[channel] = np.std(pix_std_charge) + + return mean_charge, std_charge, std_err, npixels + + + + + +class ToMContainer(NectarCAMContainer): + """ + Attributes: + run_number (np.uint16): The run number associated with the waveforms. + npixels (np.uint16): The number of effective pixels. + pixels_id (np.ndarray[np.uint16]): The pixel IDs. + ucts_timestamp (np.ndarray[np.uint64]): The UCTS timestamps of the events. + event_type (np.ndarray[np.uint8]): The trigger event types. + event_id (np.ndarray[np.uint32]): The event IDs. + charge_hg (np.ndarray[np.float64]): The mean high gain charge per event. + tom_no_fit (np.ndarray[np.float64]): The time of maximum from the data (no fitting). + good_evts (np.ndarray[np.uint32]): The IDs of the good (non-cosmic ray) events. + """ + + run_number = Field( + type=np.uint16, + description="run number associated to the waveforms", + ) + npixels = Field( + type=np.uint16, + description="number of effective pixels", + ) + pixels_id = Field(type=np.ndarray, dtype=np.uint16, ndim=1, description="pixel ids") + ucts_timestamp = Field( + type=np.ndarray, dtype=np.uint64, ndim=1, description="events ucts timestamp" + ) + event_type = Field( + type=np.ndarray, dtype=np.uint8, ndim=1, description="trigger event type" + ) + event_id = Field(type=np.ndarray, dtype=np.uint32, ndim=1, description="event ids") + + + charge_hg = Field( + type=np.ndarray, dtype=np.float64, ndim=2, description="The mean high gain charge per event" + ) + + # tom_mu = Field( + # type=np.ndarray, dtype=np.float64, ndim=2, description="Time of maximum of signal fitted with gaussian" + # ) + + # tom_sigma = Field( + # type=np.ndarray, dtype=np.float64, ndim=2, description="Time of fitted maximum sigma" + # ) + tom_no_fit = Field( + type=np.ndarray, dtype=np.float64, ndim=2, description="Time of maximum from data (no fitting)" + ) + good_evts = Field(type=np.ndarray, dtype=np.uint32, ndim=1, description="good (non cosmic ray) event ids") + + + +class ToMComp(NectarCAMComponent): + """ + This class, `ToMComp`, is a component of the NectarCAM system that is responsible for processing waveform data. It has several configurable parameters, including the width and shift before the peak of the time window for charge extraction, the peak height threshold. + + The `__init__` method initializes some important component members, such as timestamps, event type, event ids, pedestal and charge values for both gain channels. + + The `__call__` method is the main entry point for processing an event. It extracts the waveform data, calculates the pedestal, charge, and time of maximum (ToM) for each pixel, and filters out events that do not meet the peak height threshold. The results are stored in various member variables, which are then returned in the `finish` method. + + The `finish` method collects the processed data from the member variables and returns a `ToMContainer` object, which contains the run number, number of pixels, pixel IDs, UCTS timestamps, event types, event IDs, high-gain charge, ToM without fitting, and IDs of good (non-cosmic ray) events. + """ + window_shift = Integer( + default_value=6, + help="the time in ns before the peak to extract charge", + ).tag(config=True) + + window_width = Integer( + default_value=16, + help="the duration of the extraction window in ns", + ).tag(config=True) + + peak_height = Integer( + default_value=10, + help="height of peak to consider event not to be just pedestal (ADC counts)", + ).tag(config=True) + + def __init__(self, subarray, config=None, parent=None, *args, **kwargs): + super().__init__( + subarray=subarray, config=config, parent=parent, *args, **kwargs + ) + ## If you want you can add here members of MyComp, they will contain interesting quantity during the event loop process + + self.__ucts_timestamp = [] + self.__event_type = [] + self.__event_id = [] + + self.__charge_hg = [] + self.__pedestal_hg = [] + + # self.__tom_mu = [] + + # self.__tom_sigma = [] + + self.__tom_no_fit = [] + + self.__good_evts = [] + + self.__ff_event_ind = -1 + + + + + ##This method need to be defined ! + def __call__(self, event: NectarCAMDataContainer, *args, **kwargs): + self.__event_id.append(np.uint32(event.index.event_id)) + + + + + self.__event_type.append(event.trigger.event_type.value) + self.__ucts_timestamp.append(event.nectarcam.tel[0].evt.ucts_timestamp) + + + wfs = [] + wfs.append(event.r0.tel[0].waveform[constants.HIGH_GAIN][self.pixels_id]) + + self.__ff_event_ind += 1 + + + #####THE JOB IS HERE###### + + for i, (pedestal,charge,tom_no_fit) in enumerate(zip([self.__pedestal_hg],[self.__charge_hg],[self.__tom_no_fit])): + wf = np.array(wfs[i]) #waveform per gain + index_peak = np.argmax(wf,axis=1) #tom per event/pixel + #use it to first find pedestal and then will filter out no signal events + index_peak[index_peak<20]=20 + index_peak[index_peak>40]=40 + signal_start = index_peak - self.window_shift + signal_stop = index_peak + self.window_width - self.window_shift + + ped = np.array([np.mean(wf[pix, 0:signal_start[pix]]) for pix in range(len(self.pixels_id))]) + + pedestal.append(ped) + + + tom_no_fit_evt = np.zeros(len(self.pixels_id)) + chg = np.zeros(len(self.pixels_id)) + #will calculate tom & charge like federica + for pix in range(len(self.pixels_id)): + y = wf[pix] - ped[pix]#np.maximum(wf[pix] - ped[pix],np.zeros(len(wf[pix]))) + + + x = np.linspace(0, len(y), len(y)) + xi = np.linspace(0, len(y), 251) + ius = InterpolatedUnivariateSpline(x, y) + yi = ius(xi) + peaks, _ = find_peaks(yi, height=self.peak_height, ) + #print(peaks) + peaks = peaks[xi[peaks]>20] + peaks = peaks[xi[peaks]<32] + #print(peaks) + + if len(peaks)>0: + # find the max peak + # max_peak = max(yi[peaks]) + max_peak_index = np.argmax(yi[peaks], axis=0) + + # Check if there is not a peak but a plateaux + yi_rounded = np.around(yi[peaks]/max(yi[peaks]),1) + maxima_peak_index = np.argwhere(yi_rounded== np.amax(yi_rounded)) + + + #saturated events + if ((len(maxima_peak_index) > 1) and (min(xi[peaks[maxima_peak_index]]) > signal_start[pix]) and ( + max(xi[peaks[maxima_peak_index]]) < signal_stop[pix])): + # saturated event + + max_peak_index = int(np.median(maxima_peak_index)) + + + # simple sum integration + x_max = closest_value(x, xi[peaks[max_peak_index]]) # find the maximum in the not splined array + charge_sum = y[(list(x).index(x_max) - (self.window_width - 10)):(list(x).index(x_max) + (self.window_width - 6))].sum() + + # gaussian fit + # change_grad_pos_left = 3 + # change_grad_pos_right = 3 + + # mask = np.ma.masked_where(y > (4095 - 400), x) + # x_fit = np.ma.compressed(mask) + # mask = np.ma.masked_where(y > (4095 - 400), y) + # y_fit = np.ma.compressed(mask) + # mean = xi[peaks[max_peak_index]] + # sigma = change_grad_pos_right + change_grad_pos_left + + # # fit + # model = Model(gaus) + # params = model.make_params(a=yi[peaks[max_peak_index]] * 3, mu=mean, sigma=sigma) + # result = model.fit(y_fit, params, x=x_fit) + + # result_sigma = result.params['sigma'].value + # result_mu = result.params['mu'].value + + + max_position_x_prefit = xi[peaks[max_peak_index]] + + + elif ((len(maxima_peak_index) == 1) and (xi[peaks[max_peak_index]] > signal_start[pix]) and ( + xi[peaks[max_peak_index]]< signal_stop[pix])): + + # simple sum integration + x_max = closest_value(x, xi[peaks[max_peak_index]]) # find the maximum in the not splined array + charge_sum = y[(list(x).index(x_max) - (self.window_width - 10)):(list(x).index(x_max) + (self.window_width - 6))].sum() + + + + # gaussian fit + # change_grad_pos_left = 3 + # change_grad_pos_right = 3 + # mean = xi[peaks[max_peak_index]] + # sigma = change_grad_pos_right + change_grad_pos_left # define window for the gaussian fit + + # x_fit = xi[peaks[max_peak_index]-change_grad_pos_left:peaks[max_peak_index]+change_grad_pos_right] + # y_fit = yi[peaks[max_peak_index]-change_grad_pos_left:peaks[max_peak_index]+change_grad_pos_right] + # model = Model(gaus) + # params = model.make_params(a=yi[peaks[max_peak_index]],mu=mean,sigma=sigma) + # result = model.fit(y_fit, params, x=x_fit) + + max_position_x_prefit = xi[peaks[max_peak_index]] + # result_sigma = result.params['sigma'].value + # result_mu = result.params['mu'].value + + + + else: + + # index_x_window_min = list(xi).index(closest_value(xi, signal_start[pix])) + charge_sum = y[signal_start[pix]:signal_start[pix] + self.window_width].sum() + + + + max_position_x_prefit = -1 + # result_sigma = -1 + # result_mu = -1 + + + else: + + # If no maximum is found, the integration is done between 20 and 36 ns. + signal_start[pix]=20 + + # index_x_window_min = list(xi).index(closest_value(xi, signal_start[pix])) + charge_sum = y[signal_start[pix]:signal_start[pix] + self.window_width].sum() + + max_position_x_prefit = -1 + + # result_sigma = -1 + # result_mu = -1 + + + # tom_mu_evt[pix] = result_mu + # tom_sigma_evt[pix] = result_sigma + tom_no_fit_evt[pix] = max_position_x_prefit + chg[pix] = charge_sum + + + + + # tom_mu.append(tom_mu_evt) + # tom_sigma.append(tom_sigma_evt) + tom_no_fit.append(tom_no_fit_evt) + #print("tom for event", tom_no_fit_evt) + charge.append(chg) + + + #is it a good event? + if np.max(chg)<10*np.mean(chg): + #print("is good evt") + self.__good_evts.append(self.__ff_event_ind) + + + + + + + ##This method need to be defined ! + def finish(self): + + output = ToMContainer( + run_number=ToMContainer.fields["run_number"].type(self._run_number), + npixels=ToMContainer.fields["npixels"].type(self._npixels), + pixels_id=ToMContainer.fields["pixels_id"].dtype.type(self.pixels_id), + ucts_timestamp=ToMContainer.fields["ucts_timestamp"].dtype.type( + self.__ucts_timestamp + ), + event_type=ToMContainer.fields["event_type"].dtype.type(self.__event_type), + event_id=ToMContainer.fields["event_id"].dtype.type(self.__event_id), + + charge_hg=ToMContainer.fields["charge_hg"].dtype.type( + self.__charge_hg + ), + + # tom_mu=ToMContainer.fields["tom_mu"].dtype.type( + # self.__tom_mu + # ), + + # tom_sigma=ToMContainer.fields["tom_sigma"].dtype.type( + # self.__tom_sigma + # ), + + tom_no_fit=ToMContainer.fields["tom_no_fit"].dtype.type( + self.__tom_no_fit + ), + + good_evts=ToMContainer.fields["good_evts"].dtype.type( + self.__good_evts + ), + + ) + return output + + + + + +class TimingResolutionTestTool(EventsLoopNectarCAMCalibrationTool): + + """ + This class, `TimingResolutionTestTool`, is a subclass of `EventsLoopNectarCAMCalibrationTool` and is used to perform timing resolution tests on NectarCAM data. It reads the output data from the `ToMContainer` dataset and processes the charge, timing, and event information to calculate the timing resolution and mean charge in photoelectrons. + + The `finish()` method is the main entry point for this tool. It reads the output data from the HDF5 file, filters the data to remove cosmic ray events, and then calculates the timing resolution and mean charge per photoelectron. The timing resolution is calculated using a weighted mean and variance approach, with an option to use a bootstrapping method to estimate the error on the RMS value. + + The method returns the RMS of the timing resolution, the error on the RMS, and the mean charge in photoelectrons. + """ + name = "TimingResolutionTestTool" + + componentsList = ComponentNameList( + NectarCAMComponent, + default_value=["ToMComp"], + help="List of Component names to be apply, the order will be respected", + ).tag(config=True) + + + + + def _init_output_path(self): + if self.max_events is None: + filename = f"{self.name}_run{self.run_number}.h5" + else: + filename = f"{self.name}_run{self.run_number}_maxevents{self.max_events}.h5" + self.output_path = pathlib.Path( + f"{os.environ.get('NECTARCAMDATA','/tmp')}/tests/{filename}" + ) + + def finish(self, bootstrap=False, *args, **kwargs): + + super().finish(return_output_component=True, *args, **kwargs) + + #tom_mu_all= output[0].tom_mu + #tom_sigma_all= output[0].tom_sigma + output_file = h5py.File(self.output_path) + + charge_all = [] + tom_no_fit_all = [] + good_evts = [] + + for thing in output_file: + group = output_file[thing] + dataset = group['ToMContainer'] + data = dataset[:] + #print("data",data) + for tup in data: + try: + npixels = tup[1] + charge_all.extend(tup[6]) + tom_no_fit_all.extend(tup[7]) + good_evts.extend(tup[8]) + except: + break + + output_file.close() + + + tom_no_fit_all= np.array(tom_no_fit_all) + charge_all = np.array(charge_all) + #print(tom_no_fit_all) + #print(charge_all) + + + + #clean cr events + good_evts = np.array(good_evts) + #print(good_evts) + charge=charge_all[good_evts] + mean_charge_pe = np.mean(np.mean(charge,axis=0))/58. + #tom_mu = np.array(tom_mu_all[good_evts]).reshape(len(good_evts),output[0].npixels) + + #tom_sigma = np.array(tom_sigma_all[good_evts]).reshape(len(good_evts),output[0].npixels) + tom_no_fit = np.array(tom_no_fit_all[good_evts]).reshape(len(good_evts),npixels) + #print(tom_no_fit) + #print(tom_no_fit) + + #rms_mu = np.zeros(output[0].npixels) + rms_no_fit = np.zeros(npixels) + + #rms_mu_err = np.zeros(output[0].npixels) + rms_no_fit_err = np.zeros(npixels) + + + #bootstrapping method + + for pix in range(npixels): + + for tom, rms, err in (zip([tom_no_fit[:,pix]],[rms_no_fit],[rms_no_fit_err])): + tom_pos = tom[tom>20] + boot_rms = [] + + sample = tom_pos[tom_pos<32] + #print(sample) + bins = np.linspace(20,32,50) + hist_values, bin_edges = np.histogram(sample, bins=bins) + + + # Compute bin centers + bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2 + + + + if bootstrap: + #print(len(sample)) + if len(sample)!=0: + for _ in range(1000): + bootsample = np.random.choice(sample,size=int(3/4*(len(sample))), replace=True) + # print(len(bootsample), bootsample.mean(), bootsample.std()) + boot_rms.append(bootsample.std()) + # simulated mean of rms + bootrms_mean = np.mean(boot_rms) + # simulated standard deviation of rms + bootrms_std = np.std(boot_rms) + else: + bootrms_std = 0 + bootrms_mean = 0 + # print(bootrms_std) + err[pix]=bootrms_std + rms[pix] = bootrms_mean + + else: + + try: + + + weighted_mean = np.average(bin_centers, weights=hist_values) + #print("Weighted Mean:", weighted_mean) + + # Compute weighted variance + weighted_variance = np.average((bin_centers - weighted_mean)**2, weights=hist_values) + #print("Weighted Variance:", weighted_variance) + + # Compute RMS value (Standard deviation) + rms[pix] = np.sqrt(weighted_variance) + #print("RMS:", rms[pix]) + + # Compute the total number of data points (sum of histogram values, i.e. N) + N = np.sum(hist_values) + #print("Total number of events (N):", N) + + # Error on the standard deviation + err[pix] = rms[pix] / np.sqrt(2 * N) + #print("Error on RMS:", err[pix]) + except: + #no data + rms[pix] = np.nan + err[pix] = np.nan + + + + return rms_no_fit, rms_no_fit_err, mean_charge_pe + + + + +class ToMPairsTool(EventsLoopNectarCAMCalibrationTool): + """ + This class, `ToMPairsTool`, is an `EventsLoopNectarCAMCalibrationTool` that is used to process ToM (Time of maximum) data from NectarCAM. + + The `finish` method has the following functionalities: + + - It reads in ToM data from an HDF5 file and applies a transit time correction to the ToM values using a provided lookup table. + - It calculates the time difference between ToM pairs for both corrected and uncorrected ToM values. + - It returns the uncorrected ToM values, the corrected ToM values, the pixel IDs, and the time difference calculations for the uncorrected and corrected ToM values. + + The class has several configurable parameters, including the list of NectarCAM components to apply, the maximum number of events to process, and the output file path. + + """ + + name = "ToMPairsTool" + + componentsList = ComponentNameList( + NectarCAMComponent, + default_value=["ToMComp"], + help="List of Component names to be apply, the order will be respected", + ).tag(config=True) + + + + + + def finish(self, *args, **kwargs): + + super().finish(return_output_component=True, *args[1:], **kwargs) + + tt_path = args[0] + pmt_tt = pd.read_csv(tt_path)["tom_pmttt_delta_correction"].values + + tom_no_fit_all= [] + + pixels_id = [] + + + output_file = h5py.File(self.output_path) + + + for thing in output_file: + group = output_file[thing] + dataset = group['ToMContainer'] + data = dataset[:] + #print("data",data) + for tup in data: + try: + pixels_id.extend(tup[2]) + tom_no_fit_all.extend(tup[7]) + except: + break + + output_file.close() + pixels_id = np.array(pixels_id) + tom_no_fit_all = np.array(tom_no_fit_all) + + #clean cr events + #good_evts = output[0].good_evts + #charge=charge_all[good_evts] + #mean_charge_pe = np.mean(np.mean(charge,axis=0))/58. + tom_no_fit = np.array(tom_no_fit_all,dtype=np.float64) #tom(event,pixel) + #tom_no_fit = tom_no_fit[np.all(tom_no_fit>0,axis=0)] + tom_corrected = -np.ones(tom_no_fit.shape,dtype=np.float128) #-1 for the ones that have tom beyond 0-60 + + iter = enumerate(pixels_id) + + for i,pix in iter: + #print(pix, pmt_tt[pix]) + normal_values = [a and b for a,b in zip(tom_no_fit[:,i]>0,tom_no_fit[:,i]<60)] + + tom_corrected[normal_values,i] = tom_no_fit[:,i][normal_values]-pmt_tt[pix] + + #print(tom_corrected) + + pixel_ind = [i for i in range(len(pixels_id))] #dealing with indices of pixels in array + pixel_pairs = list(combinations(pixel_ind,2)) + dt_no_correction = np.zeros((len(pixel_pairs),tom_no_fit_all.shape[0])) + dt_corrected = np.zeros((len(pixel_pairs),tom_no_fit_all.shape[0])) + + for i,pair in enumerate(pixel_pairs): + pix1_ind = pixel_ind[pair[0]] + pix2_ind = pixel_ind[pair[1]] + + for event in range(tom_no_fit_all.shape[0]): + cond_no_correction = (tom_no_fit[event,pix1_ind]>0 and tom_no_fit[event,pix1_ind]<60 + and tom_no_fit[event,pix2_ind]>0 and tom_no_fit[event,pix2_ind]<60) + cond_correction = (tom_corrected[event,pix1_ind]>0 and tom_corrected[event,pix1_ind]<60 + and tom_corrected[event,pix2_ind]>0 and tom_corrected[event,pix2_ind]<60) + + if cond_no_correction: #otherwise will be nan + + dt_no_correction[i,event] = tom_no_fit[event,pix1_ind] - tom_no_fit[event,pix2_ind] + + else: + dt_no_correction[i,event] = np.nan + + if cond_correction: + + dt_corrected[i,event] = tom_corrected[event,pix1_ind] - tom_corrected[event,pix2_ind] + + else: + dt_corrected[i,event] = np.nan + + + + + # rms_no_fit = np.zeros(output[0].npixels) + # rms_no_fit_err = np.zeros(output[0].npixels) + + + return tom_no_fit,tom_corrected,pixels_id, dt_no_correction, dt_corrected, pixel_pairs + + + +class PedestalContainer(NectarCAMContainer): + + """ + Attributes of the PedestalContainer class that store various data related to the pedestal of a NectarCAM event. + + Attributes: + run_number (np.uint16): The run number associated with the waveforms. + npixels (np.uint16): The number of effective pixels. + pixels_id (np.ndarray[np.uint16]): The IDs of the pixels. + ucts_timestamp (np.ndarray[np.uint64]): The UCTS timestamp of the events. + event_type (np.ndarray[np.uint8]): The trigger event type. + event_id (np.ndarray[np.uint32]): The event IDs. + pedestal_hg (np.ndarray[np.float64]): The high gain pedestal per event. + pedestal_lg (np.ndarray[np.float64]): The low gain pedestal per event. + rms_ped_hg (np.ndarray[np.float64]): The high gain pedestal RMS per event. + rms_ped_lg (np.ndarray[np.float64]): The low gain pedestal RMS per event. + """ + run_number = Field( + type=np.uint16, + description="run number associated to the waveforms", + ) + npixels = Field( + type=np.uint16, + description="number of effective pixels", + ) + pixels_id = Field(type=np.ndarray, dtype=np.uint16, ndim=1, description="pixel ids") + ucts_timestamp = Field( + type=np.ndarray, dtype=np.uint64, ndim=1, description="events ucts timestamp" + ) + event_type = Field( + type=np.ndarray, dtype=np.uint8, ndim=1, description="trigger event type" + ) + event_id = Field(type=np.ndarray, dtype=np.uint32, ndim=1, description="event ids") + + + pedestal_hg = Field( + type=np.ndarray, dtype=np.float64, ndim=2, description="High gain pedestal per event" + ) + + pedestal_lg = Field( + type=np.ndarray, dtype=np.float64, ndim=2, description="Low gain pedestal per event" + ) + + rms_ped_hg = Field( + type=np.ndarray, dtype=np.float64, ndim=2, description="High gain pedestal rms per event" + ) + + rms_ped_lg = Field( + type=np.ndarray, dtype=np.float64, ndim=2, description="Low gain pedestal rms per event" + ) + + + + +class PedestalComp(NectarCAMComponent): + + """ + The `PedestalComp` class is a NectarCAMComponent that is responsible for processing the pedestal and RMS of the high and low gain waveforms for each event. + + The `__init__` method initializes the `PedestalComp` class. It sets up several member variables to store pedestal related data such as timestamps, event types, event IDs, pedestal and pedestal rms values for both gains. + + The `__call__` method is called for each event, and it processes the waveforms to calculate the pedestal and RMS for the high and low gain channels. The results are stored in the class attributes `__pedestal_hg`, `__pedestal_lg`, `__rms_ped_hg`, and `__rms_ped_lg`. + + The `finish` method is called at the end of processing, and it returns a `PedestalContainer` object containing the calculated pedestal and RMS values, as well as other event information. + """ + + def __init__(self, subarray, config=None, parent=None, *args, **kwargs): + super().__init__( + subarray=subarray, config=config, parent=parent, *args, **kwargs + ) + ## If you want you can add here members of MyComp, they will contain interesting quantity during the event loop process + + self.__ucts_timestamp = [] + self.__event_type = [] + self.__event_id = [] + + + self.__pedestal_hg = [] + + self.__pedestal_lg = [] + + self.__rms_ped_hg = [] + self.__rms_ped_lg = [] + + + + + + ##This method need to be defined ! + def __call__(self, event: NectarCAMDataContainer, *args, **kwargs): + self.__event_id.append(np.uint32(event.index.event_id)) + + + #print(event.trigger.event_type) + + self.__event_type.append(event.trigger.event_type.value) + self.__ucts_timestamp.append(event.nectarcam.tel[0].evt.ucts_timestamp) + + + wfs = [] + wfs.append(event.r0.tel[0].waveform[constants.HIGH_GAIN][self.pixels_id]) + wfs.append(event.r0.tel[0].waveform[constants.LOW_GAIN][self.pixels_id]) + + + + + #####THE JOB IS HERE###### + + + for i, (pedestal,rms_pedestal) in enumerate(zip([self.__pedestal_hg,self.__pedestal_lg],[self.__rms_ped_hg,self.__rms_ped_lg])): + wf = np.array(wfs[i]) #waveform per gain + + ped = np.array([np.mean(wf[pix]) for pix in range(len(self.pixels_id))]) + rms_ped = np.array([np.std(wf[pix]) for pix in range(len(self.pixels_id))]) + + ped[ped>1000]=np.nan + rms_ped[ped>1000]=np.nan + pedestal.append(ped) + rms_pedestal.append(rms_ped) + + + + + + + + + + ##This method need to be defined ! + def finish(self): + + output = PedestalContainer( + run_number=PedestalContainer.fields["run_number"].type(self._run_number), + npixels=PedestalContainer.fields["npixels"].type(self._npixels), + pixels_id=PedestalContainer.fields["pixels_id"].dtype.type(self.pixels_id), + ucts_timestamp=PedestalContainer.fields["ucts_timestamp"].dtype.type( + self.__ucts_timestamp + ), + event_type=PedestalContainer.fields["event_type"].dtype.type(self.__event_type), + event_id=PedestalContainer.fields["event_id"].dtype.type(self.__event_id), + + pedestal_hg=PedestalContainer.fields["pedestal_hg"].dtype.type( + self.__pedestal_hg + ), + + pedestal_lg=PedestalContainer.fields["pedestal_lg"].dtype.type( + self.__pedestal_lg + ), + rms_ped_hg=PedestalContainer.fields["rms_ped_hg"].dtype.type( + self.__rms_ped_hg + ), + + rms_ped_lg=PedestalContainer.fields["rms_ped_lg"].dtype.type( + self.__rms_ped_lg + ), + + + ) + return output + + + + + + +class PedestalTool(EventsLoopNectarCAMCalibrationTool): + """ + This class is a part of the PedestalTool, which is an EventsLoopNectarCAMCalibrationTool. + + The finish() method opens the output file, which is an HDF5 file, and extracts the `rms_ped_hg` (root mean square of the high gain pedestal) values from the `PedestalContainer` dataset. Finally, it closes the output file and returns the list of `rms_ped_hg` values. + + This method is used to post-process the output of the PedestalTool and extract specific information from the generated HDF5 file. + """ + name = "PedestalTool" + + componentsList = ComponentNameList( + NectarCAMComponent, + default_value=["PedestalComp"], + help="List of Component names to be apply, the order will be respected", + ).tag(config=True) + + + + + def _init_output_path(self): + if self.max_events is None: + filename = f"{self.name}_run{self.run_number}.h5" + else: + filename = f"{self.name}_run{self.run_number}_maxevents{self.max_events}.h5" + self.output_path = pathlib.Path( + f"{os.environ.get('NECTARCAMDATA','/tmp')}/tests/{filename}" + ) + + + + + def finish(self, *args, **kwargs): + + super().finish(return_output_component=False, *args, **kwargs) + #print(self.output_path) + output_file = h5py.File(self.output_path) + + rms_ped_hg = [] + + for thing in output_file: + group = output_file[thing] + dataset = group['PedestalContainer'] + data = dataset[:] + #print("data",data) + for tup in data: + try: + rms_ped_hg.extend(tup[8]) + except: + break + + output_file.close() + + + + return rms_ped_hg + + + +class UCTSContainer(NectarCAMContainer): + """ + Defines the fields for the UCTSContainer class, which is used to store various data related to UCTS events. + + The fields include: + - `run_number`: The run number associated with the waveforms. + - `npixels`: The number of effective pixels. + - `pixels_id`: The IDs of the pixels. + - `ucts_timestamp`: The UCTS timestamp of the events. + - `event_type`: The trigger event type. + - `event_id`: The IDs of the events. + - `ucts_busy_counter`: The UCTS busy counter. + - `ucts_event_counter`: The UCTS event counter. + """ + run_number = Field( + type=np.uint16, + description="run number associated to the waveforms", + ) + npixels = Field( + type=np.uint16, + description="number of effective pixels", + ) + pixels_id = Field(type=np.ndarray, dtype=np.uint16, ndim=1, description="pixel ids") + ucts_timestamp = Field( + type=np.ndarray, dtype=np.uint64, ndim=1, description="events ucts timestamp" + ) + event_type = Field( + type=np.ndarray, dtype=np.uint8, ndim=1, description="trigger event type" + ) + event_id = Field(type=np.ndarray, dtype=np.uint32, ndim=1, description="event ids") + ucts_busy_counter = Field(type=np.ndarray, dtype=np.uint32, ndim=1, description="busy counter") + ucts_event_counter = Field(type=np.ndarray, dtype=np.uint32, ndim=1, description="event counter") + + + +class UCTSComp(NectarCAMComponent): + + + """ + The `__init__` method initializes the `UCTSComp` class, which is a NectarCAMComponent. It sets up several member variables to store UCTS related data, such as timestamps, event types, event IDs, busy counters, and event counters. + + The `__call__` method is called for each event, and it appends the UCTS-related data from the event to the corresponding member variables. + + The `finish` method creates and returns a `UCTSContainer` object, which is a container for the UCTS-related data that was collected during the event loop. + """ + def __init__(self, subarray, config=None, parent=None, *args, **kwargs): + super().__init__( + subarray=subarray, config=config, parent=parent, *args, **kwargs + ) + ## If you want you can add here members of MyComp, they will contain interesting quantity during the event loop process + + self.__ucts_timestamp = [] + self.__event_type = [] + self.__event_id = [] + self.__ucts_busy_counter = [] + self.__ucts_event_counter = [] + + + + + + + + ##This method need to be defined ! + def __call__(self, event: NectarCAMDataContainer, *args, **kwargs): + self.__event_id.append(np.uint32(event.index.event_id)) + self.__event_type.append(event.trigger.event_type.value) + self.__ucts_timestamp.append(event.nectarcam.tel[0].evt.ucts_timestamp) + self.__ucts_busy_counter.append(event.nectarcam.tel[0].evt.ucts_busy_counter) + self.__ucts_event_counter.append(event.nectarcam.tel[0].evt.ucts_event_counter) + + + + + + + + + ##This method need to be defined ! + def finish(self): + + output = UCTSContainer( + run_number=UCTSContainer.fields["run_number"].type(self._run_number), + npixels=UCTSContainer.fields["npixels"].type(self._npixels), + pixels_id=UCTSContainer.fields["pixels_id"].dtype.type(self.pixels_id), + ucts_timestamp=UCTSContainer.fields["ucts_timestamp"].dtype.type( + self.__ucts_timestamp + ), + ucts_busy_counter=UCTSContainer.fields["ucts_busy_counter"].dtype.type( + self.__ucts_busy_counter + ), + ucts_event_counter=UCTSContainer.fields["ucts_event_counter"].dtype.type( + self.__ucts_event_counter + ), + event_type=UCTSContainer.fields["event_type"].dtype.type(self.__event_type), + event_id=UCTSContainer.fields["event_id"].dtype.type(self.__event_id), + + ) + return output + + + + + +class DeadtimeTestTool(EventsLoopNectarCAMCalibrationTool): + + """ + The `DeadtimeTestTool` class is an `EventsLoopNectarCAMCalibrationTool` that is used to test the deadtime of NectarCAM. + + The `finish` method is responsible for reading the data from the HDF5 file, extracting the relevant information (UCTS timestamps, event counters, and busy counters), and calculating the deadtime-related metrics. The method returns the UCTS timestamps, the time differences between consecutive UCTS timestamps, the event counters, the busy counters, the collected trigger rate, the total time, and the deadtime percentage. + """ + name = "DeadtimeTestTool" + + componentsList = ComponentNameList( + NectarCAMComponent, + default_value=["UCTSComp"], + help="List of Component names to be apply, the order will be respected", + ).tag(config=True) + + + + def finish(self, *args, **kwargs): + + super().finish(return_output_component=False, *args, **kwargs) + #print(self.output_path) + output_file = h5py.File(self.output_path) + + ucts_timestamps=[] + event_counter=[] + busy_counter=[] + + for thing in output_file: + group = output_file[thing] + dataset = group['UCTSContainer'] + data = dataset[:] + #print("data",data) + for tup in data: + try: + ucts_timestamps.extend(tup[3]) + event_counter.extend(tup[7]) + busy_counter.extend(tup[6]) + except: + break + #print(output_file.keys()) + #tom_mu_all= output[0].tom_mu + #tom_sigma_all= output[0].tom_sigma + #ucts_timestamps= np.array(output_file["ucts_timestamp"]) + ucts_timestamps = np.array(ucts_timestamps).flatten() + #print(ucts_timestamps) + event_counter = np.array(event_counter).flatten() + busy_counter = np.array(busy_counter).flatten() + #print(ucts_timestamps) + delta_t = [ucts_timestamps[i]-ucts_timestamps[i-1] for i in range(1,len(ucts_timestamps))] + # event_counter = np.array(output_file['ucts_event_counter']) + # busy_counter=np.array(output_file['ucts_busy_counter']) + output_file.close() + + time_tot = ((ucts_timestamps[-1] - ucts_timestamps[0]) * u.ns).to(u.s) + collected_trigger_rate = (event_counter[-1]+busy_counter[-1])/time_tot + deadtime_pc = busy_counter[-1]/(event_counter[-1]+busy_counter[-1])*100 + + + + return ucts_timestamps,delta_t,event_counter,busy_counter,collected_trigger_rate, time_tot, deadtime_pc + + diff --git a/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/utils.py b/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/utils.py new file mode 100644 index 00000000..5eb7c50b --- /dev/null +++ b/src/nectarchain/user_scripts/dmousadi/test_scripts/tests/utils.py @@ -0,0 +1,500 @@ +from ctapipe.utils import get_dataset_path + +import numpy as np +import matplotlib.pyplot as plt + + +from scipy.interpolate import interp1d +from traitlets.config import Config +from scipy.stats import expon, poisson +from lmfit.models import Model +from astropy import units as u +from iminuit import Minuit +import os + + + +config = Config(dict(NectarCAMEventSource=dict( + NectarCAMR0Corrections=dict( + calibration_path=None, + apply_flatfield=False, + select_gain=False, + )))) + +#(filter, optical density, transmission) +#configurations with the same order as in the document +filters=np.array( + [ + [0.00, 0.15, 0.30, 0.45, 0.60, 0.75, 0.90, 1.00, 1.15, 1.30, 1.30, 1.45, 1.50, 1.60, 1.65, 1.80, 1.90, 2.00, 2.10, 2.15, 2.30, 2.30, 2.50, 2.60, 2.80, 3.00, 3.30, 3.50], + [0.000000, 0.213853, 0.370538, 0.584391, 0.739449, 0.953302, 1.109987, 1.344491, 1.558344, 1.715029, 1.792931, 2.006784, 2.089531, 2.163469, 2.303384, 2.460068, 2.532380, 2.695312, 2.828980, 2.909165, 3.065849, 3.137422, 3.434022, 3.434761, 3.882462, 4.039803, 4.488243, 4.784842], + [1.000000, 0.611149, 0.426052, 0.260381, 0.182201, 0.111352, 0.077627, 0.045239, 0.027647, 0.019274, 0.016109, 0.009845, 0.008137, 0.006863, 0.004973, 0.003467, 0.002935, 0.002017, 0.001483, 0.001233, 0.000859, 0.000729, 0.000368, 0.000367, 0.000131, 0.000091, 0.000032, 0.000016] + ] +).T + +#the next ones are in the order of the test linearity runs by federica +#maybe more practical because the runs will be in that order again +trasmission_390ns = np.array([0.002016919, +0.045238588, +0.0276475, +0.019273981, +0.008242516, +0.000368111, +9.12426E-05, +1, +0.611148605, +0.260380956, +0.111351889, +0.004972972, +0.001232637, +0.009844997, +0.426051788, +0.077627064, +0.006863271, +0.003466823, +0.000859312, +0.016109006, +0.182201004, +0.002935077, +0.001482586, +0.000367485, +0.008137092, +0.00013108, +1.64119E-05, +3.24906E-05, +0.000728749,]) + + +optical_density_390ns = np.array([2.69531155, +1.34449096, +1.55834414, +1.71502857, +2.08394019, +3.43402173, +4.03980251, +0.00000000, +0.21385318, +0.58439078, +0.95330241, +2.30338395, +2.90916473, +2.00678442, +0.37053761, +1.10998684, +2.16346885, +2.46006838, +3.06584916, +1.79293125, +0.73944923, +2.53238048, +2.82898001, +3.43476079, +2.08953077, +3.88246202, +4.78484233, +4.48824280, +3.13742221,]) + +adc_to_pe = 58. + +plot_parameters = { + "High Gain": { + "linearity_range": [5,-5], + "text_coords":[0.1, 600], + "label_coords":[20, 400], + "color":"C0", + "initials":"HG" + }, + "Low Gain": { + "linearity_range": [11,-1], + "text_coords":[0.3e2, 5], + "label_coords":[1.5e3, 9], + "color":"C4", + "initials":"LG" + } +} + +intensity_percent = np.array([13. , 15. , 16. , 16.5, + 20.6, 22. , + 25. , 35. , 33]) +intensity_to_charge =np.array([ 1.84110824, 2.09712394, 2.24217532, 2.37545181, + 10.32825426, 28.80655155, + 102.79668643, 191.47895686, 188]) + + +source_ids_deadtime = [0 for i in range(3332,3342)] + [1 for i in range(3342,3350)] + [2 for i in range(3552,3562)] + +deadtime_labels = {0:{'source':"random generator",'color':'blue'}, + 1:{'source':"nsb source", 'color':'green'}, + 2:{'source':"laser",'color':'red'} +} + + +def pe_from_intensity_percentage (percent,percent_from_calibration=intensity_percent,known_charge=intensity_to_charge): + """ + Converts a percentage of intensity to the corresponding charge value based on a known calibration. + + Args: + percent (numpy.ndarray): The percentage of intensity to convert to charge. + percent_from_calibration (numpy.ndarray, optional): The known percentages of intensity used in the calibration. Defaults to `intensity_percent`. + known_charge (numpy.ndarray, optional): The known charge values corresponding to the calibration percentages. Defaults to `intensity_to_charge`. + + Returns: + numpy.ndarray: The charge values corresponding to the input percentages. + """ + #known values from laser calibration + + #runs are done with intensity of 15-35 percent + f = interp1d(percent_from_calibration, known_charge) + charge = np.zeros(len(percent)) + for i in range(len(percent)): + charge[i]=f(percent[i]) + + return charge + +#functions by federica +def fit_function(x,a,b): + """ + Computes a linear function of the form `a*x + b`. + + Args: + x (float): The input value. + a (float): The slope coefficient. + b (float): The y-intercept. + + Returns: + float: The result of the linear function. + """ + return a*x + b + + +def fit_function2(x,a,b,c): + """ + Computes a quadratic function of the form `a*(x**2) + b*x + c`. + + Args: + x (float): The input value. + a (float): The coefficient of the squared term. + b (float): The coefficient of the linear term. + c (float): The constant term. + + Returns: + float: The result of the quadratic function. + """ + return a*(x**2) + b*x + c + +def fit_function3(x,a,b,c, d): + """ + Computes a function of the form `(a*x + b)/(1+c) + d`. + + Args: + x (float): The input value. + a (float): The coefficient of the linear term. + b (float): The constant term in the numerator. + c (float): The coefficient of the denominator. + d (float): The constant term added to the result. + + Returns: + float: The result of the function. + """ + return (a*x + b)/(1+c) + d + +def fit_function_hv(x,a,b): + """ + Computes a function of the form `a/sqrt(x) + b`. + + Args: + x (float): The input value. + a (float): The coefficient of the term `1/sqrt(x)`. + b (float): The constant term. + + Returns: + float: The result of the function. + """ + return a/np.sqrt(x)+b + + +def err_ratio(nominator, denominator, err_norm, err_denom, cov_nom_den = 0): + """ + Computes the error ratio for a given nominator, denominator, and their respective errors. + + Args: + nominator (float): The nominator value. + denominator (float): The denominator value. + err_norm (float): The error of the nominator. + err_denom (float): The error of the denominator. + cov_nom_den (float, optional): The covariance between the nominator and denominator. Defaults to 0. + + Returns: + float: The error ratio. + """ + delta_err2 = (err_norm/nominator)**2 + (err_denom/denominator)**2 -2*cov_nom_den/(nominator*denominator) + ratio = nominator/denominator + return np.sqrt(delta_err2)*ratio + +def err_sum(err_a, err_b, cov_a_b=0): + """ + Computes the square root of the sum of the squares of `err_a` and `err_b`, plus twice the covariance `cov_a_b`. + + This function is used to calculate the combined error of two values, taking into account their individual errors and the covariance between them. + + Args: + err_a (float): The error of the first value. + err_b (float): The error of the second value. + cov_a_b (float, optional): The covariance between the two values. Defaults to 0. + + Returns: + float: The combined error. + """ + return np.sqrt(err_a**2 + err_b**2 + 2*cov_a_b) + +def closest_value(input_list, input_value): + """ + Returns the value in `input_list` that is closest to `input_value`. + + Args: + input_list (list): The list of values to search. + input_value (float): The value to find the closest match for. + + Returns: + float: The value in `input_list` that is closest to `input_value`. + """ + arr = np.asarray(input_list) + i = (np.abs(arr - input_value)).argmin() + return arr[i] + +def gaus(x,a,mu,sigma): + """ + Computes the Gaussian function with the given parameters. + + Args: + x (float): The input value. + a (float): The amplitude of the Gaussian function. + mu (float): The mean (center) of the Gaussian function. + sigma (float): The standard deviation of the Gaussian function. + + Returns: + float: The value of the Gaussian function at the given input `x`. + """ + return a*np.exp(-(x-mu)**2/(2*sigma**2)) + + +#from stackoverflow +def argmedian(x, axis=None): + """ + Returns the index of the median element in the input array `x` along the specified axis. + + If `axis` is `None`, the function returns the index of the median element in the flattened array. + Otherwise, it computes the argmedian along the specified axis and returns an array of indices. + + Args: + x (numpy.ndarray): The input array. + axis (int or None, optional): The axis along which to compute the argmedian. If `None`, the argmedian is computed on the flattened array. + + Returns: + int or numpy.ndarray: The index or indices of the median element(s) in the input array. + """ + if axis is None: + return np.argpartition(x, len(x) // 2)[len(x) // 2] + else: + # Compute argmedian along specified axis + return np.apply_along_axis( + lambda x: np.argpartition(x, len(x) // 2)[len(x) // 2], + axis=axis, arr=x + ) + + + + +def pe2photons(x): + """ + Converts the input value `x` from photons to photoelectrons (PE) by multiplying it by 4. + + Args: + x (float): The input value in photons. + + Returns: + float: The input value converted to photoelectrons. + """ + return x * 4 + + +def photons2pe(x): + """ + Converts the input value `x` from photoelectrons (PE) to photons by dividing it by 4. + + Args: + x (float): The input value in photoelectrons. + + Returns: + float: The input value converted to photons. + """ + return x / 4 + + +#from federica's notebook +class ExponentialFitter(): + + """ + Represents an exponential fitter class that computes the expected distribution and the minus 2 log likelihood for a given dataset and exponential parameters. + + Attributes: + datas (numpy.ndarray): The input data array. + bin_edges (numpy.ndarray): The bin edges for the data. + + Methods: + compute_expected_distribution(norm, loc, scale): + Computes the expected distribution given the normalization, location, and scale parameters. + expected_distribution(x): + Returns the expected distribution given the parameters in `x`. + compute_minus2loglike(x): + Computes the minus 2 log likelihood given the parameters in `x`. + """ + def __init__(self,datas,bin_edges): + self.datas = datas.copy() + self.bin_edges = bin_edges.copy() + + def compute_expected_distribution(self,norm,loc,scale): + cdf_low = expon.cdf(self.bin_edges[:-1],loc=loc,scale=scale) + cdf_up = expon.cdf(self.bin_edges[1:],loc=loc,scale=scale) + delta_cdf = cdf_up - cdf_low + return norm*delta_cdf + + def expected_distribution(self,x): + return self.compute_expected_distribution(x[0],x[1],x[2]) + + def compute_minus2loglike(self,x): + norm = x[0] + loc = x[1] + scale = x[2] + + expected = self.compute_expected_distribution(norm,loc,scale) + chi2_mask = expected>0. + minus2loglike = -2.*np.sum( poisson.logpmf(self.datas[chi2_mask],mu=expected[chi2_mask]) ) + minus2loglike0 = -2.*np.sum( poisson.logpmf(self.datas[chi2_mask],mu=self.datas[chi2_mask]) ) +# print(f'{minus2loglike0 = }') + return minus2loglike - minus2loglike0 + + +def pois(x, A, R): + """ + Computes the expected distribution for a Poisson process with rate parameter `R`. + + Args: + x (float): The input value. + A (float): The amplitude parameter. + R (float): The rate parameter. + + Returns: + float: The expected distribution for the Poisson process. + """ + '''poisson function, parameter r (rate) is the fit parameter''' + return A*np.exp(x*R) + + +def deadtime_and_expo_fit(time_tot, deadtime_us,run, output_plot = None): + """ + Computes the deadtime and exponential fit parameters for a given dataset. + + Args: + time_tot (float): The total time of the dataset. + deadtime_us (float): The deadtime of the dataset in microseconds. + run (int): The run number. + output_plot (str, optional): The path to save the output plot. + + Returns: + tuple: A tuple containing the following values: + - deadtime (float): The deadtime of the dataset. + - deadtime_bin (float): The bin edge corresponding to the deadtime. + - deadtime_err (float): The error on the deadtime. + - deadtime_bin_length (float): The length of the deadtime bin. + - total_delta_t_for_busy_time (float): The total time of the dataset. + - parameter_A_new (float): The amplitude parameter of the exponential fit. + - parameter_R_new (float): The rate parameter of the exponential fit. + - parameter_A_err_new (float): The error on the amplitude parameter. + - parameter_R_err_new (float): The error on the rate parameter. + - first_bin_length (float): The length of the first bin. + - tot_nr_events_histo (int): The total number of events in the histogram. + """ + #function by federica + + total_delta_t_for_busy_time = time_tot + data = deadtime_us + pucts = np.histogram(data, bins=100, range=(1e-2, 50)) + + deadtime = (min(data[~np.isnan(data)])) + + first_nonemptybin = np.where(pucts[0] > 0)[0][0] + deadtime_err = (pucts[1][first_nonemptybin] - pucts[1][first_nonemptybin - 1]) /\ + np.sqrt(pucts[0][first_nonemptybin]) + deadtime_bin = (pucts[1][first_nonemptybin]) + deadtime_bin_length = (pucts[1][first_nonemptybin] - pucts[1][first_nonemptybin - 1]) + + # the bins should be of integer width, because poisson is an integer distribution + bins = np.arange(100000) - 0.5 + entries, bin_edges = np.histogram(data, bins=bins, range=[0, total_delta_t_for_busy_time]) + bin_middles = 0.5 * (bin_edges[1:] + bin_edges[:-1]) + + first_bin_length = (bin_edges[1] - bin_edges[0]) + tot_nr_events_histo = (np.sum(entries)) + + # second fit + model = Model(pois) + params = model.make_params(A=2e3, R=-0.001) + result = model.fit(entries, params, x=bin_middles) + #print('**** FIT RESULTS RUN {} ****'.format(run)) + #print(result.fit_report()) + #print('chisqr = {}, redchi = {}'.format(result.chisqr, result.redchi)) + + parameter_A_new = (result.params['A'].value) + parameter_R_new = (-1 * result.params['R'].value) + parameter_A_err_new = (result.params['A'].stderr) + parameter_R_err_new = (result.params['R'].stderr) + + plt.close() + + if output_plot: + + fig, ax = plt.subplots(1, 1, figsize=(10, 10 / 1.61)) + plt.hist(data, bins= np.arange(1000) - 0.5, + alpha=0.8, + histtype='step', density=0, lw=2, label='Run {}'.format(run)) + + # # plot poisson-deviation with fitted parameter + x_plot = np.arange(0, 1000) + + + plt.plot( + x_plot, + pois(x_plot, result.params['A'].value, result.params['R'].value), + marker='', linestyle='-', lw=3, color='C3', alpha=0.85, + label=r'Fit result: %2.3f $\exp^{-{x}/({%2.0f}~\mu\mathrm{s})}$' % ( + result.params['A'].value, abs(1 / result.params['R'].value)), + ) + + R = ((-1 * result.params['R'].value) * (1 / u.us)).to(u.kHz).value + + + R_stderr = ((result.params['R'].stderr) * (1 / u.us)).to(u.kHz).value + + ax.text(600, entries[1] / 1 , '$y = A \cdot \exp({-R \cdot x})$\n' + # + r'$A=%2.2f \pm %2.2f$'%(as_si((result.params['A'].value/1000)*1e3,2), as_si((result.params['A'].stderr/1000)*1e3,2)) + + r'$A=%2.2f \pm %2.2f$' % (result.params['A'].value, result.params['A'].stderr) + + '\n' + r'$R=(%2.2f \pm %2.2f)$ kHz' % (R, R_stderr) + + '\n' + r'$\chi^2_\nu = %2.2f$' % ((result.redchi)) + + '\n' + r'$\delta_{\mathrm{deadtime}} = %2.3f \, \mu$s' % (deadtime), + backgroundcolor='white', bbox=dict(facecolor='white', edgecolor='C3', lw=2, + boxstyle='round,pad=0.3'), + ha='left', va='top', fontsize=17, color='k', alpha=0.9) + + + plt.xlabel(r'$\Delta$T [$\mu$s]') + plt.ylabel('Entries') + plt.title('Run {}'.format(run), y=1.02) + plt.savefig(os.path.join(output_plot,"deadtime_and_expo_fit_{}.png".format(run))) + + + return deadtime, deadtime_bin, deadtime_err, deadtime_bin_length, \ + total_delta_t_for_busy_time, parameter_A_new, parameter_R_new, \ + parameter_A_err_new, parameter_R_err_new, \ + first_bin_length, tot_nr_events_histo +