From 73e64ba6581178df4514a9160c670a8835ef2792 Mon Sep 17 00:00:00 2001 From: Luca Baldini Date: Thu, 11 Apr 2024 16:51:09 +0200 Subject: [PATCH 01/52] Initial implementation of new digitization classes. --- hexsample/digi.py | 296 ++++++++++++++++++++++++++++++++++++++++++++- tests/test_digi.py | 67 ++++++++-- 2 files changed, 351 insertions(+), 12 deletions(-) diff --git a/hexsample/digi.py b/hexsample/digi.py index 89c3585..bdf4d8a 100644 --- a/hexsample/digi.py +++ b/hexsample/digi.py @@ -20,6 +20,7 @@ """Digitization facilities. """ +from collections import Counter from dataclasses import dataclass from typing import Tuple @@ -34,6 +35,299 @@ +@dataclass +class DigiEventBase: + + """Base class for a digitized event. + + This includes alll the timing information that is common to the different digitized + event structure, as well as the PHA content of the pixels in the event, but no + information about the physical location on the latter within the readout chip. + It is responsibility of the derived classes to provide this information, in the + way that is more conventient, depending on the particular readout strategy. + + Arguments + --------- + trigger_id : int + The trigger identifier. + + seconds : int + The integer part of the timestamp. + + microseconds : int + The fractional part of the timestamp. + + pha : np.ndarray + The pixel content of the event, in the form of a 1-dimensional array. + """ + + trigger_id: int + seconds: int + microseconds: int + livetime: int + pha: np.ndarray + + def __eq__(self, other) -> bool: + """Overloaded comparison operator. + """ + return (self.trigger_id, self.seconds, self.microseconds, self.livetime) == \ + (other.trigger_id, other.seconds, other.microseconds, other.livetime) and \ + np.allclose(self.pha, other.pha) + + def timestamp(self) -> float: + """Return the timestamp of the event, that is, the sum of the second and + microseconds parts of the DigiEvent contributions as a floating point number. + """ + return self.seconds + 1.e-6 * self.microseconds + + def ascii(self, pha_width: int = 5): + """ + """ + raise NotImplementedError + + + +@dataclass +class DigiEventSparse(DigiEventBase): + + """Sparse digitized event. + + In this particular incarnation of a digitized event we have no ROI structure, + nor any rule as to the shape or morphology of the particular set of pixels + being read out. The event represent an arbitrary collection of pixels, and we + carry over the arrays of their row and column identifiers, in the form of + two arrays whose length must match that of the PHA values. + + Arguments + --------- + columns : np.ndarray + The column identifier of the pixels in the event (must have the same length + of the pha class member). + + rows : np.ndarray + The row identifier of the pixels in the event (must have the same length + of the pha class member). + """ + + columns: np.ndarray + rows: np.ndarray + + def __post_init__(self) -> None: + """Post-initialization code. + """ + if not len(self.rows) == len(self.columns) == len(self.pha): + raise RuntimeError(f'{self.__class__.__name__} has {len(self.rows)} rows' + f', {len(self.columns)} columns, and {len(self.pha)} PHA values') + + def as_dict(self) -> dict: + """Return the pixel content of the event in the form of a {(col, row): pha} + value. + + This is useful to address the pixel content at a given coordinate. We refrain, + for the moment, from implementing a __call__() hook on top of this, as the + semantics would be fundamentally different from that implemented with a + rectangular ROI, where the access happen in ROI coordinates, and not in chip + coordinates, but we might want to come back to this. + """ + return {(col, row): pha for col, row, pha in zip(self.columns, self.rows, self.pha)} + + def ascii(self, pha_width: int = 5) -> str: + """Ascii representation. + """ + pha_dict = self.as_dict() + fmt = f'%{pha_width}d' + cols = np.arange(self.columns.min(), self.columns.max() + 1) + num_cols = cols[-1] - cols[0] + 1 + rows = np.arange(self.rows.min(), self.rows.max() + 1) + big_space = space(2 * pha_width + 1) + text = f'\n{big_space}' + text += ''.join([fmt % col for col in cols]) + text += f'\n{big_space}+{line(pha_width * num_cols)}\n' + for row in rows: + text += f' {fmt % row} |' + for col in cols: + try: + pha = fmt % pha_dict[(col, row)] + except KeyError: + pha = ' ' * pha_width + text += pha + text += f'\n{big_space}|\n' + return text + + + +@dataclass +class DigiEventRectangular(DigiEventBase): + + """ + """ + + roi: RegionOfInterest + + + +@dataclass +class DigiEventCircular(DigiEventBase): + + """ + """ + + row: int + column: int + + + +class HexagonalReadoutBase(HexagonalGrid): + + """Description of a pixel readout chip on a hexagonal matrix. + + Arguments + --------- + layout : HexagonalLayout + The layout of the hexagonal matrix. + + num_cols : int + The number of columns in the readout. + + num_rows : int + The number of rows in the readout. + + pitch : float + The readout pitch in cm. + + enc : float + The equivalent noise charge in electrons. + + gain : float + The readout gain in ADC counts per electron. + """ + + def __init__(self, layout: HexagonalLayout, num_cols: int, num_rows: int, + pitch: float, enc: float, gain: float) -> None: + """Constructor. + """ + # pylint: disable=too-many-arguments + super().__init__(layout, num_cols, num_rows, pitch) + self.enc = enc + self.gain = gain + self.shape = (self.num_rows, self.num_cols) + self.trigger_id = -1 + + + +class HexagonalReadoutSparse(HexagonalReadoutBase): + + """ + """ + + @staticmethod + def zero_suppress(array: np.ndarray, threshold: float) -> None: + """Utility function to zero-suppress an generic array. + + This is returning an array of the same shape of the input where all the + values lower or equal than the zero suppression threshold are set to zero. + + Arguments + --------- + array : array_like + The input array. + + threshold : float + The zero suppression threshold. + """ + array[array <= threshold] = 0 + + @staticmethod + def latch_timestamp(timestamp: float) -> Tuple[int, int, int]: + """Latch the event timestamp and return the corresponding fields of the + digi event contribution: seconds, microseconds and livetime. + + .. warning:: + The livetime calculation is not implemented, yet. + + Arguments + --------- + timestamp : float + The ground-truth event timestamp from the event generator. + """ + microseconds, seconds = np.modf(timestamp) + livetime = 0 + return int(seconds), int(1000000 * microseconds), livetime + + def digitize(self, pha: np.ndarray, zero_sup_threshold: int = 0, + offset: int = 0) -> np.ndarray: + """Digitize the actual signal within a given ROI. + + Arguments + --------- + signal : array_like + The input array of pixel signals to be digitized. + + roi : RegionOfInterest + The target ROI. + + zero_sup_threshold : int + Zero-suppression threshold in ADC counts. + + offset : int + Optional offset in ADC counts to be applied before the zero suppression. + """ + # Add the noise. + if self.enc > 0: + pha += rng.generator.normal(0., self.enc, size=pha.shape) + # ... apply the conversion between electrons and ADC counts... + pha *= self.gain + # ... round to the neirest integer... + pha = np.round(pha).astype(int) + # ... if necessary, add the offset for diagnostic events... + pha += offset + # ... zero suppress the thing... + self.zero_suppress(pha, zero_sup_threshold) + # ... flatten the array to simulate the serial readout and return the + # array as the BEE would have. + return pha.flatten() + + def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: float, + zero_sup_threshold: int = 0, offset: int = 0) -> DigiEventSparse: + """Readout an event. + + Arguments + --------- + timestamp : float + The event timestamp. + + x : array_like + The physical x coordinates of the input charge. + + y : array_like + The physical x coordinates of the input charge. + + trg_threshold : float + Trigger threshold in electron equivalent. + + zero_sup_threshold : int + Zero suppression threshold in ADC counts. + + offset : int + Optional offset in ADC counts to be applied before the zero suppression. + """ + signal = Counter((col, row) for col, row in zip(*self.world_to_pixel(x, y))) + columns, rows, pha = np.array([[*key, value] for key, value in signal.items()]).T + # Trigger missing here! + pha = self.digitize(pha, zero_sup_threshold, offset) + seconds, microseconds, livetime = self.latch_timestamp(timestamp) + + + + + + + + + + +# The stuff below is going away as soon as we have the new class hierarchy in place. + @dataclass class DigiEvent: @@ -110,7 +404,7 @@ def from_digi(cls, row: np.ndarray, pha: np.ndarray): pad_top, pad_right, pad_bottom, pad_left = row padding = Padding(pad_top, pad_right, pad_bottom, pad_left) roi = RegionOfInterest(min_col, max_col, min_row, max_row, padding) - return cls(trigger_id, seconds, microseconds, livetime, roi, pha) + return cls(trigger_id, seconds, microseconds, livetime, pha, roi) def __call__(self, col: int, row: int) -> int: """Retrieve the pha content of the event for a given column and row. diff --git a/tests/test_digi.py b/tests/test_digi.py index 991c428..92e3971 100644 --- a/tests/test_digi.py +++ b/tests/test_digi.py @@ -19,19 +19,64 @@ import numpy as np from loguru import logger +import pytest -from hexsample.digi import DigiEvent, HexagonalReadout +from hexsample import digi from hexsample.hexagon import HexagonalLayout from hexsample.roi import Padding, RegionOfInterest -def test_digi_event(min_col : int = 106, max_col : int = 113, min_row : int = 15, - max_row : int = 22, padding : Padding = Padding(1, 2, 3, 4)): + +def test_digi_event_base(): + """Test for the base event class. + """ + pha = np.full(3, 100.) + event = digi.DigiEventBase(0, 0, 0, 0, pha) + print(event) + print(event.timestamp()) + +def test_digi_event_sparse(): + """ + """ + pha = np.array([50., 150., 25.]) + rows = np.array([1, 2, 3]) + columns = np.array([11, 12, 12]) + event = digi.DigiEventSparse(0, 0, 0, 0, pha, rows, columns) + print(event) + print(event.timestamp()) + print(event.ascii()) + # Make sure that the check on the dimensions of the row and column arrays is + # at work + with pytest.raises(RuntimeError): + rows = np.array([1, 2, 3]) + columns = np.array([11, 12, 12, 12]) + event = digi.DigiEventSparse(0, 0, 0, 0, pha, rows, columns) + with pytest.raises(RuntimeError): + rows = np.array([1, 2, 3, 4]) + columns = np.array([11, 12, 12]) + event = digi.DigiEventSparse(0, 0, 0, 0, pha, rows, columns) + +def test_digitization_sparse(layout: HexagonalLayout = HexagonalLayout.ODD_R, + num_cols: int = 100, num_rows: int = 100, pitch: float = 0.1, enc: float = 0., + gain: float = 0.5, num_pairs: int = 1000, trg_threshold: float = 200.): + """ + """ + readout = digi.HexagonalReadoutSparse(layout, num_cols, num_rows, pitch, enc, gain) + # Pick out a particular pixel... + col, row = num_cols // 3, num_rows // 4 + logger.debug(f'Testing pixel ({col}, {row})...') + # ... create the x and y arrays of the pair positions in the center of the pixel. + x0, y0 = readout.pixel_to_world(col, row) + x, y = np.full(num_pairs, x0), np.full(num_pairs, y0) + a = readout.read(0., x, y, 100.) + +def test_digi_event(min_col: int = 106, max_col: int = 113, min_row: int = 15, + max_row: int = 22, padding: Padding = Padding(1, 2, 3, 4)): """Build a digi event and make sure it makes sense. """ roi = RegionOfInterest(min_col, max_col, min_row, max_row, padding) # The pha is basically the serial readout order, here. pha = np.arange(roi.size) - evt = DigiEvent(0, 0, 0, 0, roi, pha) + evt = digi.DigiEvent(0, 0, 0, 0, roi, pha) print(evt.ascii()) i, j = 0, 2 assert evt.pha[i, j] == 2 @@ -44,18 +89,18 @@ def test_digi_event_comparison(): padding = Padding(2) roi = RegionOfInterest(10, 23, 10, 23, padding) pha = np.full(roi.size, 2) - evt1 = DigiEvent(0, 0, 0, 0, roi, pha) - evt2 = DigiEvent(0, 0, 0, 0, roi, 1. * pha) - evt3 = DigiEvent(0, 0, 0, 0, roi, 2. * pha) + evt1 = digi.DigiEvent(0, 0, 0, 0, roi, pha) + evt2 = digi.DigiEvent(0, 0, 0, 0, roi, 1. * pha) + evt3 = digi.DigiEvent(0, 0, 0, 0, roi, 2. * pha) assert evt1 == evt2 assert evt1 != evt3 -def test_digitization(layout : HexagonalLayout = HexagonalLayout.ODD_R, num_cols : int = 100, - num_rows : int = 100, pitch : float = 0.1, enc : float = 0., gain : float = 0.5, - num_pairs : int = 1000, trg_threshold : float = 200., padding : Padding = Padding(1)): +def test_digitization(layout: HexagonalLayout = HexagonalLayout.ODD_R, num_cols: int = 100, + num_rows: int = 100, pitch: float = 0.1, enc: float = 0., gain: float = 0.5, + num_pairs: int = 1000, trg_threshold: float = 200., padding: Padding = Padding(1)): """Create a fake digi event and test all the steps of the digitization. """ - readout = HexagonalReadout(layout, num_cols, num_rows, pitch, enc, gain) + readout = digi.HexagonalReadout(layout, num_cols, num_rows, pitch, enc, gain) # Pick out a particular pixel... col, row = num_cols // 3, num_rows // 4 logger.debug(f'Testing pixel ({col}, {row})...') From cb257323108a95f68c722ce0709ae6d2a1ff85af Mon Sep 17 00:00:00 2001 From: chiara Date: Thu, 11 Apr 2024 18:27:03 +0200 Subject: [PATCH 02/52] Continued implementation of sparse readout class --- hexsample/digi.py | 46 +++++++++++++++++++++++++++++++++++----------- tests/test_digi.py | 18 +++++++++++++----- 2 files changed, 48 insertions(+), 16 deletions(-) diff --git a/hexsample/digi.py b/hexsample/digi.py index bdf4d8a..500942c 100644 --- a/hexsample/digi.py +++ b/hexsample/digi.py @@ -212,17 +212,10 @@ def __init__(self, layout: HexagonalLayout, num_cols: int, num_rows: int, self.gain = gain self.shape = (self.num_rows, self.num_cols) self.trigger_id = -1 - - - -class HexagonalReadoutSparse(HexagonalReadoutBase): - - """ - """ - + @staticmethod def zero_suppress(array: np.ndarray, threshold: float) -> None: - """Utility function to zero-suppress an generic array. + """Utility function to zero-suppress a generic array. This is returning an array of the same shape of the input where all the values lower or equal than the zero suppression threshold are set to zero. @@ -276,7 +269,7 @@ def digitize(self, pha: np.ndarray, zero_sup_threshold: int = 0, if self.enc > 0: pha += rng.generator.normal(0., self.enc, size=pha.shape) # ... apply the conversion between electrons and ADC counts... - pha *= self.gain + pha = pha*self.gain # ... round to the neirest integer... pha = np.round(pha).astype(int) # ... if necessary, add the offset for diagnostic events... @@ -287,9 +280,39 @@ def digitize(self, pha: np.ndarray, zero_sup_threshold: int = 0, # array as the BEE would have. return pha.flatten() + + +class HexagonalReadoutSparse(HexagonalReadoutBase): + + """Description of a pixel sparse readout chip on a hexagonal matrix. + In the following readout, no ROI is formed, every (and only) triggered pixel of + the event is kept with its positional information in (col, row) format on the + hexagonal grid. + + Arguments + --------- + layout : HexagonalLayout + The layout of the hexagonal matrix. + + num_cols : int + The number of columns in the readout. + + num_rows : int + The number of rows in the readout. + + pitch : float + The readout pitch in cm. + + enc : float + The equivalent noise charge in electrons. + + gain : float + The readout gain in ADC counts per electron. + """ + def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: float, zero_sup_threshold: int = 0, offset: int = 0) -> DigiEventSparse: - """Readout an event. + """Sparse readout an event. Arguments --------- @@ -316,6 +339,7 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl # Trigger missing here! pha = self.digitize(pha, zero_sup_threshold, offset) seconds, microseconds, livetime = self.latch_timestamp(timestamp) + return DigiEventSparse(self.trigger_id, seconds, microseconds, livetime, pha, columns, rows) diff --git a/tests/test_digi.py b/tests/test_digi.py index 92e3971..bb48ab4 100644 --- a/tests/test_digi.py +++ b/tests/test_digi.py @@ -35,7 +35,7 @@ def test_digi_event_base(): print(event.timestamp()) def test_digi_event_sparse(): - """ + """Test for sparse digi event class. """ pha = np.array([50., 150., 25.]) rows = np.array([1, 2, 3]) @@ -58,16 +58,24 @@ def test_digi_event_sparse(): def test_digitization_sparse(layout: HexagonalLayout = HexagonalLayout.ODD_R, num_cols: int = 100, num_rows: int = 100, pitch: float = 0.1, enc: float = 0., gain: float = 0.5, num_pairs: int = 1000, trg_threshold: float = 200.): - """ + """Test for sparse event digitalization class. """ readout = digi.HexagonalReadoutSparse(layout, num_cols, num_rows, pitch, enc, gain) # Pick out a particular pixel... col, row = num_cols // 3, num_rows // 4 - logger.debug(f'Testing pixel ({col}, {row})...') + col1, row1= num_cols // 6, num_rows // 3 + logger.debug(f'Testing pixel ({col}, {row}) and ({col1}, {row1})...') # ... create the x and y arrays of the pair positions in the center of the pixel. x0, y0 = readout.pixel_to_world(col, row) - x, y = np.full(num_pairs, x0), np.full(num_pairs, y0) - a = readout.read(0., x, y, 100.) + x1, y1 = readout.pixel_to_world(col1, row1) + x, y = np.full(int(num_pairs), x0), np.full(int(num_pairs), y0) + print(len(y)) + x = np.append(x, np.full(int(num_pairs), x1)) + y = np.append(y, np.full(int(num_pairs), y1)) + print(x,y) + print(len(x), len(y)) + a = readout.read(0., x, y, 100.) #this is a DigiEventSparse + print(a.ascii()) def test_digi_event(min_col: int = 106, max_col: int = 113, min_row: int = 15, max_row: int = 22, padding: Padding = Padding(1, 2, 3, 4)): From 921e568a115403928a6f1baa42fe34515ee731c2 Mon Sep 17 00:00:00 2001 From: Luca Baldini Date: Fri, 12 Apr 2024 13:24:59 +0200 Subject: [PATCH 03/52] Minor. --- hexsample/digi.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/hexsample/digi.py b/hexsample/digi.py index 500942c..cd64432 100644 --- a/hexsample/digi.py +++ b/hexsample/digi.py @@ -212,7 +212,7 @@ def __init__(self, layout: HexagonalLayout, num_cols: int, num_rows: int, self.gain = gain self.shape = (self.num_rows, self.num_cols) self.trigger_id = -1 - + @staticmethod def zero_suppress(array: np.ndarray, threshold: float) -> None: """Utility function to zero-suppress a generic array. @@ -269,7 +269,7 @@ def digitize(self, pha: np.ndarray, zero_sup_threshold: int = 0, if self.enc > 0: pha += rng.generator.normal(0., self.enc, size=pha.shape) # ... apply the conversion between electrons and ADC counts... - pha = pha*self.gain + pha = pha * self.gain # ... round to the neirest integer... pha = np.round(pha).astype(int) # ... if necessary, add the offset for diagnostic events... @@ -286,7 +286,7 @@ class HexagonalReadoutSparse(HexagonalReadoutBase): """Description of a pixel sparse readout chip on a hexagonal matrix. In the following readout, no ROI is formed, every (and only) triggered pixel of - the event is kept with its positional information in (col, row) format on the + the event is kept with its positional information in (col, row) format on the hexagonal grid. Arguments @@ -340,7 +340,7 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl pha = self.digitize(pha, zero_sup_threshold, offset) seconds, microseconds, livetime = self.latch_timestamp(timestamp) return DigiEventSparse(self.trigger_id, seconds, microseconds, livetime, pha, columns, rows) - + From 76933510c0b14474e048f7609c64033a6f159e73 Mon Sep 17 00:00:00 2001 From: Luca Baldini Date: Fri, 12 Apr 2024 13:57:43 +0200 Subject: [PATCH 04/52] First complete implementation of the sparse readout. --- hexsample/digi.py | 32 +++++++++++++++++++++++++++++--- tests/test_digi.py | 29 +++++++++++++++++------------ 2 files changed, 46 insertions(+), 15 deletions(-) diff --git a/hexsample/digi.py b/hexsample/digi.py index cd64432..155faae 100644 --- a/hexsample/digi.py +++ b/hexsample/digi.py @@ -213,6 +213,19 @@ def __init__(self, layout: HexagonalLayout, num_cols: int, num_rows: int, self.shape = (self.num_rows, self.num_cols) self.trigger_id = -1 + @staticmethod + def discriminate(array: np.ndarray, threshold: float) -> np.ndarray: + """Utility function acting as a simple constant-threshold discriminator + over a generic array. This returns a boolean mask with True for all the + array elements larger than the threshold. + + This is intented to avoid possible confusion between strict and loose + comparison operators (e.g., < vs <=) when comparing the content of an array + with a threshold, and all functions downstream doing this (e.g., zero_suppress) + should use this and refrain from re-implementing their own logic. + """ + return array > threshold + @staticmethod def zero_suppress(array: np.ndarray, threshold: float) -> None: """Utility function to zero-suppress a generic array. @@ -228,7 +241,8 @@ def zero_suppress(array: np.ndarray, threshold: float) -> None: threshold : float The zero suppression threshold. """ - array[array <= threshold] = 0 + mask = np.logical_not(HexagonalReadoutBase.discriminate(array, threshold)) + array[mask] = 0 @staticmethod def latch_timestamp(timestamp: float) -> Tuple[int, int, int]: @@ -265,9 +279,17 @@ def digitize(self, pha: np.ndarray, zero_sup_threshold: int = 0, offset : int Optional offset in ADC counts to be applied before the zero suppression. """ + # Note that the array type of the input pha argument is not guaranteed, here. + # Over the course of the calculation the pha is bound to be a float (the noise + # and the gain are floating-point numbere) before it is rounded to the neirest + # integer. In order to take advantage of the automatic type casting that + # numpy implements in multiplication and addition, we use the pha = pha +/* + # over the pha +/*= form. + # See https://stackoverflow.com/questions/38673531 + # # Add the noise. if self.enc > 0: - pha += rng.generator.normal(0., self.enc, size=pha.shape) + pha = pha + rng.generator.normal(0., self.enc, size=pha.shape) # ... apply the conversion between electrons and ADC counts... pha = pha * self.gain # ... round to the neirest integer... @@ -334,9 +356,13 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl offset : int Optional offset in ADC counts to be applied before the zero suppression. """ + # Sample the input positions over the readout... signal = Counter((col, row) for col, row in zip(*self.world_to_pixel(x, y))) columns, rows, pha = np.array([[*key, value] for key, value in signal.items()]).T - # Trigger missing here! + # ...apply the trigger... + trigger_mask = self.discriminate(pha, trg_threshold) + columns, rows, pha = columns[trigger_mask], rows[trigger_mask], pha[trigger_mask] + # .. and digitize the pha values. pha = self.digitize(pha, zero_sup_threshold, offset) seconds, microseconds, livetime = self.latch_timestamp(timestamp) return DigiEventSparse(self.trigger_id, seconds, microseconds, livetime, pha, columns, rows) diff --git a/tests/test_digi.py b/tests/test_digi.py index bb48ab4..00304d4 100644 --- a/tests/test_digi.py +++ b/tests/test_digi.py @@ -62,20 +62,25 @@ def test_digitization_sparse(layout: HexagonalLayout = HexagonalLayout.ODD_R, """ readout = digi.HexagonalReadoutSparse(layout, num_cols, num_rows, pitch, enc, gain) # Pick out a particular pixel... - col, row = num_cols // 3, num_rows // 4 - col1, row1= num_cols // 6, num_rows // 3 - logger.debug(f'Testing pixel ({col}, {row}) and ({col1}, {row1})...') + col1, row1 = num_cols // 3, num_rows // 4 + col2, row2 = col1 + 8, row1 + 5 + col3, row3 = col1 + 4, row1 + 2 + logger.debug(f'Testing pixel ({col1}, {row1}) and ({col2}, {row2})...') # ... create the x and y arrays of the pair positions in the center of the pixel. - x0, y0 = readout.pixel_to_world(col, row) x1, y1 = readout.pixel_to_world(col1, row1) - x, y = np.full(int(num_pairs), x0), np.full(int(num_pairs), y0) - print(len(y)) - x = np.append(x, np.full(int(num_pairs), x1)) - y = np.append(y, np.full(int(num_pairs), y1)) - print(x,y) - print(len(x), len(y)) - a = readout.read(0., x, y, 100.) #this is a DigiEventSparse - print(a.ascii()) + x2, y2 = readout.pixel_to_world(col2, row2) + x, y = np.full(int(num_pairs), x1), np.full(int(num_pairs), y1) + x = np.append(x, np.full(num_pairs, x2)) + y = np.append(y, np.full(num_pairs, y2)) + # Add one more pixel below the trigger threshold, that we want to see disappear + # in the final event. + logger.debug(f'Adding pixel ({col3}, {row3}) below teh trigger threshold...') + x3, y3 = readout.pixel_to_world(col3, row3) + n = int(0.5 * trg_threshold) + x = np.append(x, np.full(n, x3)) + y = np.append(y, np.full(n, y3)) + event = readout.read(0., x, y, 100.) #this is a DigiEventSparse + print(event.ascii()) def test_digi_event(min_col: int = 106, max_col: int = 113, min_row: int = 15, max_row: int = 22, padding: Padding = Padding(1, 2, 3, 4)): From aca48161f2a434b9c8af5bbd02a77e6ea5cefe2d Mon Sep 17 00:00:00 2001 From: chiara Date: Fri, 12 Apr 2024 13:59:35 +0200 Subject: [PATCH 05/52] Started implementing circular readout --- hexsample/digi.py | 115 ++++++++++++++++++++++++++++++++++++++++----- hexsample/roi.py | 27 +++++++++++ tests/test_digi.py | 15 ++++++ tests/test_roi.py | 13 ++++- 4 files changed, 157 insertions(+), 13 deletions(-) diff --git a/hexsample/digi.py b/hexsample/digi.py index 500942c..8832c69 100644 --- a/hexsample/digi.py +++ b/hexsample/digi.py @@ -153,6 +153,24 @@ def ascii(self, pha_width: int = 5) -> str: text += pha text += f'\n{big_space}|\n' return text + + def highest_pixel(self, absolute: bool = True) -> Tuple[int, int]: + """Return the coordinates (col, row) of the highest pixel. + + Arguments + --------- + absolute : bool + If true, the absolute coordinates (i.e., those referring to the readout + chip) are returned; otherwise the coordinates are intended relative + to the readout window (i.e., they can be used to index the pha array). + """ + # Note col and row are swapped, here, due to how the numpy array are indexed. + # pylint: disable = unbalanced-tuple-unpacking + row, col = np.unravel_index(np.argmax(self.pha), self.pha.shape) + #if absolute: + # col += self.roi.min_col + # row += self.roi.min_row + return col, row @@ -167,15 +185,29 @@ class DigiEventRectangular(DigiEventBase): @dataclass -class DigiEventCircular(DigiEventBase): +class DigiEventCircular(DigiEventSparse): - """ - """ + """Circular digitized event. - row: int - column: int + In this particular incarnation of a digitized event the ROI is built around + a central pixel, that is the one corresponding to maximum PHA. The ROI is then + always (except in border-pixel cases) composed by 7 pixels: the central one and + its 6 neighbours. + + Arguments + --------- + column : int + The column identifier of the maximum PHA pixel in the event in pixel + coordinates. + row : int + The column identifier of the maximum PHA pixel in the event in pixel + coordinates. + """ + def center_coordinates(self): + column, row = self.highest_pixel(self) + return column, row class HexagonalReadoutBase(HexagonalGrid): @@ -249,16 +281,13 @@ def latch_timestamp(timestamp: float) -> Tuple[int, int, int]: def digitize(self, pha: np.ndarray, zero_sup_threshold: int = 0, offset: int = 0) -> np.ndarray: - """Digitize the actual signal within a given ROI. + """Digitize the actual signal. Arguments --------- - signal : array_like + pha : array_like The input array of pixel signals to be digitized. - roi : RegionOfInterest - The target ROI. - zero_sup_threshold : int Zero-suppression threshold in ADC counts. @@ -269,7 +298,7 @@ def digitize(self, pha: np.ndarray, zero_sup_threshold: int = 0, if self.enc > 0: pha += rng.generator.normal(0., self.enc, size=pha.shape) # ... apply the conversion between electrons and ADC counts... - pha = pha*self.gain + pha *= self.gain # ... round to the neirest integer... pha = np.round(pha).astype(int) # ... if necessary, add the offset for diagnostic events... @@ -340,10 +369,72 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl pha = self.digitize(pha, zero_sup_threshold, offset) seconds, microseconds, livetime = self.latch_timestamp(timestamp) return DigiEventSparse(self.trigger_id, seconds, microseconds, livetime, pha, columns, rows) - + +class HexagonalReadoutCircular(HexagonalReadoutBase): + """Description of a pixel circular readout chip on a hexagonal matrix. + In the following readout, the maximum PHA pixel is found and the ROI + formed by that pixel and its 6 adjacent neighbours. + The standard shape of columns, rows and pha array is then 7, except + for events on border, that will have len<7. + Arguments + --------- + layout : HexagonalLayout + The layout of the hexagonal matrix. + num_cols : int + The number of columns in the readout. + num_rows : int + The number of rows in the readout. + + pitch : float + The readout pitch in cm. + + enc : float + The equivalent noise charge in electrons. + + gain : float + The readout gain in ADC counts per electron. + + """ + + def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: float, + zero_sup_threshold: int = 0, offset: int = 0) -> DigiEventCircular: + """Circular readout an event. + + Arguments + --------- + timestamp : float + The event timestamp. + + x : float + The physical x coordinate of the highest pha pixel. + + y : float + The physical y coordinate of the highest pha pixel. + + trg_threshold : float + Trigger threshold in electron equivalent. + + zero_sup_threshold : int + Zero suppression threshold in ADC counts. + + offset : int + Optional offset in ADC counts to be applied before the zero suppression. + """ + signal = Counter((col, row) for col, row in zip(*self.world_to_pixel(x, y))) + coords = signal.key[np.argmax(signal.values)] + print(coords) + #column, row, pha = np.array([[*key, value] for key, value in signal.items()]).T + #Checking if there is any px outside the boundaries + #and constructing the circular ROI + #circular_px_coords = [(x,y), (x, y-1), (x, y+1), (x-1, y), (x+1, y), (x-1, y-1), (x-1, y+1)] + #columns, rows, pha = np.array([[*key, value] for key, value in signal.items()]).T + # Trigger missing here! + #pha = self.digitize(pha, zero_sup_threshold, offset) + seconds, microseconds, livetime = self.latch_timestamp(timestamp) + #return DigiEventSparse(self.trigger_id, seconds, microseconds, livetime, pha, columns, rows) diff --git a/hexsample/roi.py b/hexsample/roi.py index cd35b05..167bc29 100644 --- a/hexsample/roi.py +++ b/hexsample/roi.py @@ -206,3 +206,30 @@ def in_rot(self, col: np.array, row: np.array) -> np.array: row >= self.min_row + self.padding.top, row <= self.max_row - self.padding.bottom )) + +@dataclass +class CircularRegionOfInterest: + """Class describing a circular region of interest. + + A region of interest is the datum of the logical coorinates of its center + It contains a method that return a boolean, indicating if the center is at + chip borders (True) or not (False). Those information should be exhaustive + for finding the circular ROI pixels. + """ + + col: int + row: int + + def at_border(self, chip_size: Tuple[int, int]): + """Return True if the highest PHA pixel is on the border for a given chip_size. + Note that differently from RegionOfInterest class, here the absolute coordinates + of the center pixel are the actual members of the class, so the comparison with + the number of columns and row does not need any +/- 1. + + We should consider making the chip size a class member, because it looks + kind of awkward to pass it as an argument here. + """ + num_cols, num_rows = chip_size + return self.col == 0 or self.col == num_cols or\ + self.row == 0 or self.row == num_rows + diff --git a/tests/test_digi.py b/tests/test_digi.py index bb48ab4..fc66a6f 100644 --- a/tests/test_digi.py +++ b/tests/test_digi.py @@ -42,6 +42,7 @@ def test_digi_event_sparse(): columns = np.array([11, 12, 12]) event = digi.DigiEventSparse(0, 0, 0, 0, pha, rows, columns) print(event) + #print(event.highest_pixel()) print(event.timestamp()) print(event.ascii()) # Make sure that the check on the dimensions of the row and column arrays is @@ -77,6 +78,19 @@ def test_digitization_sparse(layout: HexagonalLayout = HexagonalLayout.ODD_R, a = readout.read(0., x, y, 100.) #this is a DigiEventSparse print(a.ascii()) +def test_digi_event_circular(): + """Test for circular digi event class. + """ + pha = np.array([50., 150., 25.]) + rows = np.array([1, 2, 3]) + columns = np.array([11, 12, 12]) + event = digi.DigiEventCircular(0, 0, 0, 0, pha, rows, columns) + print(event) + print(event.timestamp()) + print(event.ascii()) + # Make sure that the check on the dimensions of the row and column arrays is + # at work + def test_digi_event(min_col: int = 106, max_col: int = 113, min_row: int = 15, max_row: int = 22, padding: Padding = Padding(1, 2, 3, 4)): """Build a digi event and make sure it makes sense. @@ -85,6 +99,7 @@ def test_digi_event(min_col: int = 106, max_col: int = 113, min_row: int = 15, # The pha is basically the serial readout order, here. pha = np.arange(roi.size) evt = digi.DigiEvent(0, 0, 0, 0, roi, pha) + print(evt.highest_pixel()) print(evt.ascii()) i, j = 0, 2 assert evt.pha[i, j] == 2 diff --git a/tests/test_roi.py b/tests/test_roi.py index ef77209..8af3dc1 100644 --- a/tests/test_roi.py +++ b/tests/test_roi.py @@ -17,7 +17,7 @@ """ -from hexsample.roi import Padding, RegionOfInterest +from hexsample.roi import Padding, RegionOfInterest, CircularRegionOfInterest def test_padding(top : int = 2, right : int = 4, bottom : int = 3, left : int = 5) -> None: @@ -89,6 +89,17 @@ def test_roi(min_col : int = 0, max_col : int = 5, min_row : int = 25, print(roi.rot_slice()) print(roi.rot_mask()) +def test_circular_roi(): + """Unit test for CircularRegionOfInterest class. + """ + circ_roi = CircularRegionOfInterest(1, 1) + print(circ_roi) + print(circ_roi.at_border((5,5))) + circ_roi_at_border = CircularRegionOfInterest(5, 5) + print(circ_roi_at_border) + print(circ_roi_at_border.at_border((5,5))) + + def test_roi_comparison(): """Test the equality operator for ROI objects. """ From b49776a0faf45d9074246b069fc060e9e6d3168d Mon Sep 17 00:00:00 2001 From: Luca Baldini Date: Fri, 12 Apr 2024 14:04:22 +0200 Subject: [PATCH 06/52] Minor. --- hexsample/digi.py | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/hexsample/digi.py b/hexsample/digi.py index 5570d67..79bab9d 100644 --- a/hexsample/digi.py +++ b/hexsample/digi.py @@ -153,7 +153,7 @@ def ascii(self, pha_width: int = 5) -> str: text += pha text += f'\n{big_space}|\n' return text - + def highest_pixel(self, absolute: bool = True) -> Tuple[int, int]: """Return the coordinates (col, row) of the highest pixel. @@ -193,7 +193,7 @@ class DigiEventCircular(DigiEventSparse): a central pixel, that is the one corresponding to maximum PHA. The ROI is then always (except in border-pixel cases) composed by 7 pixels: the central one and its 6 neighbours. - + Arguments --------- column : int @@ -209,6 +209,8 @@ def center_coordinates(self): column, row = self.highest_pixel(self) return column, row + + class HexagonalReadoutBase(HexagonalGrid): """Description of a pixel readout chip on a hexagonal matrix. @@ -320,11 +322,7 @@ def digitize(self, pha: np.ndarray, zero_sup_threshold: int = 0, if self.enc > 0: pha = pha + rng.generator.normal(0., self.enc, size=pha.shape) # ... apply the conversion between electrons and ADC counts... -<<<<<<< HEAD - pha *= self.gain -======= pha = pha * self.gain ->>>>>>> 76933510c0b14474e048f7609c64033a6f159e73 # ... round to the neirest integer... pha = np.round(pha).astype(int) # ... if necessary, add the offset for diagnostic events... @@ -399,17 +397,16 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl pha = self.digitize(pha, zero_sup_threshold, offset) seconds, microseconds, livetime = self.latch_timestamp(timestamp) return DigiEventSparse(self.trigger_id, seconds, microseconds, livetime, pha, columns, rows) -<<<<<<< HEAD - + + + class HexagonalReadoutCircular(HexagonalReadoutBase): + """Description of a pixel circular readout chip on a hexagonal matrix. In the following readout, the maximum PHA pixel is found and the ROI formed by that pixel and its 6 adjacent neighbours. The standard shape of columns, rows and pha array is then 7, except for events on border, that will have len<7. -======= - ->>>>>>> 76933510c0b14474e048f7609c64033a6f159e73 Arguments --------- @@ -430,7 +427,6 @@ class HexagonalReadoutCircular(HexagonalReadoutBase): gain : float The readout gain in ADC counts per electron. - """ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: float, From 1cf5489dbac2123db7c372b140ffa599332587b5 Mon Sep 17 00:00:00 2001 From: Luca Baldini Date: Fri, 12 Apr 2024 14:27:45 +0200 Subject: [PATCH 07/52] Started working on the refactoring of the fileio. --- hexsample/fileio.py | 93 ++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 92 insertions(+), 1 deletion(-) diff --git a/hexsample/fileio.py b/hexsample/fileio.py index 84eb41a..bbb3a36 100644 --- a/hexsample/fileio.py +++ b/hexsample/fileio.py @@ -32,7 +32,7 @@ from hexsample import __version__, __tagdate__ from hexsample.mc import MonteCarloEvent -from hexsample.digi import DigiEvent +from hexsample.digi import DigiEventBase, DigiEvent from hexsample.recon import ReconEvent @@ -68,6 +68,97 @@ def _fill_mc_row(row: tables.tableextension.Row, event: MonteCarloEvent) -> None row.append() +class DigiDescriptionBase(tables.IsDescription): + + """ + """ + + trigger_id = tables.Int32Col(pos=0) + seconds = tables.Int32Col(pos=1) + microseconds = tables.Int32Col(pos=2) + livetime = tables.Int32Col(pos=3) + + def _fill_digi_row(row: tables.tableextension.Row, event: DigiEventBase) -> None: + """Helper function to fill an output table row, given a DigiEvent object. + + Note that this method of the base class is not calling the row.appen() hook, + which is delegated to the implementations in derived classes. + + .. note:: + This would have naturally belonged to the DigiDescription class as + a @staticmethod, but doing so is apparently breaking something into the + tables internals, and all of a sudden you get an exception due to the + fact that a staticmethod cannot be pickled. + """ + row['trigger_id'] = event.trigger_id + row['seconds'] = event.seconds + row['microseconds'] = event.microseconds + row['livetime'] = event.livetime + + + +class DigiDescriptionSparse(DigiDescriptionBase): + + """ + """ + + def _fill_digi_row(row: tables.tableextension.Row, event: DigiEventBase) -> None: + """Overloaded method. + """ + DigiDescriptionBase._fill_digi_row(row, event) + row.append() + + + +class DigiDescriptionRectangular(DigiDescriptionBase): + + """ + """ + + min_col = tables.Int16Col(pos=4) + max_col = tables.Int16Col(pos=5) + min_row = tables.Int16Col(pos=6) + max_row = tables.Int16Col(pos=7) + padding_top = tables.Int8Col(pos=8) + padding_right = tables.Int8Col(pos=9) + padding_bottom = tables.Int8Col(pos=10) + padding_left = tables.Int8Col(pos=11) + + def _fill_digi_row(row: tables.tableextension.Row, event: DigiEventBase) -> None: + """Overloaded method. + """ + DigiDescriptionBase._fill_digi_row(row, event) + row['min_col'] = event.roi.min_col + row['max_col'] = event.roi.max_col + row['min_row'] = event.roi.min_row + row['max_row'] = event.roi.max_row + row['padding_top'] = event.roi.padding.top + row['padding_right'] = event.roi.padding.right + row['padding_bottom'] = event.roi.padding.bottom + row['padding_left'] = event.roi.padding.left + row.append() + + + +class DigiDescriptionCircular(DigiDescriptionBase): + + """ + """ + + column = tables.Int16Col(pos=4) + row = tables.Int16Col(pos=5) + pha = tables.Int16Col(shape=7, pos=6) + + def _fill_digi_row(row: tables.tableextension.Row, event: DigiEventBase) -> None: + """Overloaded method. + """ + DigiDescriptionBase._fill_digi_row(row, event) + row['column'] = event.column + row['row'] = event.row + row['pha'] = event.pha + row.append() + + class DigiDescription(tables.IsDescription): From 1dc9b44bbb4e8f25cc346458ed7948cc14f22871 Mon Sep 17 00:00:00 2001 From: Luca Baldini Date: Fri, 12 Apr 2024 15:20:23 +0200 Subject: [PATCH 08/52] Fix for trigger id. --- hexsample/digi.py | 1 + 1 file changed, 1 insertion(+) diff --git a/hexsample/digi.py b/hexsample/digi.py index 79bab9d..31f29b2 100644 --- a/hexsample/digi.py +++ b/hexsample/digi.py @@ -396,6 +396,7 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl # .. and digitize the pha values. pha = self.digitize(pha, zero_sup_threshold, offset) seconds, microseconds, livetime = self.latch_timestamp(timestamp) + self.trigger_id += 1 return DigiEventSparse(self.trigger_id, seconds, microseconds, livetime, pha, columns, rows) From ca99f4bb89d14b2ce814dd1c4d6431e652e2f2a9 Mon Sep 17 00:00:00 2001 From: Luca Baldini Date: Fri, 12 Apr 2024 15:20:38 +0200 Subject: [PATCH 09/52] Initial import. --- scripts/test_bench.py | 64 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) create mode 100644 scripts/test_bench.py diff --git a/scripts/test_bench.py b/scripts/test_bench.py new file mode 100644 index 0000000..40c7bf2 --- /dev/null +++ b/scripts/test_bench.py @@ -0,0 +1,64 @@ +# Copyright (C) 2024 the hexample team. +# +# For the license terms see the file LICENSE, distributed along with this +# software. +# +# This program is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 2 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# General Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +"""Create a simple simulation with a sparse readout for execising the test bench. +""" + +from loguru import logger +from tqdm import tqdm + +from hexsample import rng, HEXSAMPLE_DATA +from hexsample.digi import HexagonalReadoutSparse +from hexsample.hexagon import HexagonalLayout +from hexsample.mc import PhotonList +from hexsample.source import LineForest, GaussianBeam, Source +from hexsample.sensor import Material, Sensor + + +output_file_path = HEXSAMPLE_DATA / 'test_bench_cu.txt' + +kwargs = dict(seed=None, srcelement='Cu', srclevel='K', srcposx=0., srcposy=0., + srcsigma=0.02, actmedium='Si', fano=0.116, thickness=0.03, transdiffsigma=40., + numevents=1000, layout='ODD_R', numcolumns=32, numrows=32, pitch=0.005, noise=0., + gain=1., trgthreshold=200., zsupthreshold=0., offset=0) + + +rng.initialize(seed=kwargs['seed']) +spectrum = LineForest(kwargs['srcelement'], kwargs['srclevel']) +beam = GaussianBeam(kwargs['srcposx'], kwargs['srcposy'], kwargs['srcsigma']) +source = Source(spectrum, beam) +material = Material(kwargs['actmedium'], kwargs['fano']) +sensor = Sensor(material, kwargs['thickness'], kwargs['transdiffsigma']) +photon_list = PhotonList(source, sensor, kwargs['numevents']) +args = HexagonalLayout(kwargs['layout']), kwargs['numcolumns'], kwargs['numrows'],\ + kwargs['pitch'], kwargs['noise'], kwargs['gain'] +readout = HexagonalReadoutSparse(*args) +logger.info(f'Readout chip: {readout}') +readout_args = kwargs['trgthreshold'], kwargs['zsupthreshold'], kwargs['offset'] +logger.info(f'Opening output file {output_file_path}...') +output_file = open(output_file_path, 'w') +logger.info('Starting the event loop...') +for mc_event in tqdm(photon_list): + x, y = mc_event.propagate(sensor.trans_diffusion_sigma) + event = readout.read(mc_event.timestamp, x, y, *readout_args) + output_file.write(f'{event.trigger_id} {event.timestamp():.9f} {len(event.pha)}\n') + for col, row, pha in zip(event.columns, event.rows, event.pha): + output_file.write(f'{col:2} {row:2} {pha}\n') +output_file.close() +logger.info('Output file closed.') From 2274f01593af2b2bb8cf27a8dadf9a1065a5b807 Mon Sep 17 00:00:00 2001 From: Luca Baldini Date: Fri, 12 Apr 2024 15:52:09 +0200 Subject: [PATCH 10/52] Old DigiEvent class removed and fixed unit tests --- hexsample/clustering.py | 6 +- hexsample/digi.py | 261 ++++++++++++++----------------------- hexsample/display.py | 4 +- hexsample/fileio.py | 16 +-- tests/test_digi.py | 13 +- tests/test_fileio.py | 6 +- tests/test_optimize_sim.py | 6 +- 7 files changed, 127 insertions(+), 185 deletions(-) diff --git a/hexsample/clustering.py b/hexsample/clustering.py index a6a7fe2..c1f154c 100644 --- a/hexsample/clustering.py +++ b/hexsample/clustering.py @@ -25,7 +25,7 @@ import numpy as np -from hexsample.digi import DigiEvent +from hexsample.digi import DigiEventRectangular from hexsample.hexagon import HexagonalGrid @@ -80,7 +80,7 @@ def zero_suppress(self, array: np.ndarray) -> np.ndarray: out[out <= self.zero_sup_threshold] = 0 return out - def run(self, event: DigiEvent) -> Cluster: + def run(self, event: DigiEventRectangular) -> Cluster: """Workhorse method to be reimplemented by derived classes. """ raise NotImplementedError @@ -104,7 +104,7 @@ class ClusteringNN(ClusteringBase): num_neighbors: int - def run(self, event: DigiEvent) -> Cluster: + def run(self, event: DigiEventRectangular) -> Cluster: """Overladed method. .. warning:: diff --git a/hexsample/digi.py b/hexsample/digi.py index 31f29b2..2fe40f8 100644 --- a/hexsample/digi.py +++ b/hexsample/digi.py @@ -80,8 +80,14 @@ def timestamp(self) -> float: """ return self.seconds + 1.e-6 * self.microseconds - def ascii(self, pha_width: int = 5): + @classmethod + def from_digi(cls, *args, **kwargs): + """Build an event object from a row in a digitized file. """ + raise NotImplementedError + + def ascii(self, *args, **kwargs): + """Ascii representation of the event. """ raise NotImplementedError @@ -154,6 +160,70 @@ def ascii(self, pha_width: int = 5) -> str: text += f'\n{big_space}|\n' return text + + +@dataclass +class DigiEventRectangular(DigiEventBase): + + """Specialized class for a digitized event based on a rectangular ROI. + + This implements the basic legacy machinery of the XPOL-I and XPOL-III readout chips. + """ + + roi: RegionOfInterest + + def __post_init__(self) -> None: + """Post-initialization code. + + Here we reshape the one-dimensional PHA array coming from the serial + readout to the proper ROI shape for all subsequent operations. + """ + try: + self.pha = self.pha.reshape(self.roi.shape()) + except ValueError as error: + logger.error(f'Error in {self.__class__.__name__} post-initializaion.') + print(self.roi) + print(f'ROI size: {self.roi.size}') + print(f'pha size: {self.pha.size}') + logger.error(error) + + def __eq__(self, other) -> bool: + """Overloaded comparison operator. + """ + return DigiEventBase.__eq__(self, other) and self.roi == other.roi + + @classmethod + def from_digi(cls, row: np.ndarray, pha: np.ndarray): + """Alternative constructor rebuilding an object from a row on a digi file. + + This is used internally when we access event data in a digi file, and + we need to reassemble a DigiEvent object from a given row of a digi table. + """ + # pylint: disable=too-many-locals + trigger_id, seconds, microseconds, livetime, min_col, max_col, min_row, max_row,\ + pad_top, pad_right, pad_bottom, pad_left = row + padding = Padding(pad_top, pad_right, pad_bottom, pad_left) + roi = RegionOfInterest(min_col, max_col, min_row, max_row, padding) + return cls(trigger_id, seconds, microseconds, livetime, pha, roi) + + def __call__(self, col: int, row: int) -> int: + """Retrieve the pha content of the event for a given column and row. + + Internally this is subtracting the proper offset to the column and row + indexes in order to convert from chip coordinates to indexes of the + underlying PHA array. Note that, due to the way arrays are indexed in numpy, + we do need to swap the column and the row. + + Arguments + --------- + col : int + The column number (in chip coordinates). + + row : int + The row number (in chip coordinates). + """ + return self.pha[row - self.roi.min_row, col - self.roi.min_col] + def highest_pixel(self, absolute: bool = True) -> Tuple[int, int]: """Return the coordinates (col, row) of the highest pixel. @@ -167,20 +237,33 @@ def highest_pixel(self, absolute: bool = True) -> Tuple[int, int]: # Note col and row are swapped, here, due to how the numpy array are indexed. # pylint: disable = unbalanced-tuple-unpacking row, col = np.unravel_index(np.argmax(self.pha), self.pha.shape) - #if absolute: - # col += self.roi.min_col - # row += self.roi.min_row + if absolute: + col += self.roi.min_col + row += self.roi.min_row return col, row - - -@dataclass -class DigiEventRectangular(DigiEventBase): - - """ - """ - - roi: RegionOfInterest + def ascii(self, pha_width: int = 5): + """Ascii representation. + """ + fmt = f'%{pha_width}d' + cols = self.roi.col_indexes() + rows = self.roi.row_indexes() + big_space = space(2 * pha_width + 1) + text = f'\n{big_space}' + text += ''.join([fmt % col for col in cols]) + text += f'\n{big_space}' + text += ''.join([fmt % (col - self.roi.min_col) for col in cols]) + text += f'\n{big_space}+{line(pha_width * self.roi.num_cols)}\n' + for row in rows: + text += f'{fmt % row} {fmt % (row - self.roi.min_row)}|' + for col in cols: + pha = fmt % self(col, row) + if not self.roi.in_rot(col, row): + pha = ansi_format(pha, AnsiFontEffect.FG_BRIGHT_BLUE) + text += pha + text += f'\n{big_space}|\n' + text += f'{self.roi}\n' + return text @@ -205,6 +288,9 @@ class DigiEventCircular(DigiEventSparse): coordinates. """ + column: int + row: int + def center_coordinates(self): column, row = self.highest_pixel(self) return column, row @@ -474,151 +560,6 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl # The stuff below is going away as soon as we have the new class hierarchy in place. -@dataclass -class DigiEvent: - - """Descriptor for a digitized event. - - A digitized event is the datum of a RegionOfInterest object, a trigger - identifier, a timestamp and a 1-dimensional array of ADC counts, in the - readout order. Note that the length of the pha array must match the size of - the ROI. - - Note that, since in the simulated digitization process we typically create - an ROI first, and only in a second moment a fully fledged event, we prefer - composition over inheritance, and deemed more convenient to have a - :class:`RegionOfInterest ` object as a class - member, rather than inherit from :class:`RegionOfInterest `. - - Arguments - --------- - trigger_id : int - The trigger identifier. - - seconds : int - The integer part of the timestamp. - - microseconds : int - The fractional part of the timestamp. - - roi : RegionOfInterest - The region of interest for the event. - - pha : np.ndarray - The pixel content of the event, in the form of a 1-dimensional array - whose length matches the size of the ROI. - """ - - trigger_id: int - seconds: int - microseconds: int - livetime: int - roi: RegionOfInterest - pha: np.ndarray - - def __post_init__(self) -> None: - """Post-initialization code. - - Here we reshape the one-dimensional PHA array coming from the serial - readout to the proper ROI shape for all subsequent operations. - """ - try: - self.pha = self.pha.reshape(self.roi.shape()) - except ValueError as error: - logger.error(f'Error in {self.__class__.__name__} post-initializaion.') - print(self.roi) - print(f'ROI size: {self.roi.size}') - print(f'pha size: {self.pha.size}') - logger.error(error) - - def __eq__(self, other) -> bool: - """Overloaded comparison operator. - """ - return (self.trigger_id, self.seconds, self.microseconds, self.livetime) == \ - (other.trigger_id, other.seconds, other.microseconds, other.livetime) and \ - self.roi == other.roi and np.allclose(self.pha, other.pha) - - @classmethod - def from_digi(cls, row: np.ndarray, pha: np.ndarray): - """Alternative constructor rebuilding an object from a row on a digi file. - - This is used internally when we access event data in a digi file, and - we need to reassemble a DigiEvent object from a given row of a digi table. - """ - # pylint: disable=too-many-locals - trigger_id, seconds, microseconds, livetime, min_col, max_col, min_row, max_row,\ - pad_top, pad_right, pad_bottom, pad_left = row - padding = Padding(pad_top, pad_right, pad_bottom, pad_left) - roi = RegionOfInterest(min_col, max_col, min_row, max_row, padding) - return cls(trigger_id, seconds, microseconds, livetime, pha, roi) - - def __call__(self, col: int, row: int) -> int: - """Retrieve the pha content of the event for a given column and row. - - Internally this is subtracting the proper offset to the column and row - indexes in order to convert from chip coordinates to indexes of the - underlying PHA array. Note that, due to the way arrays are indexed in numpy, - we do need to swap the column and the row. - - Arguments - --------- - col : int - The column number (in chip coordinates). - - row : int - The row number (in chip coordinates). - """ - return self.pha[row - self.roi.min_row, col - self.roi.min_col] - - def highest_pixel(self, absolute: bool = True) -> Tuple[int, int]: - """Return the coordinates (col, row) of the highest pixel. - - Arguments - --------- - absolute : bool - If true, the absolute coordinates (i.e., those referring to the readout - chip) are returned; otherwise the coordinates are intended relative - to the readout window (i.e., they can be used to index the pha array). - """ - # Note col and row are swapped, here, due to how the numpy array are indexed. - # pylint: disable = unbalanced-tuple-unpacking - row, col = np.unravel_index(np.argmax(self.pha), self.pha.shape) - if absolute: - col += self.roi.min_col - row += self.roi.min_row - return col, row - - def timestamp(self) -> float: - """Return the timestamp of the event, that is, the sum of the second and - microseconds parts of the DigiEvent contributions as a floating point number. - """ - return self.seconds + 1.e-6 * self.microseconds - - def ascii(self, pha_width: int = 5): - """Ascii representation. - """ - fmt = f'%{pha_width}d' - cols = self.roi.col_indexes() - rows = self.roi.row_indexes() - big_space = space(2 * pha_width + 1) - text = f'\n{big_space}' - text += ''.join([fmt % col for col in cols]) - text += f'\n{big_space}' - text += ''.join([fmt % (col - self.roi.min_col) for col in cols]) - text += f'\n{big_space}+{line(pha_width * self.roi.num_cols)}\n' - for row in rows: - text += f'{fmt % row} {fmt % (row - self.roi.min_row)}|' - for col in cols: - pha = fmt % self(col, row) - if not self.roi.in_rot(col, row): - pha = ansi_format(pha, AnsiFontEffect.FG_BRIGHT_BLUE) - text += pha - text += f'\n{big_space}|\n' - text += f'{self.roi}\n' - return text - - - class HexagonalReadout(HexagonalGrid): """Description of a pixel readout chip on a hexagonal matrix. @@ -841,7 +782,7 @@ def latch_timestamp(timestamp: float) -> Tuple[int, int, int]: return int(seconds), int(1000000 * microseconds), livetime def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: float, - padding: Padding, zero_sup_threshold: int = 0, offset: int = 0) -> DigiEvent: + padding: Padding, zero_sup_threshold: int = 0, offset: int = 0) -> DigiEventRectangular: """Readout an event. Arguments @@ -872,7 +813,7 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl roi, pha = self.trigger(signal, trg_threshold, min_col, min_row, padding) pha = self.digitize(pha, zero_sup_threshold, offset) seconds, microseconds, livetime = self.latch_timestamp(timestamp) - return DigiEvent(self.trigger_id, seconds, microseconds, livetime, roi, pha) + return DigiEventRectangular(self.trigger_id, seconds, microseconds, livetime, pha, roi) diff --git a/hexsample/display.py b/hexsample/display.py index 52a124f..1278576 100644 --- a/hexsample/display.py +++ b/hexsample/display.py @@ -27,7 +27,7 @@ from matplotlib.collections import PatchCollection import numpy as np -from hexsample.digi import DigiEvent +from hexsample.digi import DigiEventRectangular from hexsample.hexagon import HexagonalGrid from hexsample.plot import plt from hexsample.roi import RegionOfInterest @@ -172,7 +172,7 @@ def draw_roi(self, roi: RegionOfInterest, offset: Tuple[float, float] = (0., 0.) plt.text(x + dx - self._grid.pitch, y + dy, f'{row}', **fmt) return collection - def draw_digi_event(self, event: DigiEvent, offset: Tuple[float, float] = (0., 0.), + def draw_digi_event(self, event: DigiEventRectangular, offset: Tuple[float, float] = (0., 0.), indices: bool = True, padding: bool = True, zero_sup_threshold: float = 0, values: bool = True, **kwargs) -> HexagonCollection: """Draw an actual event int the parent hexagonal grid. diff --git a/hexsample/fileio.py b/hexsample/fileio.py index bbb3a36..e9fe276 100644 --- a/hexsample/fileio.py +++ b/hexsample/fileio.py @@ -32,7 +32,7 @@ from hexsample import __version__, __tagdate__ from hexsample.mc import MonteCarloEvent -from hexsample.digi import DigiEventBase, DigiEvent +from hexsample.digi import DigiEventBase, DigiEventRectangular from hexsample.recon import ReconEvent @@ -180,7 +180,7 @@ class DigiDescription(tables.IsDescription): padding_bottom = tables.Int8Col(pos=10) padding_left = tables.Int8Col(pos=11) -def _fill_digi_row(row: tables.tableextension.Row, event: DigiEvent) -> None: +def _fill_digi_row(row: tables.tableextension.Row, event: DigiEventBase) -> None: """Helper function to fill an output table row, given a DigiEvent object. .. note:: @@ -366,12 +366,12 @@ def __init__(self, file_path: str): self.mc_group = self.create_group(self.root, 'mc', 'Monte Carlo') self.mc_table = self.create_table(self.mc_group, *self.MC_TABLE_SPECS) - def add_row(self, digi_event: DigiEvent, mc_event: MonteCarloEvent) -> None: + def add_row(self, digi_event: DigiEventRectangular, mc_event: MonteCarloEvent) -> None: """Add one row to the file. Arguments --------- - digi : DigiEvent + digi : DigiEventRectangular The digitized event contribution. mc : MonteCarloEvent @@ -425,7 +425,7 @@ def add_row(self, recon_event: ReconEvent, mc_event: MonteCarloEvent) -> None: Arguments --------- - digi : DigiEvent + digi : DigiEventRectangular The digitized event contribution. mc : MonteCarloEvent @@ -508,7 +508,7 @@ def mc_column(self, name: str) -> np.ndarray: """ return self.mc_table.col(name) - def digi_event(self, row_index: int) -> DigiEvent: + def digi_event(self, row_index: int) -> DigiEventRectangular: """Random access to the DigiEvent part of the event contribution. Arguments @@ -518,7 +518,7 @@ def digi_event(self, row_index: int) -> DigiEvent: """ row = self.digi_table[row_index] pha = self.pha_array[row_index] - return DigiEvent.from_digi(row, pha) + return DigiEventRectangular.from_digi(row, pha) def mc_event(self, row_index: int) -> MonteCarloEvent: """Random access to the MonteCarloEvent part of the event contribution. @@ -537,7 +537,7 @@ def __iter__(self): self.__index = -1 return self - def __next__(self) -> DigiEvent: + def __next__(self) -> DigiEventRectangular: """Overloaded method for the implementation of the iterator protocol. """ self.__index += 1 diff --git a/tests/test_digi.py b/tests/test_digi.py index 51852eb..cc6a839 100644 --- a/tests/test_digi.py +++ b/tests/test_digi.py @@ -83,6 +83,7 @@ def test_digitization_sparse(layout: HexagonalLayout = HexagonalLayout.ODD_R, event = readout.read(0., x, y, 100.) #this is a DigiEventSparse print(event.ascii()) +@pytest.mark.skip('Under development') def test_digi_event_circular(): """Test for circular digi event class. """ @@ -96,14 +97,14 @@ def test_digi_event_circular(): # Make sure that the check on the dimensions of the row and column arrays is # at work -def test_digi_event(min_col: int = 106, max_col: int = 113, min_row: int = 15, +def test_digi_event_rectangular(min_col: int = 106, max_col: int = 113, min_row: int = 15, max_row: int = 22, padding: Padding = Padding(1, 2, 3, 4)): """Build a digi event and make sure it makes sense. """ roi = RegionOfInterest(min_col, max_col, min_row, max_row, padding) # The pha is basically the serial readout order, here. pha = np.arange(roi.size) - evt = digi.DigiEvent(0, 0, 0, 0, roi, pha) + evt = digi.DigiEventRectangular(0, 0, 0, 0, pha, roi) print(evt.highest_pixel()) print(evt.ascii()) i, j = 0, 2 @@ -111,15 +112,15 @@ def test_digi_event(min_col: int = 106, max_col: int = 113, min_row: int = 15, col, row = j + evt.roi.min_col, i + evt.roi.min_row assert evt(col, row) == 2 -def test_digi_event_comparison(): +def test_digi_event_rectangular_comparison(): """ """ padding = Padding(2) roi = RegionOfInterest(10, 23, 10, 23, padding) pha = np.full(roi.size, 2) - evt1 = digi.DigiEvent(0, 0, 0, 0, roi, pha) - evt2 = digi.DigiEvent(0, 0, 0, 0, roi, 1. * pha) - evt3 = digi.DigiEvent(0, 0, 0, 0, roi, 2. * pha) + evt1 = digi.DigiEventRectangular(0, 0, 0, 0, pha, roi) + evt2 = digi.DigiEventRectangular(0, 0, 0, 0, 1. * pha, roi) + evt3 = digi.DigiEventRectangular(0, 0, 0, 0, 2. * pha, roi) assert evt1 == evt2 assert evt1 != evt3 diff --git a/tests/test_fileio.py b/tests/test_fileio.py index ac026c4..33881f2 100644 --- a/tests/test_fileio.py +++ b/tests/test_fileio.py @@ -20,7 +20,7 @@ import numpy as np from hexsample import HEXSAMPLE_DATA -from hexsample.digi import DigiEvent +from hexsample.digi import DigiEventRectangular from hexsample.fileio import DigiInputFile, DigiOutputFile, ReconInputFile, ReconOutputFile,\ FileType, peek_file_type, open_input_file from hexsample.mc import MonteCarloEvent @@ -33,12 +33,12 @@ def _mc_event(index : int) -> MonteCarloEvent: """ return MonteCarloEvent(0.1 * index, 5.9, 0., 0., 0.02, 1000 + index) -def _digi_event(index : int) -> DigiEvent: +def _digi_event(index : int) -> DigiEventRectangular: """Create a bogus DigiEvent object with index-dependent properties. """ roi = RegionOfInterest(100, 107, 150, 155 + index * 2, Padding(2)) pha = np.full(roi.size, index) - return DigiEvent(index, index, index, 0, roi, pha) + return DigiEventRectangular(index, index, index, 0, pha, roi) def _test_write(file_path, num_events : int = 10): """Small test writing a bunch of toy event strcutures to file. diff --git a/tests/test_optimize_sim.py b/tests/test_optimize_sim.py index 8d365db..2a90739 100644 --- a/tests/test_optimize_sim.py +++ b/tests/test_optimize_sim.py @@ -29,7 +29,7 @@ from hexsample import xpol from hexsample.digi import HexagonalReadout from hexsample.hexagon import HexagonalLayout -from hexsample.fileio import DigiEvent +from hexsample.fileio import DigiEventRectangular from hexsample.mc import PhotonList from hexsample.sensor import SiliconSensor from hexsample.source import LineForest, GaussianBeam, Source @@ -175,7 +175,7 @@ def digitize(self, signal : np.ndarray, roi : RegionOfInterest, return pha.flatten() def read(self, timestamp : float, x : np.ndarray, y : np.ndarray, trg_threshold : float, - padding : Padding, zero_sup_threshold : int = 0, offset : int = 0) -> DigiEvent: + padding : Padding, zero_sup_threshold : int = 0, offset : int = 0) -> DigiEventRectangular: """Readout an event. Arguments @@ -207,7 +207,7 @@ def read(self, timestamp : float, x : np.ndarray, y : np.ndarray, trg_threshold roi = self.calculate_roi(trg, padding) pha = self.digitize(signal, roi, zero_sup_threshold, offset) seconds, microseconds, livetime = self.latch_timestamp(timestamp) - return DigiEvent(self.trigger_id, seconds, microseconds, livetime, roi, pha) + return DigiEventRectangular(self.trigger_id, seconds, microseconds, livetime, pha, roi) # XPOL-III like readouts, with the old and the new, streamlined implementetion. From b0f714c32449316b74d21367b4508c257532f3e8 Mon Sep 17 00:00:00 2001 From: Luca Baldini Date: Fri, 12 Apr 2024 15:59:17 +0200 Subject: [PATCH 11/52] Old generic classes removed. --- hexsample/bin/hxdisplay.py | 4 +- hexsample/bin/hxrecon.py | 4 +- hexsample/bin/hxsim.py | 4 +- hexsample/digi.py | 109 +------------------------------------ tests/test_digi.py | 2 +- tests/test_mc.py | 4 +- tests/test_optimize_sim.py | 6 +- 7 files changed, 15 insertions(+), 118 deletions(-) diff --git a/hexsample/bin/hxdisplay.py b/hexsample/bin/hxdisplay.py index 67892cf..980f77f 100755 --- a/hexsample/bin/hxdisplay.py +++ b/hexsample/bin/hxdisplay.py @@ -24,7 +24,7 @@ from hexsample import logger from hexsample.app import ArgumentParser -from hexsample.digi import HexagonalReadout +from hexsample.digi import HexagonalReadoutRectangular from hexsample.display import HexagonalGridDisplay from hexsample.fileio import DigiInputFile from hexsample.hexagon import HexagonalLayout @@ -48,7 +48,7 @@ def hxdisplay(**kwargs): header = input_file.header args = HexagonalLayout(header['layout']), header['numcolumns'], header['numrows'],\ header['pitch'], header['noise'], header['gain'] - readout = HexagonalReadout(*args) + readout = HexagonalReadoutRectangular(*args) logger.info(f'Readout chip: {readout}') display = HexagonalGridDisplay(readout) for event in input_file: diff --git a/hexsample/bin/hxrecon.py b/hexsample/bin/hxrecon.py index 3084ee6..88714e7 100755 --- a/hexsample/bin/hxrecon.py +++ b/hexsample/bin/hxrecon.py @@ -27,7 +27,7 @@ from hexsample import logger from hexsample.app import ArgumentParser, check_required_args from hexsample.clustering import ClusteringNN -from hexsample.digi import HexagonalReadout +from hexsample.digi import HexagonalReadoutRectangular from hexsample.fileio import DigiInputFile, ReconOutputFile from hexsample.hexagon import HexagonalLayout from hexsample.recon import ReconEvent @@ -56,7 +56,7 @@ def hxrecon(**kwargs): header = input_file.header args = HexagonalLayout(header['layout']), header['numcolumns'], header['numrows'],\ header['pitch'], header['noise'], header['gain'] - readout = HexagonalReadout(*args) + readout = HexagonalReadoutRectangular(*args) logger.info(f'Readout chip: {readout}') clustering = ClusteringNN(readout, kwargs['zsupthreshold'], kwargs['nneighbors']) suffix = kwargs['suffix'] diff --git a/hexsample/bin/hxsim.py b/hexsample/bin/hxsim.py index 2227100..2580b37 100755 --- a/hexsample/bin/hxsim.py +++ b/hexsample/bin/hxsim.py @@ -29,7 +29,7 @@ from hexsample import rng from hexsample import HEXSAMPLE_DATA from hexsample.app import ArgumentParser -from hexsample.digi import HexagonalReadout +from hexsample.digi import HexagonalReadoutRectangular from hexsample.fileio import DigiOutputFile from hexsample.hexagon import HexagonalLayout from hexsample.mc import PhotonList @@ -65,7 +65,7 @@ def hxsim(**kwargs): photon_list = PhotonList(source, sensor, kwargs['numevents']) args = HexagonalLayout(kwargs['layout']), kwargs['numcolumns'], kwargs['numrows'],\ kwargs['pitch'], kwargs['noise'], kwargs['gain'] - readout = HexagonalReadout(*args) + readout = HexagonalReadoutRectangular(*args) logger.info(f'Readout chip: {readout}') output_file_path = kwargs.get('outfile') output_file = DigiOutputFile(output_file_path) diff --git a/hexsample/digi.py b/hexsample/digi.py index 2fe40f8..623695d 100644 --- a/hexsample/digi.py +++ b/hexsample/digi.py @@ -555,47 +555,11 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl - - - -# The stuff below is going away as soon as we have the new class hierarchy in place. - -class HexagonalReadout(HexagonalGrid): +class HexagonalReadoutRectangular(HexagonalReadoutBase): """Description of a pixel readout chip on a hexagonal matrix. - - Arguments - --------- - layout : HexagonalLayout - The layout of the hexagonal matrix. - - num_cols : int - The number of columns in the readout. - - num_rows : int - The number of rows in the readout. - - pitch : float - The readout pitch in cm. - - enc : float - The equivalent noise charge in electrons. - - gain : float - The readout gain in ADC counts per electron. """ - def __init__(self, layout: HexagonalLayout, num_cols: int, num_rows: int, - pitch: float, enc: float, gain: float) -> None: - """Constructor. - """ - # pylint: disable=too-many-arguments - super().__init__(layout, num_cols, num_rows, pitch) - self.enc = enc - self.gain = gain - self.shape = (self.num_rows, self.num_cols) - self.trigger_id = -1 - @staticmethod def sum_miniclusters(array: np.ndarray) -> np.ndarray: """Sum the values in a given numpy array over its 2 x 2 trigger miniclusters. @@ -606,23 +570,6 @@ def sum_miniclusters(array: np.ndarray) -> np.ndarray: num_rows, num_cols = array.shape return array.reshape((num_rows // 2, 2, num_cols // 2, 2)).sum(-1).sum(1) - @staticmethod - def zero_suppress(array: np.ndarray, threshold: float) -> None: - """Utility function to zero-suppress an generic array. - - This is returning an array of the same shape of the input where all the - values lower or equal than the zero suppression threshold are set to zero. - - Arguments - --------- - array : array_like - The input array. - - threshold : float - The zero suppression threshold. - """ - array[array <= threshold] = 0 - @staticmethod def is_odd(value: int) -> bool: """Return whether the input integer is odd. @@ -636,7 +583,7 @@ def is_odd(value: int) -> bool: def is_even(value: int) -> bool: """Return whether the input integer is even. """ - return not HexagonalReadout.is_odd(value) + return not HexagonalReadoutRectangular.is_odd(value) def sample(self, x: np.ndarray, y: np.ndarray) -> Tuple[Tuple[int, int], np.ndarray]: """Spatially sample a pair of arrays of x and y coordinates in physical @@ -731,56 +678,6 @@ def trigger(self, signal: np.ndarray, trg_threshold, min_col: int, min_row: int, self.trigger_id += 1 return roi, pha - def digitize(self, pha: np.ndarray, zero_sup_threshold: int = 0, - offset: int = 0) -> np.ndarray: - """Digitize the actual signal within a given ROI. - - Arguments - --------- - signal : array_like - The input array of pixel signals to be digitized. - - roi : RegionOfInterest - The target ROI. - - zero_sup_threshold : int - Zero-suppression threshold in ADC counts. - - offset : int - Optional offset in ADC counts to be applied before the zero suppression. - """ - # Add the noise. - if self.enc > 0: - pha += rng.generator.normal(0., self.enc, size=pha.shape) - # ... apply the conversion between electrons and ADC counts... - pha *= self.gain - # ... round to the neirest integer... - pha = np.round(pha).astype(int) - # ... if necessary, add the offset for diagnostic events... - pha += offset - # ... zero suppress the thing... - self.zero_suppress(pha, zero_sup_threshold) - # ... flatten the array to simulate the serial readout and return the - # array as the BEE would have. - return pha.flatten() - - @staticmethod - def latch_timestamp(timestamp: float) -> Tuple[int, int, int]: - """Latch the event timestamp and return the corresponding fields of the - digi event contribution: seconds, microseconds and livetime. - - .. warning:: - The livetime calculation is not implemented, yet. - - Arguments - --------- - timestamp : float - The ground-truth event timestamp from the event generator. - """ - microseconds, seconds = np.modf(timestamp) - livetime = 0 - return int(seconds), int(1000000 * microseconds), livetime - def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: float, padding: Padding, zero_sup_threshold: int = 0, offset: int = 0) -> DigiEventRectangular: """Readout an event. @@ -817,7 +714,7 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl -class Xpol3(HexagonalReadout): +class Xpol3(HexagonalReadoutRectangular): """Derived class representing the XPOL-III readout chip. """ diff --git a/tests/test_digi.py b/tests/test_digi.py index cc6a839..ec39b34 100644 --- a/tests/test_digi.py +++ b/tests/test_digi.py @@ -129,7 +129,7 @@ def test_digitization(layout: HexagonalLayout = HexagonalLayout.ODD_R, num_cols: num_pairs: int = 1000, trg_threshold: float = 200., padding: Padding = Padding(1)): """Create a fake digi event and test all the steps of the digitization. """ - readout = digi.HexagonalReadout(layout, num_cols, num_rows, pitch, enc, gain) + readout = digi.HexagonalReadoutRectangular(layout, num_cols, num_rows, pitch, enc, gain) # Pick out a particular pixel... col, row = num_cols // 3, num_rows // 4 logger.debug(f'Testing pixel ({col}, {row})...') diff --git a/tests/test_mc.py b/tests/test_mc.py index b0eb5e5..449afa4 100644 --- a/tests/test_mc.py +++ b/tests/test_mc.py @@ -17,7 +17,7 @@ """ -from hexsample.digi import HexagonalReadout, Padding +from hexsample.digi import HexagonalReadoutRectangular, Padding from hexsample.display import HexagonalGridDisplay from hexsample.hexagon import HexagonalGrid, HexagonalLayout from hexsample.mc import MonteCarloEvent @@ -30,7 +30,7 @@ def test_diffusion(diff_sigma=40.): grid = HexagonalGrid(HexagonalLayout.ODD_R, 2, 2, 0.005) evt = MonteCarloEvent(0., 8000., 0., 0., 0.05, 3000) x, y = evt.propagate(diff_sigma) - readout = HexagonalReadout(HexagonalLayout.ODD_R, 10, 10, 0.005, 40., 1.) + readout = HexagonalReadoutRectangular(HexagonalLayout.ODD_R, 10, 10, 0.005, 40., 1.) padding = Padding(1) digi_event = readout.read(evt.timestamp, x, y, 500., padding, 80, 0) print(digi_event.ascii()) diff --git a/tests/test_optimize_sim.py b/tests/test_optimize_sim.py index 2a90739..9530b4d 100644 --- a/tests/test_optimize_sim.py +++ b/tests/test_optimize_sim.py @@ -27,7 +27,7 @@ from hexsample import logger from hexsample import xpol -from hexsample.digi import HexagonalReadout +from hexsample.digi import HexagonalReadoutRectangular from hexsample.hexagon import HexagonalLayout from hexsample.fileio import DigiEventRectangular from hexsample.mc import PhotonList @@ -37,7 +37,7 @@ -class HexagonalReadoutCompat(HexagonalReadout): +class HexagonalReadoutCompat(HexagonalReadoutRectangular): """Compatibility class implementing the readout behavior up to hexsample 0.3.1, i.e., before the simulation optimization described in @@ -214,7 +214,7 @@ def read(self, timestamp : float, x : np.ndarray, y : np.ndarray, trg_threshold # Note that we set the noise to 0. in order to allow for a deterministic # comparison among the two readouts. OLD_READOUT = HexagonalReadoutCompat(xpol.XPOL1_LAYOUT, *xpol.XPOL3_SIZE, xpol.XPOL_PITCH, 0., 1.) -NEW_READOUT = HexagonalReadout(xpol.XPOL1_LAYOUT, *xpol.XPOL3_SIZE, xpol.XPOL_PITCH, 0., 1.) +NEW_READOUT = HexagonalReadoutRectangular(xpol.XPOL1_LAYOUT, *xpol.XPOL3_SIZE, xpol.XPOL_PITCH, 0., 1.) def _compare_readouts(x, y, trg_threshold=200., padding=Padding(2)): From 1dd58a87e27fe5ed3005d27f2e0127a790857cc2 Mon Sep 17 00:00:00 2001 From: Luca Baldini Date: Fri, 12 Apr 2024 16:13:56 +0200 Subject: [PATCH 12/52] Initial import. --- docs/readout.rst | 15 ++ hexsample/readout.py | 464 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 479 insertions(+) create mode 100644 docs/readout.rst create mode 100644 hexsample/readout.py diff --git a/docs/readout.rst b/docs/readout.rst new file mode 100644 index 0000000..08da658 --- /dev/null +++ b/docs/readout.rst @@ -0,0 +1,15 @@ +:mod:`hexsample.readout` --- Readout facilities +=============================================== + +This module contains... + + +Digitization process +-------------------- + + + +Module documentation +-------------------- + +.. automodule:: hexsample.readout diff --git a/hexsample/readout.py b/hexsample/readout.py new file mode 100644 index 0000000..a574a80 --- /dev/null +++ b/hexsample/readout.py @@ -0,0 +1,464 @@ +# Copyright (C) 2023 luca.baldini@pi.infn.it +# +# For the license terms see the file LICENSE, distributed along with this +# software. +# +# This program is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 2 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# General Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +"""Readout facilities facilities. +""" + +from collections import Counter +from dataclasses import dataclass +from typing import Tuple + +from loguru import logger +import numpy as np + +from hexsample import rng +from hexsample.digi import DigiEventSparse, DigiEventRectangular, DigiEventCircular +from hexsample.hexagon import HexagonalGrid, HexagonalLayout +from hexsample.pprint import AnsiFontEffect, ansi_format, space, line +from hexsample.roi import Padding, RegionOfInterest +from hexsample import xpol + + + +class HexagonalReadoutBase(HexagonalGrid): + + """Description of a pixel readout chip on a hexagonal matrix. + + Arguments + --------- + layout : HexagonalLayout + The layout of the hexagonal matrix. + + num_cols : int + The number of columns in the readout. + + num_rows : int + The number of rows in the readout. + + pitch : float + The readout pitch in cm. + + enc : float + The equivalent noise charge in electrons. + + gain : float + The readout gain in ADC counts per electron. + """ + + def __init__(self, layout: HexagonalLayout, num_cols: int, num_rows: int, + pitch: float, enc: float, gain: float) -> None: + """Constructor. + """ + # pylint: disable=too-many-arguments + super().__init__(layout, num_cols, num_rows, pitch) + self.enc = enc + self.gain = gain + self.shape = (self.num_rows, self.num_cols) + self.trigger_id = -1 + + @staticmethod + def discriminate(array: np.ndarray, threshold: float) -> np.ndarray: + """Utility function acting as a simple constant-threshold discriminator + over a generic array. This returns a boolean mask with True for all the + array elements larger than the threshold. + + This is intented to avoid possible confusion between strict and loose + comparison operators (e.g., < vs <=) when comparing the content of an array + with a threshold, and all functions downstream doing this (e.g., zero_suppress) + should use this and refrain from re-implementing their own logic. + """ + return array > threshold + + @staticmethod + def zero_suppress(array: np.ndarray, threshold: float) -> None: + """Utility function to zero-suppress a generic array. + + This is returning an array of the same shape of the input where all the + values lower or equal than the zero suppression threshold are set to zero. + + Arguments + --------- + array : array_like + The input array. + + threshold : float + The zero suppression threshold. + """ + mask = np.logical_not(HexagonalReadoutBase.discriminate(array, threshold)) + array[mask] = 0 + + @staticmethod + def latch_timestamp(timestamp: float) -> Tuple[int, int, int]: + """Latch the event timestamp and return the corresponding fields of the + digi event contribution: seconds, microseconds and livetime. + + .. warning:: + The livetime calculation is not implemented, yet. + + Arguments + --------- + timestamp : float + The ground-truth event timestamp from the event generator. + """ + microseconds, seconds = np.modf(timestamp) + livetime = 0 + return int(seconds), int(1000000 * microseconds), livetime + + def digitize(self, pha: np.ndarray, zero_sup_threshold: int = 0, + offset: int = 0) -> np.ndarray: + """Digitize the actual signal. + + Arguments + --------- + pha : array_like + The input array of pixel signals to be digitized. + + zero_sup_threshold : int + Zero-suppression threshold in ADC counts. + + offset : int + Optional offset in ADC counts to be applied before the zero suppression. + """ + # Note that the array type of the input pha argument is not guaranteed, here. + # Over the course of the calculation the pha is bound to be a float (the noise + # and the gain are floating-point numbere) before it is rounded to the neirest + # integer. In order to take advantage of the automatic type casting that + # numpy implements in multiplication and addition, we use the pha = pha +/* + # over the pha +/*= form. + # See https://stackoverflow.com/questions/38673531 + # + # Add the noise. + if self.enc > 0: + pha = pha + rng.generator.normal(0., self.enc, size=pha.shape) + # ... apply the conversion between electrons and ADC counts... + pha = pha * self.gain + # ... round to the neirest integer... + pha = np.round(pha).astype(int) + # ... if necessary, add the offset for diagnostic events... + pha += offset + # ... zero suppress the thing... + self.zero_suppress(pha, zero_sup_threshold) + # ... flatten the array to simulate the serial readout and return the + # array as the BEE would have. + return pha.flatten() + + + +class HexagonalReadoutSparse(HexagonalReadoutBase): + + """Description of a pixel sparse readout chip on a hexagonal matrix. + In the following readout, no ROI is formed, every (and only) triggered pixel of + the event is kept with its positional information in (col, row) format on the + hexagonal grid. + + Arguments + --------- + layout : HexagonalLayout + The layout of the hexagonal matrix. + + num_cols : int + The number of columns in the readout. + + num_rows : int + The number of rows in the readout. + + pitch : float + The readout pitch in cm. + + enc : float + The equivalent noise charge in electrons. + + gain : float + The readout gain in ADC counts per electron. + """ + + def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: float, + zero_sup_threshold: int = 0, offset: int = 0) -> DigiEventSparse: + """Sparse readout an event. + + Arguments + --------- + timestamp : float + The event timestamp. + + x : array_like + The physical x coordinates of the input charge. + + y : array_like + The physical x coordinates of the input charge. + + trg_threshold : float + Trigger threshold in electron equivalent. + + zero_sup_threshold : int + Zero suppression threshold in ADC counts. + + offset : int + Optional offset in ADC counts to be applied before the zero suppression. + """ + # Sample the input positions over the readout... + signal = Counter((col, row) for col, row in zip(*self.world_to_pixel(x, y))) + columns, rows, pha = np.array([[*key, value] for key, value in signal.items()]).T + # ...apply the trigger... + trigger_mask = self.discriminate(pha, trg_threshold) + columns, rows, pha = columns[trigger_mask], rows[trigger_mask], pha[trigger_mask] + # .. and digitize the pha values. + pha = self.digitize(pha, zero_sup_threshold, offset) + seconds, microseconds, livetime = self.latch_timestamp(timestamp) + self.trigger_id += 1 + return DigiEventSparse(self.trigger_id, seconds, microseconds, livetime, pha, columns, rows) + + + +class HexagonalReadoutRectangular(HexagonalReadoutBase): + + """Description of a pixel readout chip on a hexagonal matrix. + """ + + @staticmethod + def sum_miniclusters(array: np.ndarray) -> np.ndarray: + """Sum the values in a given numpy array over its 2 x 2 trigger miniclusters. + + Note that the shape of the target 2-dimensional array must be even in + both dimensions for the thing to work. + """ + num_rows, num_cols = array.shape + return array.reshape((num_rows // 2, 2, num_cols // 2, 2)).sum(-1).sum(1) + + @staticmethod + def is_odd(value: int) -> bool: + """Return whether the input integer is odd. + + See https://stackoverflow.com/questions/14651025/ for some metrics about + the speed of this particular implementation. + """ + return value & 0x1 + + @staticmethod + def is_even(value: int) -> bool: + """Return whether the input integer is even. + """ + return not HexagonalReadoutRectangular.is_odd(value) + + def sample(self, x: np.ndarray, y: np.ndarray) -> Tuple[Tuple[int, int], np.ndarray]: + """Spatially sample a pair of arrays of x and y coordinates in physical + space onto logical (hexagonal) coordinates in logical space. + + This is achieved by converting the (x, y) physical coordinates into the + corresponding (col, row) logical coordinates on the hexagonal grid, and + then filling a two-dimensional histogram in logical space. + + .. note:: + + The output two-dimensional histogram is restricted to the pixels with + a physical signal, in order to avoid having to deal with large sparse + arrays downstream. See https://github.com/lucabaldini/hexsample/issues/12 + for more details about the reasoning behind this. + + Arguments + --------- + x : array_like + The physical x coordinates to sample. + + y : array_like + The physical y coordinates to sample. + + Returns + ------- + min_col, min_row, signal : 3-element tuple (2 integers and an array) + The coordinates of the bottom-left corner of the smallest rectangle + containing all the signal, and the corresponding histogram of the + signal itself, in electron equivalent. + """ + # pylint: disable=invalid-name + col, row = self.world_to_pixel(x, y) + # Determine the corners of the relevant rectangle where the signal histogram + # should be built. Reminder: in our trigger minicluster arrangement the minimum + # column and row coordinates are always even and the maximum column and + # row coordinates are always odd. + min_col, max_col, min_row, max_row = col.min(), col.max(), row.min(), row.max() + if self.is_odd(min_col): + min_col -= 1 + if self.is_even(max_col): + max_col += 1 + if self.is_odd(min_row): + min_row -= 1 + if self.is_even(max_row): + max_row += 1 + # Streamlined version of a two-dimensional histogram. As obscure as it + # might seem, this four-liner is significantly faster than a call to + # np.histogram2d and allows for a substantial speedup in the event generation. + num_cols = max_col - min_col + 1 + num_rows = max_row - min_row + 1 + index = num_cols * (row - min_row) + (col - min_col) + signal = np.bincount(index, minlength=num_cols * num_rows).reshape((num_rows, num_cols)) + return min_col, min_row, signal + + def trigger(self, signal: np.ndarray, trg_threshold, min_col: int, min_row: int, + padding: Padding) -> Tuple[RegionOfInterest, np.ndarray]: + """Apply the trigger, calculate the region of interest, and pad the + signal array to the proper dimension. + + .. warning:: + This is still incorrect at the edges of the readout chip, as we are + not trimming the ROI (and the corresponding arrays) to the physical + dimensions of the chip. + """ + # pylint: disable=too-many-arguments, too-many-locals + # Sum the sampled signal into the 2 x 2 trigger miniclusters. + trg_signal = self.sum_miniclusters(signal) + # Zero-suppress the trigger signal below the trigger threshold. + self.zero_suppress(trg_signal, trg_threshold) + # This is tricky, and needs to be documented properly---basically we + # build arrays with all the (minicluster) columns and rows for which + # at least one minicluster is above threshold. The multiplicative factor + # of 2 serves the purpose of converting minicluster to pixel coordinates. + trg_cols = 2 * np.nonzero(trg_signal.sum(axis=0))[0] + trg_rows = 2 * np.nonzero(trg_signal.sum(axis=1))[0] + # Build the actual ROI in chip coordinates and initialize the RegionOfInterest + # object. + roi_min_col = min_col + trg_cols.min() - padding.left + roi_max_col = min_col + trg_cols.max() + 1 + padding.right + roi_min_row = min_row + trg_rows.min() - padding.top + roi_max_row = min_row + trg_rows.max() + 1 + padding.bottom + roi = RegionOfInterest(roi_min_col, roi_max_col, roi_min_row, roi_max_row, padding) + # And now the actual PHA array: we start with all zeroes... + pha = np.full(roi.shape(), 0.) + # ...and then we patch the original signal array into the proper submask. + num_rows, num_cols = signal.shape + start_row = padding.bottom - trg_rows.min() + start_col = padding.left - trg_cols.min() + pha[start_row:start_row + num_rows, start_col:start_col + num_cols] = signal + # And do not forget to increment the trigger identifier! + self.trigger_id += 1 + return roi, pha + + def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: float, + padding: Padding, zero_sup_threshold: int = 0, offset: int = 0) -> DigiEventRectangular: + """Readout an event. + + Arguments + --------- + timestamp : float + The event timestamp. + + x : array_like + The physical x coordinates of the input charge. + + y : array_like + The physical x coordinates of the input charge. + + trg_threshold : float + Trigger threshold in electron equivalent. + + padding : Padding + The padding to be applied to the ROT. + + zero_sup_threshold : int + Zero suppression threshold in ADC counts. + + offset : int + Optional offset in ADC counts to be applied before the zero suppression. + """ + # pylint: disable=invalid-name, too-many-arguments + min_col, min_row, signal = self.sample(x, y) + roi, pha = self.trigger(signal, trg_threshold, min_col, min_row, padding) + pha = self.digitize(pha, zero_sup_threshold, offset) + seconds, microseconds, livetime = self.latch_timestamp(timestamp) + return DigiEventRectangular(self.trigger_id, seconds, microseconds, livetime, pha, roi) + + + +class Xpol3(HexagonalReadoutRectangular): + + """Derived class representing the XPOL-III readout chip. + """ + + def __init__(self, enc: float = 20., gain: float = 1.) -> None: + """Constructor. + """ + super().__init__(xpol.XPOL1_LAYOUT, *xpol.XPOL3_SIZE, xpol.XPOL_PITCH, enc, gain) + + + +class HexagonalReadoutCircular(HexagonalReadoutBase): + + """Description of a pixel circular readout chip on a hexagonal matrix. + In the following readout, the maximum PHA pixel is found and the ROI + formed by that pixel and its 6 adjacent neighbours. + The standard shape of columns, rows and pha array is then 7, except + for events on border, that will have len<7. + + Arguments + --------- + layout : HexagonalLayout + The layout of the hexagonal matrix. + + num_cols : int + The number of columns in the readout. + + num_rows : int + The number of rows in the readout. + + pitch : float + The readout pitch in cm. + + enc : float + The equivalent noise charge in electrons. + + gain : float + The readout gain in ADC counts per electron. + """ + + def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: float, + zero_sup_threshold: int = 0, offset: int = 0) -> DigiEventCircular: + """Circular readout an event. + + Arguments + --------- + timestamp : float + The event timestamp. + + x : float + The physical x coordinate of the highest pha pixel. + + y : float + The physical y coordinate of the highest pha pixel. + + trg_threshold : float + Trigger threshold in electron equivalent. + + zero_sup_threshold : int + Zero suppression threshold in ADC counts. + + offset : int + Optional offset in ADC counts to be applied before the zero suppression. + """ + signal = Counter((col, row) for col, row in zip(*self.world_to_pixel(x, y))) + coords = signal.key[np.argmax(signal.values)] + print(coords) + #column, row, pha = np.array([[*key, value] for key, value in signal.items()]).T + #Checking if there is any px outside the boundaries + #and constructing the circular ROI + #circular_px_coords = [(x,y), (x, y-1), (x, y+1), (x-1, y), (x+1, y), (x-1, y-1), (x-1, y+1)] + #columns, rows, pha = np.array([[*key, value] for key, value in signal.items()]).T + # Trigger missing here! + #pha = self.digitize(pha, zero_sup_threshold, offset) + seconds, microseconds, livetime = self.latch_timestamp(timestamp) + #return DigiEventSparse(self.trigger_id, seconds, microseconds, livetime, pha, columns, rows) From 06c8030056693f6c03c93e88230cca446a56ab2c Mon Sep 17 00:00:00 2001 From: Luca Baldini Date: Fri, 12 Apr 2024 16:14:13 +0200 Subject: [PATCH 13/52] Old digi module split into digi and readout. --- docs/api.rst | 1 + docs/digi.rst | 8 +- hexsample/bin/hxrecon.py | 2 +- hexsample/bin/hxsim.py | 2 +- hexsample/digi.py | 431 +------------------------------------ scripts/test_bench.py | 2 +- tests/test_digi.py | 25 +-- tests/test_mc.py | 3 +- tests/test_optimize_sim.py | 2 +- 9 files changed, 23 insertions(+), 453 deletions(-) diff --git a/docs/api.rst b/docs/api.rst index 5e8e6f2..5e14482 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -27,6 +27,7 @@ List of modules mc modeling pprint + readout recon rng roi diff --git a/docs/digi.rst b/docs/digi.rst index 098779f..ccdb70b 100644 --- a/docs/digi.rst +++ b/docs/digi.rst @@ -1,20 +1,16 @@ :mod:`hexsample.digi` --- Digitization facilities ================================================= -This module contains all the digitization facilities, including the event structure -for a digitized events (:class:`DigiEvent `) and a generic -event readout structure (:class:`HexagonalReadout `). +This module contains... Digi event structure -------------------- .. literalinclude:: ../hexsample/digi.py - :pyobject: DigiEvent + :pyobject: DigiEventRectangular :end-before: __post_init__ -Digitization process --------------------- diff --git a/hexsample/bin/hxrecon.py b/hexsample/bin/hxrecon.py index 88714e7..67a7d90 100755 --- a/hexsample/bin/hxrecon.py +++ b/hexsample/bin/hxrecon.py @@ -27,7 +27,7 @@ from hexsample import logger from hexsample.app import ArgumentParser, check_required_args from hexsample.clustering import ClusteringNN -from hexsample.digi import HexagonalReadoutRectangular +from hexsample.readout import HexagonalReadoutRectangular from hexsample.fileio import DigiInputFile, ReconOutputFile from hexsample.hexagon import HexagonalLayout from hexsample.recon import ReconEvent diff --git a/hexsample/bin/hxsim.py b/hexsample/bin/hxsim.py index 2580b37..be9bce6 100755 --- a/hexsample/bin/hxsim.py +++ b/hexsample/bin/hxsim.py @@ -29,7 +29,7 @@ from hexsample import rng from hexsample import HEXSAMPLE_DATA from hexsample.app import ArgumentParser -from hexsample.digi import HexagonalReadoutRectangular +from hexsample.readout import HexagonalReadoutRectangular from hexsample.fileio import DigiOutputFile from hexsample.hexagon import HexagonalLayout from hexsample.mc import PhotonList diff --git a/hexsample/digi.py b/hexsample/digi.py index 623695d..267ed93 100644 --- a/hexsample/digi.py +++ b/hexsample/digi.py @@ -17,7 +17,7 @@ # with this program; if not, write to the Free Software Foundation Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. -"""Digitization facilities. +"""Digi event structures. """ from collections import Counter @@ -294,432 +294,3 @@ class DigiEventCircular(DigiEventSparse): def center_coordinates(self): column, row = self.highest_pixel(self) return column, row - - - -class HexagonalReadoutBase(HexagonalGrid): - - """Description of a pixel readout chip on a hexagonal matrix. - - Arguments - --------- - layout : HexagonalLayout - The layout of the hexagonal matrix. - - num_cols : int - The number of columns in the readout. - - num_rows : int - The number of rows in the readout. - - pitch : float - The readout pitch in cm. - - enc : float - The equivalent noise charge in electrons. - - gain : float - The readout gain in ADC counts per electron. - """ - - def __init__(self, layout: HexagonalLayout, num_cols: int, num_rows: int, - pitch: float, enc: float, gain: float) -> None: - """Constructor. - """ - # pylint: disable=too-many-arguments - super().__init__(layout, num_cols, num_rows, pitch) - self.enc = enc - self.gain = gain - self.shape = (self.num_rows, self.num_cols) - self.trigger_id = -1 - - @staticmethod - def discriminate(array: np.ndarray, threshold: float) -> np.ndarray: - """Utility function acting as a simple constant-threshold discriminator - over a generic array. This returns a boolean mask with True for all the - array elements larger than the threshold. - - This is intented to avoid possible confusion between strict and loose - comparison operators (e.g., < vs <=) when comparing the content of an array - with a threshold, and all functions downstream doing this (e.g., zero_suppress) - should use this and refrain from re-implementing their own logic. - """ - return array > threshold - - @staticmethod - def zero_suppress(array: np.ndarray, threshold: float) -> None: - """Utility function to zero-suppress a generic array. - - This is returning an array of the same shape of the input where all the - values lower or equal than the zero suppression threshold are set to zero. - - Arguments - --------- - array : array_like - The input array. - - threshold : float - The zero suppression threshold. - """ - mask = np.logical_not(HexagonalReadoutBase.discriminate(array, threshold)) - array[mask] = 0 - - @staticmethod - def latch_timestamp(timestamp: float) -> Tuple[int, int, int]: - """Latch the event timestamp and return the corresponding fields of the - digi event contribution: seconds, microseconds and livetime. - - .. warning:: - The livetime calculation is not implemented, yet. - - Arguments - --------- - timestamp : float - The ground-truth event timestamp from the event generator. - """ - microseconds, seconds = np.modf(timestamp) - livetime = 0 - return int(seconds), int(1000000 * microseconds), livetime - - def digitize(self, pha: np.ndarray, zero_sup_threshold: int = 0, - offset: int = 0) -> np.ndarray: - """Digitize the actual signal. - - Arguments - --------- - pha : array_like - The input array of pixel signals to be digitized. - - zero_sup_threshold : int - Zero-suppression threshold in ADC counts. - - offset : int - Optional offset in ADC counts to be applied before the zero suppression. - """ - # Note that the array type of the input pha argument is not guaranteed, here. - # Over the course of the calculation the pha is bound to be a float (the noise - # and the gain are floating-point numbere) before it is rounded to the neirest - # integer. In order to take advantage of the automatic type casting that - # numpy implements in multiplication and addition, we use the pha = pha +/* - # over the pha +/*= form. - # See https://stackoverflow.com/questions/38673531 - # - # Add the noise. - if self.enc > 0: - pha = pha + rng.generator.normal(0., self.enc, size=pha.shape) - # ... apply the conversion between electrons and ADC counts... - pha = pha * self.gain - # ... round to the neirest integer... - pha = np.round(pha).astype(int) - # ... if necessary, add the offset for diagnostic events... - pha += offset - # ... zero suppress the thing... - self.zero_suppress(pha, zero_sup_threshold) - # ... flatten the array to simulate the serial readout and return the - # array as the BEE would have. - return pha.flatten() - - - -class HexagonalReadoutSparse(HexagonalReadoutBase): - - """Description of a pixel sparse readout chip on a hexagonal matrix. - In the following readout, no ROI is formed, every (and only) triggered pixel of - the event is kept with its positional information in (col, row) format on the - hexagonal grid. - - Arguments - --------- - layout : HexagonalLayout - The layout of the hexagonal matrix. - - num_cols : int - The number of columns in the readout. - - num_rows : int - The number of rows in the readout. - - pitch : float - The readout pitch in cm. - - enc : float - The equivalent noise charge in electrons. - - gain : float - The readout gain in ADC counts per electron. - """ - - def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: float, - zero_sup_threshold: int = 0, offset: int = 0) -> DigiEventSparse: - """Sparse readout an event. - - Arguments - --------- - timestamp : float - The event timestamp. - - x : array_like - The physical x coordinates of the input charge. - - y : array_like - The physical x coordinates of the input charge. - - trg_threshold : float - Trigger threshold in electron equivalent. - - zero_sup_threshold : int - Zero suppression threshold in ADC counts. - - offset : int - Optional offset in ADC counts to be applied before the zero suppression. - """ - # Sample the input positions over the readout... - signal = Counter((col, row) for col, row in zip(*self.world_to_pixel(x, y))) - columns, rows, pha = np.array([[*key, value] for key, value in signal.items()]).T - # ...apply the trigger... - trigger_mask = self.discriminate(pha, trg_threshold) - columns, rows, pha = columns[trigger_mask], rows[trigger_mask], pha[trigger_mask] - # .. and digitize the pha values. - pha = self.digitize(pha, zero_sup_threshold, offset) - seconds, microseconds, livetime = self.latch_timestamp(timestamp) - self.trigger_id += 1 - return DigiEventSparse(self.trigger_id, seconds, microseconds, livetime, pha, columns, rows) - - - -class HexagonalReadoutCircular(HexagonalReadoutBase): - - """Description of a pixel circular readout chip on a hexagonal matrix. - In the following readout, the maximum PHA pixel is found and the ROI - formed by that pixel and its 6 adjacent neighbours. - The standard shape of columns, rows and pha array is then 7, except - for events on border, that will have len<7. - - Arguments - --------- - layout : HexagonalLayout - The layout of the hexagonal matrix. - - num_cols : int - The number of columns in the readout. - - num_rows : int - The number of rows in the readout. - - pitch : float - The readout pitch in cm. - - enc : float - The equivalent noise charge in electrons. - - gain : float - The readout gain in ADC counts per electron. - """ - - def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: float, - zero_sup_threshold: int = 0, offset: int = 0) -> DigiEventCircular: - """Circular readout an event. - - Arguments - --------- - timestamp : float - The event timestamp. - - x : float - The physical x coordinate of the highest pha pixel. - - y : float - The physical y coordinate of the highest pha pixel. - - trg_threshold : float - Trigger threshold in electron equivalent. - - zero_sup_threshold : int - Zero suppression threshold in ADC counts. - - offset : int - Optional offset in ADC counts to be applied before the zero suppression. - """ - signal = Counter((col, row) for col, row in zip(*self.world_to_pixel(x, y))) - coords = signal.key[np.argmax(signal.values)] - print(coords) - #column, row, pha = np.array([[*key, value] for key, value in signal.items()]).T - #Checking if there is any px outside the boundaries - #and constructing the circular ROI - #circular_px_coords = [(x,y), (x, y-1), (x, y+1), (x-1, y), (x+1, y), (x-1, y-1), (x-1, y+1)] - #columns, rows, pha = np.array([[*key, value] for key, value in signal.items()]).T - # Trigger missing here! - #pha = self.digitize(pha, zero_sup_threshold, offset) - seconds, microseconds, livetime = self.latch_timestamp(timestamp) - #return DigiEventSparse(self.trigger_id, seconds, microseconds, livetime, pha, columns, rows) - - - -class HexagonalReadoutRectangular(HexagonalReadoutBase): - - """Description of a pixel readout chip on a hexagonal matrix. - """ - - @staticmethod - def sum_miniclusters(array: np.ndarray) -> np.ndarray: - """Sum the values in a given numpy array over its 2 x 2 trigger miniclusters. - - Note that the shape of the target 2-dimensional array must be even in - both dimensions for the thing to work. - """ - num_rows, num_cols = array.shape - return array.reshape((num_rows // 2, 2, num_cols // 2, 2)).sum(-1).sum(1) - - @staticmethod - def is_odd(value: int) -> bool: - """Return whether the input integer is odd. - - See https://stackoverflow.com/questions/14651025/ for some metrics about - the speed of this particular implementation. - """ - return value & 0x1 - - @staticmethod - def is_even(value: int) -> bool: - """Return whether the input integer is even. - """ - return not HexagonalReadoutRectangular.is_odd(value) - - def sample(self, x: np.ndarray, y: np.ndarray) -> Tuple[Tuple[int, int], np.ndarray]: - """Spatially sample a pair of arrays of x and y coordinates in physical - space onto logical (hexagonal) coordinates in logical space. - - This is achieved by converting the (x, y) physical coordinates into the - corresponding (col, row) logical coordinates on the hexagonal grid, and - then filling a two-dimensional histogram in logical space. - - .. note:: - - The output two-dimensional histogram is restricted to the pixels with - a physical signal, in order to avoid having to deal with large sparse - arrays downstream. See https://github.com/lucabaldini/hexsample/issues/12 - for more details about the reasoning behind this. - - Arguments - --------- - x : array_like - The physical x coordinates to sample. - - y : array_like - The physical y coordinates to sample. - - Returns - ------- - min_col, min_row, signal : 3-element tuple (2 integers and an array) - The coordinates of the bottom-left corner of the smallest rectangle - containing all the signal, and the corresponding histogram of the - signal itself, in electron equivalent. - """ - # pylint: disable=invalid-name - col, row = self.world_to_pixel(x, y) - # Determine the corners of the relevant rectangle where the signal histogram - # should be built. Reminder: in our trigger minicluster arrangement the minimum - # column and row coordinates are always even and the maximum column and - # row coordinates are always odd. - min_col, max_col, min_row, max_row = col.min(), col.max(), row.min(), row.max() - if self.is_odd(min_col): - min_col -= 1 - if self.is_even(max_col): - max_col += 1 - if self.is_odd(min_row): - min_row -= 1 - if self.is_even(max_row): - max_row += 1 - # Streamlined version of a two-dimensional histogram. As obscure as it - # might seem, this four-liner is significantly faster than a call to - # np.histogram2d and allows for a substantial speedup in the event generation. - num_cols = max_col - min_col + 1 - num_rows = max_row - min_row + 1 - index = num_cols * (row - min_row) + (col - min_col) - signal = np.bincount(index, minlength=num_cols * num_rows).reshape((num_rows, num_cols)) - return min_col, min_row, signal - - def trigger(self, signal: np.ndarray, trg_threshold, min_col: int, min_row: int, - padding: Padding) -> Tuple[RegionOfInterest, np.ndarray]: - """Apply the trigger, calculate the region of interest, and pad the - signal array to the proper dimension. - - .. warning:: - This is still incorrect at the edges of the readout chip, as we are - not trimming the ROI (and the corresponding arrays) to the physical - dimensions of the chip. - """ - # pylint: disable=too-many-arguments, too-many-locals - # Sum the sampled signal into the 2 x 2 trigger miniclusters. - trg_signal = self.sum_miniclusters(signal) - # Zero-suppress the trigger signal below the trigger threshold. - self.zero_suppress(trg_signal, trg_threshold) - # This is tricky, and needs to be documented properly---basically we - # build arrays with all the (minicluster) columns and rows for which - # at least one minicluster is above threshold. The multiplicative factor - # of 2 serves the purpose of converting minicluster to pixel coordinates. - trg_cols = 2 * np.nonzero(trg_signal.sum(axis=0))[0] - trg_rows = 2 * np.nonzero(trg_signal.sum(axis=1))[0] - # Build the actual ROI in chip coordinates and initialize the RegionOfInterest - # object. - roi_min_col = min_col + trg_cols.min() - padding.left - roi_max_col = min_col + trg_cols.max() + 1 + padding.right - roi_min_row = min_row + trg_rows.min() - padding.top - roi_max_row = min_row + trg_rows.max() + 1 + padding.bottom - roi = RegionOfInterest(roi_min_col, roi_max_col, roi_min_row, roi_max_row, padding) - # And now the actual PHA array: we start with all zeroes... - pha = np.full(roi.shape(), 0.) - # ...and then we patch the original signal array into the proper submask. - num_rows, num_cols = signal.shape - start_row = padding.bottom - trg_rows.min() - start_col = padding.left - trg_cols.min() - pha[start_row:start_row + num_rows, start_col:start_col + num_cols] = signal - # And do not forget to increment the trigger identifier! - self.trigger_id += 1 - return roi, pha - - def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: float, - padding: Padding, zero_sup_threshold: int = 0, offset: int = 0) -> DigiEventRectangular: - """Readout an event. - - Arguments - --------- - timestamp : float - The event timestamp. - - x : array_like - The physical x coordinates of the input charge. - - y : array_like - The physical x coordinates of the input charge. - - trg_threshold : float - Trigger threshold in electron equivalent. - - padding : Padding - The padding to be applied to the ROT. - - zero_sup_threshold : int - Zero suppression threshold in ADC counts. - - offset : int - Optional offset in ADC counts to be applied before the zero suppression. - """ - # pylint: disable=invalid-name, too-many-arguments - min_col, min_row, signal = self.sample(x, y) - roi, pha = self.trigger(signal, trg_threshold, min_col, min_row, padding) - pha = self.digitize(pha, zero_sup_threshold, offset) - seconds, microseconds, livetime = self.latch_timestamp(timestamp) - return DigiEventRectangular(self.trigger_id, seconds, microseconds, livetime, pha, roi) - - - -class Xpol3(HexagonalReadoutRectangular): - - """Derived class representing the XPOL-III readout chip. - """ - - def __init__(self, enc: float = 20., gain: float = 1.) -> None: - """Constructor. - """ - super().__init__(xpol.XPOL1_LAYOUT, *xpol.XPOL3_SIZE, xpol.XPOL_PITCH, enc, gain) diff --git a/scripts/test_bench.py b/scripts/test_bench.py index 40c7bf2..d7e9cc3 100644 --- a/scripts/test_bench.py +++ b/scripts/test_bench.py @@ -24,7 +24,7 @@ from tqdm import tqdm from hexsample import rng, HEXSAMPLE_DATA -from hexsample.digi import HexagonalReadoutSparse +from hexsample.readout import HexagonalReadoutSparse from hexsample.hexagon import HexagonalLayout from hexsample.mc import PhotonList from hexsample.source import LineForest, GaussianBeam, Source diff --git a/tests/test_digi.py b/tests/test_digi.py index ec39b34..a7bbf1e 100644 --- a/tests/test_digi.py +++ b/tests/test_digi.py @@ -21,8 +21,9 @@ from loguru import logger import pytest -from hexsample import digi +from hexsample.digi import DigiEventBase, DigiEventSparse, DigiEventRectangular, DigiEventCircular from hexsample.hexagon import HexagonalLayout +from hexsample.readout import HexagonalReadoutRectangular, HexagonalReadoutSparse from hexsample.roi import Padding, RegionOfInterest @@ -30,7 +31,7 @@ def test_digi_event_base(): """Test for the base event class. """ pha = np.full(3, 100.) - event = digi.DigiEventBase(0, 0, 0, 0, pha) + event = DigiEventBase(0, 0, 0, 0, pha) print(event) print(event.timestamp()) @@ -40,7 +41,7 @@ def test_digi_event_sparse(): pha = np.array([50., 150., 25.]) rows = np.array([1, 2, 3]) columns = np.array([11, 12, 12]) - event = digi.DigiEventSparse(0, 0, 0, 0, pha, rows, columns) + event = DigiEventSparse(0, 0, 0, 0, pha, rows, columns) print(event) #print(event.highest_pixel()) print(event.timestamp()) @@ -50,18 +51,18 @@ def test_digi_event_sparse(): with pytest.raises(RuntimeError): rows = np.array([1, 2, 3]) columns = np.array([11, 12, 12, 12]) - event = digi.DigiEventSparse(0, 0, 0, 0, pha, rows, columns) + event = DigiEventSparse(0, 0, 0, 0, pha, rows, columns) with pytest.raises(RuntimeError): rows = np.array([1, 2, 3, 4]) columns = np.array([11, 12, 12]) - event = digi.DigiEventSparse(0, 0, 0, 0, pha, rows, columns) + event = DigiEventSparse(0, 0, 0, 0, pha, rows, columns) def test_digitization_sparse(layout: HexagonalLayout = HexagonalLayout.ODD_R, num_cols: int = 100, num_rows: int = 100, pitch: float = 0.1, enc: float = 0., gain: float = 0.5, num_pairs: int = 1000, trg_threshold: float = 200.): """Test for sparse event digitalization class. """ - readout = digi.HexagonalReadoutSparse(layout, num_cols, num_rows, pitch, enc, gain) + readout = HexagonalReadoutSparse(layout, num_cols, num_rows, pitch, enc, gain) # Pick out a particular pixel... col1, row1 = num_cols // 3, num_rows // 4 col2, row2 = col1 + 8, row1 + 5 @@ -90,7 +91,7 @@ def test_digi_event_circular(): pha = np.array([50., 150., 25.]) rows = np.array([1, 2, 3]) columns = np.array([11, 12, 12]) - event = digi.DigiEventCircular(0, 0, 0, 0, pha, rows, columns) + event = DigiEventCircular(0, 0, 0, 0, pha, rows, columns) print(event) print(event.timestamp()) print(event.ascii()) @@ -104,7 +105,7 @@ def test_digi_event_rectangular(min_col: int = 106, max_col: int = 113, min_row: roi = RegionOfInterest(min_col, max_col, min_row, max_row, padding) # The pha is basically the serial readout order, here. pha = np.arange(roi.size) - evt = digi.DigiEventRectangular(0, 0, 0, 0, pha, roi) + evt = DigiEventRectangular(0, 0, 0, 0, pha, roi) print(evt.highest_pixel()) print(evt.ascii()) i, j = 0, 2 @@ -118,9 +119,9 @@ def test_digi_event_rectangular_comparison(): padding = Padding(2) roi = RegionOfInterest(10, 23, 10, 23, padding) pha = np.full(roi.size, 2) - evt1 = digi.DigiEventRectangular(0, 0, 0, 0, pha, roi) - evt2 = digi.DigiEventRectangular(0, 0, 0, 0, 1. * pha, roi) - evt3 = digi.DigiEventRectangular(0, 0, 0, 0, 2. * pha, roi) + evt1 = DigiEventRectangular(0, 0, 0, 0, pha, roi) + evt2 = DigiEventRectangular(0, 0, 0, 0, 1. * pha, roi) + evt3 = DigiEventRectangular(0, 0, 0, 0, 2. * pha, roi) assert evt1 == evt2 assert evt1 != evt3 @@ -129,7 +130,7 @@ def test_digitization(layout: HexagonalLayout = HexagonalLayout.ODD_R, num_cols: num_pairs: int = 1000, trg_threshold: float = 200., padding: Padding = Padding(1)): """Create a fake digi event and test all the steps of the digitization. """ - readout = digi.HexagonalReadoutRectangular(layout, num_cols, num_rows, pitch, enc, gain) + readout = HexagonalReadoutRectangular(layout, num_cols, num_rows, pitch, enc, gain) # Pick out a particular pixel... col, row = num_cols // 3, num_rows // 4 logger.debug(f'Testing pixel ({col}, {row})...') diff --git a/tests/test_mc.py b/tests/test_mc.py index 449afa4..ac07ac0 100644 --- a/tests/test_mc.py +++ b/tests/test_mc.py @@ -17,11 +17,12 @@ """ -from hexsample.digi import HexagonalReadoutRectangular, Padding from hexsample.display import HexagonalGridDisplay from hexsample.hexagon import HexagonalGrid, HexagonalLayout from hexsample.mc import MonteCarloEvent from hexsample.plot import plt +from hexsample.readout import HexagonalReadoutRectangular +from hexsample.roi import Padding def test_diffusion(diff_sigma=40.): diff --git a/tests/test_optimize_sim.py b/tests/test_optimize_sim.py index 9530b4d..e212a96 100644 --- a/tests/test_optimize_sim.py +++ b/tests/test_optimize_sim.py @@ -27,7 +27,7 @@ from hexsample import logger from hexsample import xpol -from hexsample.digi import HexagonalReadoutRectangular +from hexsample.readout import HexagonalReadoutRectangular from hexsample.hexagon import HexagonalLayout from hexsample.fileio import DigiEventRectangular from hexsample.mc import PhotonList From b3b7fb38fd6a0c0b53c79db95a94abe8a41dd35a Mon Sep 17 00:00:00 2001 From: Luca Baldini Date: Fri, 12 Apr 2024 16:18:28 +0200 Subject: [PATCH 14/52] pylinted. --- hexsample/digi.py | 8 -------- hexsample/readout.py | 3 --- 2 files changed, 11 deletions(-) diff --git a/hexsample/digi.py b/hexsample/digi.py index 267ed93..cbe97f0 100644 --- a/hexsample/digi.py +++ b/hexsample/digi.py @@ -20,18 +20,14 @@ """Digi event structures. """ -from collections import Counter from dataclasses import dataclass from typing import Tuple from loguru import logger import numpy as np -from hexsample import rng -from hexsample.hexagon import HexagonalGrid, HexagonalLayout from hexsample.pprint import AnsiFontEffect, ansi_format, space, line from hexsample.roi import Padding, RegionOfInterest -from hexsample import xpol @@ -290,7 +286,3 @@ class DigiEventCircular(DigiEventSparse): column: int row: int - - def center_coordinates(self): - column, row = self.highest_pixel(self) - return column, row diff --git a/hexsample/readout.py b/hexsample/readout.py index a574a80..bb9f77a 100644 --- a/hexsample/readout.py +++ b/hexsample/readout.py @@ -21,16 +21,13 @@ """ from collections import Counter -from dataclasses import dataclass from typing import Tuple -from loguru import logger import numpy as np from hexsample import rng from hexsample.digi import DigiEventSparse, DigiEventRectangular, DigiEventCircular from hexsample.hexagon import HexagonalGrid, HexagonalLayout -from hexsample.pprint import AnsiFontEffect, ansi_format, space, line from hexsample.roi import Padding, RegionOfInterest from hexsample import xpol From 6d667ddafe20610391931e2db54672ed725f88e0 Mon Sep 17 00:00:00 2001 From: chiara Date: Mon, 15 Apr 2024 18:03:15 +0200 Subject: [PATCH 15/52] Started implementing circular readout --- hexsample/digi.py | 66 +++++++++++++++++++++++++++++++++++++++++++- hexsample/readout.py | 39 ++++++++++++++++++-------- hexsample/roi.py | 62 +++++++++++++++++++++++++++++++---------- tests/test_digi.py | 39 ++++++++++++++++++++------ 4 files changed, 171 insertions(+), 35 deletions(-) diff --git a/hexsample/digi.py b/hexsample/digi.py index cbe97f0..f2d682f 100644 --- a/hexsample/digi.py +++ b/hexsample/digi.py @@ -26,6 +26,7 @@ from loguru import logger import numpy as np +from hexsample.hexagon import HexagonalGrid, HexagonalLayout from hexsample.pprint import AnsiFontEffect, ansi_format, space, line from hexsample.roi import Padding, RegionOfInterest @@ -264,7 +265,7 @@ def ascii(self, pha_width: int = 5): @dataclass -class DigiEventCircular(DigiEventSparse): +class DigiEventCircular(DigiEventBase): """Circular digitized event. @@ -282,7 +283,70 @@ class DigiEventCircular(DigiEventSparse): row : int The column identifier of the maximum PHA pixel in the event in pixel coordinates. + + layout : HexagonalLayout + The layout of the hexagonal grid of the chip. In a circular digi event it is needed + because the map of the coordinates of the neighbors of a central pixel depend on the + layout of the hexagonal grid. """ column: int row: int + grid: HexagonalGrid + + def get_neighbors(self): + """This function returns the offset coordinates of the signal pixel, that + is the highest pixel of the event (as first element), and its neighbors + (starting from the upper left, proceiding clockwisely). + The order of the coordinates corresponds in the arrays of pha to the right + value, in order to properly reconstruct the event. + """ + # Identifying the 6 neighbours of the central pixel and saving the signal pixels + # prepending the cooridnates of the highest one... + neighbors_coords = list(self.grid.neighbors(self.column, self.row)) #returns a list of tuples + # ...and prepending the highest pha pixel to the list... + neighbors_coords.insert(0, (self.column, self.row)) + # ...dividing column and row coordinates in different arrays... + columns, rows = zip(*neighbors_coords) + columns, rows = np.array(columns), np.array(rows) + return columns, rows + + def as_dict(self) -> dict: + """Return the pixel content of the event in the form of a {(col, row): pha} + value. + + This is useful to address the pixel content at a given coordinate. We refrain, + for the moment, from implementing a __call__() hook on top of this, as the + semantics would be fundamentally different from that implemented with a + rectangular ROI, where the access happen in ROI coordinates, and not in chip + coordinates, but we might want to come back to this. + """ + columns, rows = self.get_neighbors() + return {(col, row): pha for col, row, pha in zip(columns, rows, self.pha)} + + + def ascii(self, pha_width: int = 5) -> str: + """Ascii representation. + """ + columns, rows = self.get_neighbors() + pha_dict = self.as_dict() + fmt = f'%{pha_width}d' + cols = np.arange(columns.min(), columns.max() + 1) + num_cols = cols[-1] - cols[0] + 1 + rows = np.arange(rows.min(), rows.max() + 1) + big_space = space(2 * pha_width + 1) + text = f'\n{big_space}' + text += ''.join([fmt % col for col in cols]) + text += f'\n{big_space}+{line(pha_width * num_cols)}\n' + for row in rows: + text += f' {fmt % row} |' + for col in cols: + try: + pha = fmt % pha_dict[(col, row)] + except KeyError: + pha = ' ' * pha_width + text += pha + text += f'\n{big_space}|\n' + return text + + diff --git a/hexsample/readout.py b/hexsample/readout.py index bb9f77a..34979de 100644 --- a/hexsample/readout.py +++ b/hexsample/readout.py @@ -447,15 +447,32 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl offset : int Optional offset in ADC counts to be applied before the zero suppression. """ - signal = Counter((col, row) for col, row in zip(*self.world_to_pixel(x, y))) - coords = signal.key[np.argmax(signal.values)] - print(coords) - #column, row, pha = np.array([[*key, value] for key, value in signal.items()]).T - #Checking if there is any px outside the boundaries - #and constructing the circular ROI - #circular_px_coords = [(x,y), (x, y-1), (x, y+1), (x-1, y), (x+1, y), (x-1, y-1), (x-1, y+1)] - #columns, rows, pha = np.array([[*key, value] for key, value in signal.items()]).T - # Trigger missing here! - #pha = self.digitize(pha, zero_sup_threshold, offset) + + # Sample the input positions over the readout... + sparse_signal = Counter((col, row) for col, row in zip(*self.world_to_pixel(x, y))) + # ...sampling the input position of the highest PHA pixel over the readout... + # Seeking key corresponding to a certain value: + # https://stackoverflow.com/questions/8023306/get-key-by-value-in-dictionary + col_max, row_max = list(sparse_signal.keys())[list(sparse_signal.values()).\ + index(max(sparse_signal.values()))] + # ... identifying the 6 neighbours of the central pixel and saving the signal pixels + # prepending the cooridnates of the highest one... + neighbors_coords = list(self.neighbors(col_max, row_max)) #returns a list of tuples + # ...and prepending the highest pha pixel to the list... + neighbors_coords.insert(0, (col_max, row_max)) + # ...saving pha corresponding to every neighbor (and to central pixel)... + pha = np.array([sparse_signal[px] if px in sparse_signal.keys() else 0 for px in neighbors_coords], dtype=int) + # ...apply the trigger... + # Not sure the trigger is needed, the highest px passed + # necessarily the trigger or there is no event + #trigger_mask = self.discriminate(pha, trg_threshold) + # .. and digitize the pha values. + pha = self.digitize(pha, zero_sup_threshold, offset) seconds, microseconds, livetime = self.latch_timestamp(timestamp) - #return DigiEventSparse(self.trigger_id, seconds, microseconds, livetime, pha, columns, rows) + # Do not forget to update the trigger_id! + self.trigger_id += 1 + #The correspondance between the order in the pha array and the position + #of this element is guaranteed, because the function for seeking neighbors + #is hard-coded and so returns the tuples always in the same order + return DigiEventCircular(self.trigger_id, seconds, microseconds, livetime, pha, + col_max, row_max, HexagonalGrid(self.layout, self.num_cols, self.num_rows, self.pitch)) diff --git a/hexsample/roi.py b/hexsample/roi.py index 167bc29..cc8e668 100644 --- a/hexsample/roi.py +++ b/hexsample/roi.py @@ -211,25 +211,57 @@ def in_rot(self, col: np.array, row: np.array) -> np.array: class CircularRegionOfInterest: """Class describing a circular region of interest. - A region of interest is the datum of the logical coorinates of its center - It contains a method that return a boolean, indicating if the center is at - chip borders (True) or not (False). Those information should be exhaustive - for finding the circular ROI pixels. + A region of interest is the datum of the logical coordinates of its center. + The center is, inside the whole event, the pixel with the highest ROI. + The configuration """ col: int row: int - def at_border(self, chip_size: Tuple[int, int]): - """Return True if the highest PHA pixel is on the border for a given chip_size. - Note that differently from RegionOfInterest class, here the absolute coordinates - of the center pixel are the actual members of the class, so the comparison with - the number of columns and row does not need any +/- 1. + maps = {0 : {'ul' : 5, + 'ur' : 6, + 'cl' : 4, + 'cr' : 1, + 'll' : 2, + 'lr' : 3}, + 1 : {'ul' : 6, + 'ur' : 2, + 'cl' : 0, + 'cr' : 5, + 'll' : 3, + 'lr' : 4}, + 2 : {'ul' : 4, + 'ur' : 0, + 'cl' : 6, + 'cr' : 3, + 'll' : 1, + 'lr' : 5}, + 3 : {'ul' : 0, + 'ur' : 1, + 'cl' : 2, + 'cr' : 4, + 'll' : 5, + 'lr' : 6}, + 4 : {'ul' : 1, + 'ur' : 5, + 'cl' : 3, + 'cr' : 0, + 'll' : 6, + 'lr' : 2}, + 5 : {'ul' : 2, + 'ur' : 3, + 'cl' : 1, + 'cr' : 6, + 'll' : 4, + 'lr' : 0}, + 6 : {'ul' : 3, + 'ur' : 4, + 'cl' : 5, + 'cr' : 2, + 'll' : 0, + 'lr' : 1} + } + - We should consider making the chip size a class member, because it looks - kind of awkward to pass it as an argument here. - """ - num_cols, num_rows = chip_size - return self.col == 0 or self.col == num_cols or\ - self.row == 0 or self.row == num_rows diff --git a/tests/test_digi.py b/tests/test_digi.py index a7bbf1e..8618e61 100644 --- a/tests/test_digi.py +++ b/tests/test_digi.py @@ -22,8 +22,8 @@ import pytest from hexsample.digi import DigiEventBase, DigiEventSparse, DigiEventRectangular, DigiEventCircular -from hexsample.hexagon import HexagonalLayout -from hexsample.readout import HexagonalReadoutRectangular, HexagonalReadoutSparse +from hexsample.hexagon import HexagonalGrid, HexagonalLayout +from hexsample.readout import HexagonalReadoutRectangular, HexagonalReadoutSparse, HexagonalReadoutCircular from hexsample.roi import Padding, RegionOfInterest @@ -63,7 +63,7 @@ def test_digitization_sparse(layout: HexagonalLayout = HexagonalLayout.ODD_R, """Test for sparse event digitalization class. """ readout = HexagonalReadoutSparse(layout, num_cols, num_rows, pitch, enc, gain) - # Pick out a particular pixel... + # Pick out some particular pixels... col1, row1 = num_cols // 3, num_rows // 4 col2, row2 = col1 + 8, row1 + 5 col3, row3 = col1 + 4, row1 + 2 @@ -84,20 +84,43 @@ def test_digitization_sparse(layout: HexagonalLayout = HexagonalLayout.ODD_R, event = readout.read(0., x, y, 100.) #this is a DigiEventSparse print(event.ascii()) -@pytest.mark.skip('Under development') +#@pytest.mark.skip('Under development') def test_digi_event_circular(): """Test for circular digi event class. """ - pha = np.array([50., 150., 25.]) - rows = np.array([1, 2, 3]) - columns = np.array([11, 12, 12]) - event = DigiEventCircular(0, 0, 0, 0, pha, rows, columns) + pha = np.array([550., 15., 0., 5., 72., 88, 100]) + row = 1 + column = 11 + # Defining the hexagonal grid + grid = HexagonalGrid(HexagonalLayout.ODD_R, 100, 100, 0.1) + event = DigiEventCircular(0, 0, 0, 0, pha, row, column, grid) print(event) print(event.timestamp()) print(event.ascii()) # Make sure that the check on the dimensions of the row and column arrays is # at work +#@pytest.mark.skip('Under development') +def test_digitization_circular(layout: HexagonalLayout = HexagonalLayout.ODD_R, + num_cols: int = 100, num_rows: int = 100, pitch: float = 0.1, enc: float = 0., + gain: float = 0.5, num_pairs: int = 1000, trg_threshold: float = 200.): + """Test for circular event digitalization class. + """ + readout = HexagonalReadoutCircular(layout, num_cols, num_rows, pitch, enc, gain) + # Pick out some particular pixels, we expect only the one with higher PHA + # to be saved in the DigiEventCircular. + col1, row1 = 2, 4 + col2, row2 = col1 + 8, row1 + 5 + logger.debug(f'Testing pixel ({col1}, {row1}) and ({col2}, {row2})...') + # ... create the x and y arrays of the pair positions in the center of the pixel. + x1, y1 = readout.pixel_to_world(col1, row1) + x2, y2 = readout.pixel_to_world(col2, row2) + x, y = np.full(int(num_pairs), x1), np.full(int(num_pairs), y1) + x = np.append(x, np.full(num_pairs, x2)) + y = np.append(y, np.full(num_pairs, y2)) + event = readout.read(0., x, y, 100.) #this is a DigiEventCircular + print(event.ascii()) + def test_digi_event_rectangular(min_col: int = 106, max_col: int = 113, min_row: int = 15, max_row: int = 22, padding: Padding = Padding(1, 2, 3, 4)): """Build a digi event and make sure it makes sense. From e406e34e37f92af5bef9ea5c4fbe28423c99b719 Mon Sep 17 00:00:00 2001 From: Luca Baldini Date: Tue, 16 Apr 2024 11:11:29 +0200 Subject: [PATCH 16/52] Small refactoring. --- hexsample/display.py | 4 +--- hexsample/hexagon.py | 23 +++++++++++++++++++++++ tests/test_hexagon.py | 33 +++++++++++++++++++++++++++++---- 3 files changed, 53 insertions(+), 7 deletions(-) diff --git a/hexsample/display.py b/hexsample/display.py index 1278576..5ab07a3 100644 --- a/hexsample/display.py +++ b/hexsample/display.py @@ -127,9 +127,7 @@ def draw(self, offset: Tuple[float, float] = (0., 0.), pixel_labels: bool = Fals """Draw the full grid display. """ # pylint: disable = invalid-name, too-many-locals - col = np.tile(np.arange(self._grid.num_cols), self._grid.num_rows) - row = np.repeat(np.arange(self._grid.num_rows), self._grid.num_cols) - x, y = self._grid.pixel_to_world(col, row) + col, row, x, y = self._grid.pixel_physical_coordinates() dx, dy = offset collection = HexagonCollection(x + dx, y + dy, 0.5 * self._grid.pitch, self._grid.hexagon_orientation(), **kwargs) diff --git a/hexsample/hexagon.py b/hexsample/hexagon.py index e2edaf2..3f53ee0 100644 --- a/hexsample/hexagon.py +++ b/hexsample/hexagon.py @@ -329,6 +329,29 @@ def world_to_pixel(self, x: np.array, y: np.array) -> Tuple[np.array, np.array]: col, row = self._axial_to_offset(q, r) return col, row + def pixel_logical_coordinates(self) -> Tuple[np.array, np.array]: + """Return a 2-element tuple (cols, rows) of numpy arrays containing all the + column and row indices that allow to loop over the full matrix. The specific + order in which the pixels are looped upon is completely arbitrary and, for + the sake of this function, we assume that, e.g., for a 2 x 2 grid we return + >>> cols = [0, 1, 0, 1] + >>> rows = [0, 0, 1, 1] + i.e., we loop with the following order + >>> (0, 0), (1, 0), (0, 1), (1, 1) + """ + cols = np.tile(np.arange(self.num_cols), self.num_rows) + rows = np.repeat(np.arange(self.num_rows), self.num_cols) + return cols, rows + + def pixel_physical_coordinates(self) -> Tuple[np.array, np.array, np.array, np.array]: + """Return a 4-element tuple (cols, rows, x, y) containing the logical + coordinates returned by pixel_logical_coordinates(), along with the + corresponding physical coordinates. + """ + cols, rows = self.pixel_logical_coordinates() + x, y = self.pixel_to_world(cols, rows) + return cols, rows, x, y + def __str__(self): """String formatting. """ diff --git a/tests/test_hexagon.py b/tests/test_hexagon.py index 4c122a8..5cdb651 100644 --- a/tests/test_hexagon.py +++ b/tests/test_hexagon.py @@ -23,7 +23,7 @@ from hexsample.plot import plt -def test_parity(nside : int = 10, pitch : float = 0.1): +def test_parity(nside: int = 10, pitch: float = 0.1): """Test the HexagonalGrid._parity_offset() method. """ for layout in (HexagonalLayout.EVEN_R, HexagonalLayout.EVEN_Q): @@ -35,7 +35,16 @@ def test_parity(nside : int = 10, pitch : float = 0.1): assert grid._parity_offset(0) == 0 assert grid._parity_offset(1) == -1 -def test_coordinate_transform(nside : int = 10, pitch : float = 0.1): +def test_coordinates(nside: int = 2, pitch: float = 0.1): + """Test the full list of logical and physical coordinates that we use for + looping over the matrix. + """ + for layout in HexagonalLayout: + grid = HexagonalGrid(layout, nside, nside, pitch) + print(grid.pixel_logical_coordinates()) + print(grid.pixel_physical_coordinates()) + +def test_coordinate_transform(nside: int = 10, pitch: float = 0.1): """Simple test of the coordinate transformations: we pick the four corner pixels and verify that pixel_to_world() and word_to_pixel() roundtrip. (This not really an exahustive test, and all the points are at the center @@ -48,7 +57,7 @@ def test_coordinate_transform(nside : int = 10, pitch : float = 0.1): x, y = grid.pixel_to_world(col, row) assert grid.world_to_pixel(x, y) == (col, row) -def test_display(nside : int = 10, pitch : float = 0.1): +def test_display(nside: int = 10, pitch: float = 0.1): """Display all the four possible layout in a small arrangement. """ target_col = 5 @@ -69,7 +78,7 @@ def test_display(nside : int = 10, pitch : float = 0.1): plt.scatter(x[mask], y[mask], color='b', s=4.) display.setup_gca() -def test_neighbors(nside : int = 10, pitch : float = 0.1) -> None: +def test_neighbors(nside: int = 10, pitch: float = 0.1) -> None: """ """ target_pixels = (2, 3), (7, 4) @@ -87,9 +96,25 @@ def test_neighbors(nside : int = 10, pitch : float = 0.1) -> None: plt.text(x, y, f'{i + 1}', **fmt) display.setup_gca() +def test_routing_7(nside: int = 10, pitch: float = 0.1) -> None: + """Test the routing from the pixel matrix to the ADCs on the periphery in the + 7-pixel readout strategy. + """ + fmt = dict(size='xx-small', ha='center', va='center') + for layout in HexagonalLayout: + plt.figure(f'Hexagonal routing 7 {layout}') + grid = HexagonalGrid(layout, nside, nside, pitch) + display = HexagonalGridDisplay(grid) + display.draw(pixel_labels=False) + col, row, x, y = grid.pixel_physical_coordinates() + for (_col, _row, _x, _y) in zip(col, row, x, y): + plt.text(_x, _y, f'{_col + _row}', **fmt) + display.setup_gca() + if __name__ == '__main__': test_display() test_neighbors() + test_routing_7() plt.show() From 37a95db104d00baa46daf86fe8eb316a914d4158 Mon Sep 17 00:00:00 2001 From: Luca Baldini Date: Tue, 16 Apr 2024 15:37:25 +0200 Subject: [PATCH 17/52] Minor. --- hexsample/hexagon.py | 43 +++++++++++++++++++++++++++++++++++++++++++ tests/test_hexagon.py | 5 +++-- 2 files changed, 46 insertions(+), 2 deletions(-) diff --git a/hexsample/hexagon.py b/hexsample/hexagon.py index 3f53ee0..253b79f 100644 --- a/hexsample/hexagon.py +++ b/hexsample/hexagon.py @@ -118,6 +118,48 @@ def neighbors_even_q(col: int, row: int) -> tuple: } +_ADC_SEQUENCE_SHORT = (0, 1, 5, 6, 2, 3, 4) +_ADC_SEQUENCE_SHORT_LENGTH = len(_ADC_SEQUENCE_SHORT) +_ADC_SEQUENCE_LONG = (0, 2, 5, 4, 2, 1, 4, 6, 1, 3, 6, 0, 3, 5) +_ADC_SEQUENCE_LONG_LENGTH = len(_ADC_SEQUENCE_LONG) + + +def adc_channel_odd_r(col: int, row: int) -> int: + """ + """ + return -1 + +def adc_channel_even_r(col: int, row: int) -> int: + """ + """ + start = _ADC_SEQUENCE_SHORT.index(_ADC_SEQUENCE_LONG[row % _ADC_SEQUENCE_LONG_LENGTH]) + index = (col + start) % _ADC_SEQUENCE_SHORT_LENGTH + return _ADC_SEQUENCE_SHORT[index] + +def adc_channel_odd_q(col: int, row: int) -> int: + """ + """ + return -1 + +def adc_channel_even_q(col: int, row: int) -> int: + """ + """ + start = _ADC_SEQUENCE_SHORT.index(_ADC_SEQUENCE_LONG[col % _ADC_SEQUENCE_LONG_LENGTH]) + index = (row + start) % _ADC_SEQUENCE_SHORT_LENGTH + return _ADC_SEQUENCE_SHORT[index] + + +# Lookup table for . +_ADC_PROXY_DICT = { + HexagonalLayout.ODD_R: adc_channel_odd_r, + HexagonalLayout.EVEN_R: adc_channel_even_r, + HexagonalLayout.ODD_Q: adc_channel_odd_q, + HexagonalLayout.EVEN_Q: adc_channel_even_q, +} + + + + class HexagonalGrid: # pylint: disable = too-many-instance-attributes @@ -161,6 +203,7 @@ def __init__(self, layout: HexagonalLayout, num_cols: int, num_rows: int, # Cache the proper function to retrieve the neighbor pixels from the # lookup table---this is used, e.g., for the clustering. self.neighbors = _NEIGHBORS_PROXY_DICT[self.layout] + self.adc_channel = _ADC_PROXY_DICT[self.layout] def pointy_topped(self) -> bool: """Return True if the layout is pointy-topped. diff --git a/tests/test_hexagon.py b/tests/test_hexagon.py index 5cdb651..898f824 100644 --- a/tests/test_hexagon.py +++ b/tests/test_hexagon.py @@ -18,7 +18,7 @@ import numpy as np -from hexsample.hexagon import HexagonalLayout, HexagonalGrid +from hexsample.hexagon import HexagonalLayout, HexagonalGrid, adc_channel_even_r from hexsample.display import HexagonalGridDisplay from hexsample.plot import plt @@ -108,7 +108,8 @@ def test_routing_7(nside: int = 10, pitch: float = 0.1) -> None: display.draw(pixel_labels=False) col, row, x, y = grid.pixel_physical_coordinates() for (_col, _row, _x, _y) in zip(col, row, x, y): - plt.text(_x, _y, f'{_col + _row}', **fmt) + adc = grid.adc_channel(_col, _row) + plt.text(_x, _y, f'{adc}', **fmt) display.setup_gca() From 2379b5b220b4ccb39d1b3bc2f7948fc25b3bf37a Mon Sep 17 00:00:00 2001 From: chiara Date: Tue, 16 Apr 2024 18:07:03 +0200 Subject: [PATCH 18/52] finished testing conversion from logic coords to 7-adc coords --- hexsample/hexagon.py | 69 +++++++++++++++++++++++++++++++++---------- tests/test_hexagon.py | 5 ++-- 2 files changed, 56 insertions(+), 18 deletions(-) diff --git a/hexsample/hexagon.py b/hexsample/hexagon.py index 253b79f..0c39cd2 100644 --- a/hexsample/hexagon.py +++ b/hexsample/hexagon.py @@ -118,35 +118,72 @@ def neighbors_even_q(col: int, row: int) -> tuple: } -_ADC_SEQUENCE_SHORT = (0, 1, 5, 6, 2, 3, 4) -_ADC_SEQUENCE_SHORT_LENGTH = len(_ADC_SEQUENCE_SHORT) -_ADC_SEQUENCE_LONG = (0, 2, 5, 4, 2, 1, 4, 6, 1, 3, 6, 0, 3, 5) -_ADC_SEQUENCE_LONG_LENGTH = len(_ADC_SEQUENCE_LONG) + +_N_ADC_CHANNELS = 7 +_ADC_SEQUENCE_EVEN = (0, 2, 5, 0, 3, 5, 1, 3, 6, 1, 4, 6, 2, 4) +_ADC_SEQUENCE_ODD = (0, 3, 5, 1, 3, 6, 1, 4, 6, 2, 4, 0, 2, 5) +_ADC_SEQUENCE_LENGTH = len(_ADC_SEQUENCE_EVEN) def adc_channel_odd_r(col: int, row: int) -> int: + """Transformation from offset coordinates (col, row) into 7-adc channel label, + that is an int between 0 and 6, for ODD_R grid layout. + + Arguments + --------- + col: int + column pixel logical coordinate + row: int + row pixel logical coordinate """ - """ - return -1 + start = _ADC_SEQUENCE_ODD[row % _ADC_SEQUENCE_LENGTH] + index = (col + start) % _N_ADC_CHANNELS + return index def adc_channel_even_r(col: int, row: int) -> int: + """Transformation from offset coordinates (col, row) into 7-adc channel label, + that is an int between 0 and 6, for EVEN_R grid layout. + + Arguments + --------- + col: int + column pixel logical coordinate + row: int + row pixel logical coordinate """ - """ - start = _ADC_SEQUENCE_SHORT.index(_ADC_SEQUENCE_LONG[row % _ADC_SEQUENCE_LONG_LENGTH]) - index = (col + start) % _ADC_SEQUENCE_SHORT_LENGTH - return _ADC_SEQUENCE_SHORT[index] + start = _ADC_SEQUENCE_EVEN[row % _ADC_SEQUENCE_LENGTH] + index = (col + start) % _N_ADC_CHANNELS + return index def adc_channel_odd_q(col: int, row: int) -> int: + """Transformation from offset coordinates (col, row) into 7-adc channel label, + that is an int between 0 and 6, for ODD_Q grid layout. + + Arguments + --------- + col: int + column pixel logical coordinate + row: int + row pixel logical coordinate """ - """ - return -1 + start = _ADC_SEQUENCE_ODD[col % _ADC_SEQUENCE_LENGTH] + index = (row + start) % _N_ADC_CHANNELS + return index def adc_channel_even_q(col: int, row: int) -> int: + """Transformation from offset coordinates (col, row) into 7-adc channel label, + that is an int between 0 and 6, for EVEN_Q grid layout. + + Arguments + --------- + col: int + column pixel logical coordinate + row: int + row pixel logical coordinate """ - """ - start = _ADC_SEQUENCE_SHORT.index(_ADC_SEQUENCE_LONG[col % _ADC_SEQUENCE_LONG_LENGTH]) - index = (row + start) % _ADC_SEQUENCE_SHORT_LENGTH - return _ADC_SEQUENCE_SHORT[index] + start = _ADC_SEQUENCE_EVEN[col % _ADC_SEQUENCE_LENGTH] + index = (row + start) % _N_ADC_CHANNELS + return index # Lookup table for . diff --git a/tests/test_hexagon.py b/tests/test_hexagon.py index 898f824..f85c636 100644 --- a/tests/test_hexagon.py +++ b/tests/test_hexagon.py @@ -107,6 +107,7 @@ def test_routing_7(nside: int = 10, pitch: float = 0.1) -> None: display = HexagonalGridDisplay(grid) display.draw(pixel_labels=False) col, row, x, y = grid.pixel_physical_coordinates() + print(f'layout = {grid.layout}, {grid.adc_channel(2, 1)}') for (_col, _row, _x, _y) in zip(col, row, x, y): adc = grid.adc_channel(_col, _row) plt.text(_x, _y, f'{adc}', **fmt) @@ -115,7 +116,7 @@ def test_routing_7(nside: int = 10, pitch: float = 0.1) -> None: if __name__ == '__main__': - test_display() - test_neighbors() + #test_display() + #test_neighbors() test_routing_7() plt.show() From 2f0c3a9276d7b59fd6805921ac9694fcd49dc290 Mon Sep 17 00:00:00 2001 From: chiara Date: Wed, 17 Apr 2024 16:18:55 +0200 Subject: [PATCH 19/52] Implementing circular readout --- hexsample/digi.py | 16 ++++++++-------- hexsample/readout.py | 29 ++++++++++++++--------------- tests/test_digi.py | 6 +++--- tests/test_hexagon.py | 1 - 4 files changed, 25 insertions(+), 27 deletions(-) diff --git a/hexsample/digi.py b/hexsample/digi.py index f2d682f..43a0c92 100644 --- a/hexsample/digi.py +++ b/hexsample/digi.py @@ -283,17 +283,14 @@ class DigiEventCircular(DigiEventBase): row : int The column identifier of the maximum PHA pixel in the event in pixel coordinates. - - layout : HexagonalLayout - The layout of the hexagonal grid of the chip. In a circular digi event it is needed - because the map of the coordinates of the neighbors of a central pixel depend on the - layout of the hexagonal grid. + """ column: int row: int - grid: HexagonalGrid + + ''' def get_neighbors(self): """This function returns the offset coordinates of the signal pixel, that is the highest pixel of the event (as first element), and its neighbors @@ -310,7 +307,8 @@ def get_neighbors(self): columns, rows = zip(*neighbors_coords) columns, rows = np.array(columns), np.array(rows) return columns, rows - + ''' + ''' def as_dict(self) -> dict: """Return the pixel content of the event in the form of a {(col, row): pha} value. @@ -323,8 +321,9 @@ def as_dict(self) -> dict: """ columns, rows = self.get_neighbors() return {(col, row): pha for col, row, pha in zip(columns, rows, self.pha)} + ''' - + ''' def ascii(self, pha_width: int = 5) -> str: """Ascii representation. """ @@ -348,5 +347,6 @@ def ascii(self, pha_width: int = 5) -> str: text += pha text += f'\n{big_space}|\n' return text + ''' diff --git a/hexsample/readout.py b/hexsample/readout.py index 34979de..66f74db 100644 --- a/hexsample/readout.py +++ b/hexsample/readout.py @@ -451,17 +451,19 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl # Sample the input positions over the readout... sparse_signal = Counter((col, row) for col, row in zip(*self.world_to_pixel(x, y))) # ...sampling the input position of the highest PHA pixel over the readout... - # Seeking key corresponding to a certain value: - # https://stackoverflow.com/questions/8023306/get-key-by-value-in-dictionary - col_max, row_max = list(sparse_signal.keys())[list(sparse_signal.values()).\ - index(max(sparse_signal.values()))] + # See: https://stackoverflow.com/questions/70094914/max-on-collections-counter + coord_max = max(sparse_signal, key=sparse_signal.get) + col_max, row_max = coord_max + #... and converting it in ADC channel coordinates (value from 0 to 6)... + adc_max = self.adc_channel(*coord_max) + # ... creating a 7-elements array containing the PHA of the ADC channels from 0 to 6 + # in increasing order and filling it with PHAs of the highest px and its neigbors... + pha = np.empty(7) + pha[adc_max] = sparse_signal[*coord_max] # ... identifying the 6 neighbours of the central pixel and saving the signal pixels - # prepending the cooridnates of the highest one... - neighbors_coords = list(self.neighbors(col_max, row_max)) #returns a list of tuples - # ...and prepending the highest pha pixel to the list... - neighbors_coords.insert(0, (col_max, row_max)) - # ...saving pha corresponding to every neighbor (and to central pixel)... - pha = np.array([sparse_signal[px] if px in sparse_signal.keys() else 0 for px in neighbors_coords], dtype=int) + # prepending the cooridnates of the highest one... + for coords in self.neighbors(*coord_max): + pha[self.adc_channel(*coords)] = sparse_signal[*coords] # ...apply the trigger... # Not sure the trigger is needed, the highest px passed # necessarily the trigger or there is no event @@ -471,8 +473,5 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl seconds, microseconds, livetime = self.latch_timestamp(timestamp) # Do not forget to update the trigger_id! self.trigger_id += 1 - #The correspondance between the order in the pha array and the position - #of this element is guaranteed, because the function for seeking neighbors - #is hard-coded and so returns the tuples always in the same order - return DigiEventCircular(self.trigger_id, seconds, microseconds, livetime, pha, - col_max, row_max, HexagonalGrid(self.layout, self.num_cols, self.num_rows, self.pitch)) + #The pha array is always in the order [pha(adc0), pha(adc1), pha(adc2), pha(adc3), pha(adc4), pha(adc5), pha(adc6)] + return DigiEventCircular(self.trigger_id, seconds, microseconds, livetime, pha, *coords) diff --git a/tests/test_digi.py b/tests/test_digi.py index 8618e61..0b2d42b 100644 --- a/tests/test_digi.py +++ b/tests/test_digi.py @@ -93,10 +93,10 @@ def test_digi_event_circular(): column = 11 # Defining the hexagonal grid grid = HexagonalGrid(HexagonalLayout.ODD_R, 100, 100, 0.1) - event = DigiEventCircular(0, 0, 0, 0, pha, row, column, grid) + event = DigiEventCircular(0, 0, 0, 0, pha, row, column) print(event) print(event.timestamp()) - print(event.ascii()) + #print(event.ascii()) # Make sure that the check on the dimensions of the row and column arrays is # at work @@ -119,7 +119,7 @@ def test_digitization_circular(layout: HexagonalLayout = HexagonalLayout.ODD_R, x = np.append(x, np.full(num_pairs, x2)) y = np.append(y, np.full(num_pairs, y2)) event = readout.read(0., x, y, 100.) #this is a DigiEventCircular - print(event.ascii()) + #print(event.ascii()) def test_digi_event_rectangular(min_col: int = 106, max_col: int = 113, min_row: int = 15, max_row: int = 22, padding: Padding = Padding(1, 2, 3, 4)): diff --git a/tests/test_hexagon.py b/tests/test_hexagon.py index f85c636..2fa15c8 100644 --- a/tests/test_hexagon.py +++ b/tests/test_hexagon.py @@ -107,7 +107,6 @@ def test_routing_7(nside: int = 10, pitch: float = 0.1) -> None: display = HexagonalGridDisplay(grid) display.draw(pixel_labels=False) col, row, x, y = grid.pixel_physical_coordinates() - print(f'layout = {grid.layout}, {grid.adc_channel(2, 1)}') for (_col, _row, _x, _y) in zip(col, row, x, y): adc = grid.adc_channel(_col, _row) plt.text(_x, _y, f'{adc}', **fmt) From 77b5fe401aae0f27f280f033a1c5609848a37f48 Mon Sep 17 00:00:00 2001 From: Luca Baldini Date: Fri, 19 Apr 2024 18:17:33 +0200 Subject: [PATCH 20/52] Minor. --- hexsample/fileio.py | 3 ++- hexsample/readout.py | 10 ++++++---- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/hexsample/fileio.py b/hexsample/fileio.py index e9fe276..6846800 100644 --- a/hexsample/fileio.py +++ b/hexsample/fileio.py @@ -33,6 +33,7 @@ from hexsample import __version__, __tagdate__ from hexsample.mc import MonteCarloEvent from hexsample.digi import DigiEventBase, DigiEventRectangular +from hexsample.readout import HexagonalReadoutCircular from hexsample.recon import ReconEvent @@ -147,7 +148,7 @@ class DigiDescriptionCircular(DigiDescriptionBase): column = tables.Int16Col(pos=4) row = tables.Int16Col(pos=5) - pha = tables.Int16Col(shape=7, pos=6) + pha = tables.Int16Col(shape=HexagonalReadoutCircular.NUM_PIXELS, pos=6) def _fill_digi_row(row: tables.tableextension.Row, event: DigiEventBase) -> None: """Overloaded method. diff --git a/hexsample/readout.py b/hexsample/readout.py index 66f74db..171c1f9 100644 --- a/hexsample/readout.py +++ b/hexsample/readout.py @@ -423,6 +423,8 @@ class HexagonalReadoutCircular(HexagonalReadoutBase): The readout gain in ADC counts per electron. """ + NUM_PIXELS = 7 + def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: float, zero_sup_threshold: int = 0, offset: int = 0) -> DigiEventCircular: """Circular readout an event. @@ -447,7 +449,7 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl offset : int Optional offset in ADC counts to be applied before the zero suppression. """ - + # Sample the input positions over the readout... sparse_signal = Counter((col, row) for col, row in zip(*self.world_to_pixel(x, y))) # ...sampling the input position of the highest PHA pixel over the readout... @@ -458,12 +460,12 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl adc_max = self.adc_channel(*coord_max) # ... creating a 7-elements array containing the PHA of the ADC channels from 0 to 6 # in increasing order and filling it with PHAs of the highest px and its neigbors... - pha = np.empty(7) - pha[adc_max] = sparse_signal[*coord_max] + pha = np.empty(self.NUM_PIXELS) + pha[adc_max] = sparse_signal[coord_max] # ... identifying the 6 neighbours of the central pixel and saving the signal pixels # prepending the cooridnates of the highest one... for coords in self.neighbors(*coord_max): - pha[self.adc_channel(*coords)] = sparse_signal[*coords] + pha[self.adc_channel(*coords)] = sparse_signal[coords] # ...apply the trigger... # Not sure the trigger is needed, the highest px passed # necessarily the trigger or there is no event From 28f27be0fea9bae87a4fa1db610516739a31c526 Mon Sep 17 00:00:00 2001 From: chiara Date: Fri, 19 Apr 2024 18:22:16 +0200 Subject: [PATCH 21/52] removed CircularROI class --- hexsample/roi.py | 57 ----------------------------------------------- tests/test_roi.py | 10 --------- 2 files changed, 67 deletions(-) diff --git a/hexsample/roi.py b/hexsample/roi.py index cc8e668..a7ade91 100644 --- a/hexsample/roi.py +++ b/hexsample/roi.py @@ -207,61 +207,4 @@ def in_rot(self, col: np.array, row: np.array) -> np.array: row <= self.max_row - self.padding.bottom )) -@dataclass -class CircularRegionOfInterest: - """Class describing a circular region of interest. - - A region of interest is the datum of the logical coordinates of its center. - The center is, inside the whole event, the pixel with the highest ROI. - The configuration - """ - - col: int - row: int - - maps = {0 : {'ul' : 5, - 'ur' : 6, - 'cl' : 4, - 'cr' : 1, - 'll' : 2, - 'lr' : 3}, - 1 : {'ul' : 6, - 'ur' : 2, - 'cl' : 0, - 'cr' : 5, - 'll' : 3, - 'lr' : 4}, - 2 : {'ul' : 4, - 'ur' : 0, - 'cl' : 6, - 'cr' : 3, - 'll' : 1, - 'lr' : 5}, - 3 : {'ul' : 0, - 'ur' : 1, - 'cl' : 2, - 'cr' : 4, - 'll' : 5, - 'lr' : 6}, - 4 : {'ul' : 1, - 'ur' : 5, - 'cl' : 3, - 'cr' : 0, - 'll' : 6, - 'lr' : 2}, - 5 : {'ul' : 2, - 'ur' : 3, - 'cl' : 1, - 'cr' : 6, - 'll' : 4, - 'lr' : 0}, - 6 : {'ul' : 3, - 'ur' : 4, - 'cl' : 5, - 'cr' : 2, - 'll' : 0, - 'lr' : 1} - } - - diff --git a/tests/test_roi.py b/tests/test_roi.py index 8af3dc1..551f668 100644 --- a/tests/test_roi.py +++ b/tests/test_roi.py @@ -89,16 +89,6 @@ def test_roi(min_col : int = 0, max_col : int = 5, min_row : int = 25, print(roi.rot_slice()) print(roi.rot_mask()) -def test_circular_roi(): - """Unit test for CircularRegionOfInterest class. - """ - circ_roi = CircularRegionOfInterest(1, 1) - print(circ_roi) - print(circ_roi.at_border((5,5))) - circ_roi_at_border = CircularRegionOfInterest(5, 5) - print(circ_roi_at_border) - print(circ_roi_at_border.at_border((5,5))) - def test_roi_comparison(): """Test the equality operator for ROI objects. From 398e7464c1f746699e2bd4dbcc7eec1d479eaa34 Mon Sep 17 00:00:00 2001 From: chiara Date: Fri, 19 Apr 2024 18:24:55 +0200 Subject: [PATCH 22/52] minor --- tests/test_roi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_roi.py b/tests/test_roi.py index 551f668..3e64db3 100644 --- a/tests/test_roi.py +++ b/tests/test_roi.py @@ -17,7 +17,7 @@ """ -from hexsample.roi import Padding, RegionOfInterest, CircularRegionOfInterest +from hexsample.roi import Padding, RegionOfInterest def test_padding(top : int = 2, right : int = 4, bottom : int = 3, left : int = 5) -> None: From c92f9d42023893aecd26babd2a4d4720ab64daa3 Mon Sep 17 00:00:00 2001 From: chiara Date: Sat, 20 Apr 2024 17:46:02 +0200 Subject: [PATCH 23/52] Implemented DigiEventCircular ascii function --- hexsample/digi.py | 63 ++++++++++++++------------------------------ hexsample/readout.py | 3 ++- tests/test_digi.py | 2 +- 3 files changed, 23 insertions(+), 45 deletions(-) diff --git a/hexsample/digi.py b/hexsample/digi.py index 43a0c92..ea0353e 100644 --- a/hexsample/digi.py +++ b/hexsample/digi.py @@ -26,7 +26,7 @@ from loguru import logger import numpy as np -from hexsample.hexagon import HexagonalGrid, HexagonalLayout +#from hexsample.readout import HexagonalReadoutCircular from hexsample.pprint import AnsiFontEffect, ansi_format, space, line from hexsample.roi import Padding, RegionOfInterest @@ -289,64 +289,41 @@ class DigiEventCircular(DigiEventBase): column: int row: int - - ''' - def get_neighbors(self): - """This function returns the offset coordinates of the signal pixel, that - is the highest pixel of the event (as first element), and its neighbors - (starting from the upper left, proceiding clockwisely). - The order of the coordinates corresponds in the arrays of pha to the right - value, in order to properly reconstruct the event. - """ - # Identifying the 6 neighbours of the central pixel and saving the signal pixels - # prepending the cooridnates of the highest one... - neighbors_coords = list(self.grid.neighbors(self.column, self.row)) #returns a list of tuples - # ...and prepending the highest pha pixel to the list... - neighbors_coords.insert(0, (self.column, self.row)) - # ...dividing column and row coordinates in different arrays... - columns, rows = zip(*neighbors_coords) - columns, rows = np.array(columns), np.array(rows) - return columns, rows - ''' - ''' - def as_dict(self) -> dict: - """Return the pixel content of the event in the form of a {(col, row): pha} - value. - - This is useful to address the pixel content at a given coordinate. We refrain, - for the moment, from implementing a __call__() hook on top of this, as the - semantics would be fundamentally different from that implemented with a - rectangular ROI, where the access happen in ROI coordinates, and not in chip - coordinates, but we might want to come back to this. + def __post_init__(self) -> None: + """Post-initialization code. """ - columns, rows = self.get_neighbors() - return {(col, row): pha for col, row, pha in zip(columns, rows, self.pha)} - ''' + if not len(self.pha) == 7: + raise RuntimeError(f'{self.__class__.__name__} has {len(self.pha)} PHA values' + f'instead of {7}.') + - ''' def ascii(self, pha_width: int = 5) -> str: - """Ascii representation. + """Ascii representation. + In the specific case of this class, the ascii representation is simply a px + (that is the highest PHA pixel), because the neighbor position is not accessible + by the DigiEvent. """ - columns, rows = self.get_neighbors() - pha_dict = self.as_dict() fmt = f'%{pha_width}d' - cols = np.arange(columns.min(), columns.max() + 1) + cols = np.arange(self.column - 1, self.column + 2) num_cols = cols[-1] - cols[0] + 1 - rows = np.arange(rows.min(), rows.max() + 1) + rows = np.arange(self.row - 1, self.row + 2) big_space = space(2 * pha_width + 1) text = f'\n{big_space}' text += ''.join([fmt % col for col in cols]) text += f'\n{big_space}+{line(pha_width * num_cols)}\n' - for row in rows: - text += f' {fmt % row} |' + for r in rows: + text += f' {fmt % r} |' for col in cols: try: - pha = fmt % pha_dict[(col, row)] + if col == self.column and r ==self.row: + pha = fmt % max(self.pha) + else: + pha = ' ' * pha_width except KeyError: pha = ' ' * pha_width text += pha text += f'\n{big_space}|\n' return text - ''' + diff --git a/hexsample/readout.py b/hexsample/readout.py index 171c1f9..c92d269 100644 --- a/hexsample/readout.py +++ b/hexsample/readout.py @@ -456,6 +456,7 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl # See: https://stackoverflow.com/questions/70094914/max-on-collections-counter coord_max = max(sparse_signal, key=sparse_signal.get) col_max, row_max = coord_max + print(coord_max) #... and converting it in ADC channel coordinates (value from 0 to 6)... adc_max = self.adc_channel(*coord_max) # ... creating a 7-elements array containing the PHA of the ADC channels from 0 to 6 @@ -476,4 +477,4 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl # Do not forget to update the trigger_id! self.trigger_id += 1 #The pha array is always in the order [pha(adc0), pha(adc1), pha(adc2), pha(adc3), pha(adc4), pha(adc5), pha(adc6)] - return DigiEventCircular(self.trigger_id, seconds, microseconds, livetime, pha, *coords) + return DigiEventCircular(self.trigger_id, seconds, microseconds, livetime, pha, *coord_max) diff --git a/tests/test_digi.py b/tests/test_digi.py index 0b2d42b..fe21a73 100644 --- a/tests/test_digi.py +++ b/tests/test_digi.py @@ -119,7 +119,7 @@ def test_digitization_circular(layout: HexagonalLayout = HexagonalLayout.ODD_R, x = np.append(x, np.full(num_pairs, x2)) y = np.append(y, np.full(num_pairs, y2)) event = readout.read(0., x, y, 100.) #this is a DigiEventCircular - #print(event.ascii()) + print(event.ascii()) def test_digi_event_rectangular(min_col: int = 106, max_col: int = 113, min_row: int = 15, max_row: int = 22, padding: Padding = Padding(1, 2, 3, 4)): From 688f2639e0f1b9c45934f5dff3f9223fb085a8d1 Mon Sep 17 00:00:00 2001 From: chiara Date: Mon, 22 Apr 2024 15:28:15 +0200 Subject: [PATCH 24/52] started implementing readout_strategy option for the simulation --- hexsample/app.py | 1 + hexsample/readout.py | 13 +++++++++++++ 2 files changed, 14 insertions(+) diff --git a/hexsample/app.py b/hexsample/app.py index 5e5f328..4cf82b0 100644 --- a/hexsample/app.py +++ b/hexsample/app.py @@ -164,6 +164,7 @@ def add_readout_options(self) -> None: help='trigger threshold in electron equivalent') group.add_argument('--zsupthreshold', type=int, default=0, help='zero-suppression threshold in ADC counts') + #group.add_argument('--readout_strategy', type=str, choices=) group.add_argument('--padding', type=int, nargs=4, default=(2, 2, 2, 2), help='padding on the four sides of the ROT') diff --git a/hexsample/readout.py b/hexsample/readout.py index c92d269..4bfec5a 100644 --- a/hexsample/readout.py +++ b/hexsample/readout.py @@ -21,6 +21,7 @@ """ from collections import Counter +from enum import Enum from typing import Tuple import numpy as np @@ -32,6 +33,18 @@ from hexsample import xpol +class HexagonalReadoutStrategy(Enum): + """Enum class expressing the possible readout strategies. + """ + + # Sparse readout strategy + SPARSE = 'SPARSE' + # Horizontal, pointy top, even rows are shoved right. + RECTANGULAR = 'RECTANGULAR' + # Vertical, flat top, odd columns are shoved down. + CIRCULAR = 'CIRCULAR' + + class HexagonalReadoutBase(HexagonalGrid): From c5f7013f5cc97773f2caa73ec022966b73e12d46 Mon Sep 17 00:00:00 2001 From: Luca Baldini Date: Mon, 22 Apr 2024 16:06:30 +0200 Subject: [PATCH 25/52] Minor. --- hexsample/app.py | 9 ++++++--- hexsample/bin/hxsim.py | 17 +++++++++++++---- hexsample/readout.py | 29 ++++++++++++++++++++++++----- 3 files changed, 43 insertions(+), 12 deletions(-) diff --git a/hexsample/app.py b/hexsample/app.py index 4cf82b0..a40f2d8 100644 --- a/hexsample/app.py +++ b/hexsample/app.py @@ -24,6 +24,7 @@ from hexsample import __pkgname__, __version__, __tagdate__, __url__ from hexsample.hexagon import HexagonalLayout +from hexsample.readout import HexagonalReadoutMode START_MESSAGE = f""" @@ -154,6 +155,11 @@ def add_readout_options(self) -> None: help='number of rows in the readout chip') group.add_argument('--pitch', type=float, default=0.005, help='pitch of the readout chip in cm') + modes = [item.value for item in HexagonalReadoutMode] + group.add_argument('--mode', type=str, choices=modes, default='RECTANGULAR', + help='readout mode') + group.add_argument('--padding', type=int, nargs=4, default=(2, 2, 2, 2), + help='padding on the four sides of the ROT') group.add_argument('--noise', type=float, default=20., help='equivalent noise charge rms in electrons') group.add_argument('--gain', type=float, default=1., @@ -164,9 +170,6 @@ def add_readout_options(self) -> None: help='trigger threshold in electron equivalent') group.add_argument('--zsupthreshold', type=int, default=0, help='zero-suppression threshold in ADC counts') - #group.add_argument('--readout_strategy', type=str, choices=) - group.add_argument('--padding', type=int, nargs=4, default=(2, 2, 2, 2), - help='padding on the four sides of the ROT') def add_sensor_options(self) -> None: """Add an option group for the sensor properties. diff --git a/hexsample/bin/hxsim.py b/hexsample/bin/hxsim.py index be9bce6..c5a1ede 100755 --- a/hexsample/bin/hxsim.py +++ b/hexsample/bin/hxsim.py @@ -29,7 +29,7 @@ from hexsample import rng from hexsample import HEXSAMPLE_DATA from hexsample.app import ArgumentParser -from hexsample.readout import HexagonalReadoutRectangular +from hexsample.readout import HexagonalReadoutMode, readout_chip from hexsample.fileio import DigiOutputFile from hexsample.hexagon import HexagonalLayout from hexsample.mc import PhotonList @@ -63,15 +63,24 @@ def hxsim(**kwargs): material = Material(kwargs['actmedium'], kwargs['fano']) sensor = Sensor(material, kwargs['thickness'], kwargs['transdiffsigma']) photon_list = PhotonList(source, sensor, kwargs['numevents']) + readout_mode = HexagonalReadoutMode(kwargs['mode']) + # Is there any nicer way to do this? See https://github.com/lucabaldini/hexsample/issues/51 + if readout_mode == HexagonalReadoutMode.SPARSE: + readout_args = kwargs['trgthreshold'], kwargs['zsupthreshold'], kwargs['offset'] + elif readout_mode == HexagonalReadoutMode.RECTANGULAR: + padding = Padding(*kwargs['padding']) + readout_args = kwargs['trgthreshold'], padding, kwargs['zsupthreshold'], kwargs['offset'] + elif readout_mode == HexagonalReadoutMode.CIRCULAR: + readout_args = kwargs['trgthreshold'], kwargs['zsupthreshold'], kwargs['offset'] + else: + raise RuntimeError args = HexagonalLayout(kwargs['layout']), kwargs['numcolumns'], kwargs['numrows'],\ kwargs['pitch'], kwargs['noise'], kwargs['gain'] - readout = HexagonalReadoutRectangular(*args) + readout = readout_chip(readout_mode, *args) logger.info(f'Readout chip: {readout}') output_file_path = kwargs.get('outfile') output_file = DigiOutputFile(output_file_path) output_file.update_header(**kwargs) - padding = Padding(*kwargs['padding']) - readout_args = kwargs['trgthreshold'], padding, kwargs['zsupthreshold'], kwargs['offset'] logger.info('Starting the event loop...') for mc_event in tqdm(photon_list): x, y = mc_event.propagate(sensor.trans_diffusion_sigma) diff --git a/hexsample/readout.py b/hexsample/readout.py index 4bfec5a..2c2385d 100644 --- a/hexsample/readout.py +++ b/hexsample/readout.py @@ -33,17 +33,17 @@ from hexsample import xpol -class HexagonalReadoutStrategy(Enum): +class HexagonalReadoutMode(Enum): """Enum class expressing the possible readout strategies. """ - # Sparse readout strategy + # Sparse readout strategy. SPARSE = 'SPARSE' - # Horizontal, pointy top, even rows are shoved right. + # Rectangular readout, a la XPOL. RECTANGULAR = 'RECTANGULAR' - # Vertical, flat top, odd columns are shoved down. + # Circular readout, with the highest pixel and the 6 neirest. CIRCULAR = 'CIRCULAR' - + class HexagonalReadoutBase(HexagonalGrid): @@ -491,3 +491,22 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl self.trigger_id += 1 #The pha array is always in the order [pha(adc0), pha(adc1), pha(adc2), pha(adc3), pha(adc4), pha(adc5), pha(adc6)] return DigiEventCircular(self.trigger_id, seconds, microseconds, livetime, pha, *coord_max) + + +# Mapping for the readout chip classes for each readout mode. +_READOUT_CLASS_DICT = { + HexagonalReadoutMode.SPARSE: HexagonalReadoutSparse, + HexagonalReadoutMode.RECTANGULAR: HexagonalReadoutRectangular, + HexagonalReadoutMode.CIRCULAR: HexagonalReadoutCircular +} + +def _readout_class(mode: HexagonalReadoutMode) -> type: + """Return the proper class to be used to instantiate a readout chip for a given + readout mode. + """ + return _READOUT_CLASS_DICT[mode] + +def readout_chip(mode: HexagonalReadoutMode, *args, **kwargs): + """Return an instance of the proper readout chip for a given readout mode. + """ + return _readout_class(mode)(*args, **kwargs) From 18de75e8a9d44578c9958398f40a6970d54453b6 Mon Sep 17 00:00:00 2001 From: chiara Date: Mon, 22 Apr 2024 18:06:44 +0200 Subject: [PATCH 26/52] started implementation of different DigiOutputFile(s) and started testing in tests/test_fileio.py --- hexsample/bin/hxsim.py | 4 +- hexsample/fileio.py | 293 ++++++++++++++++++++++++++++++++++------- hexsample/readout.py | 1 - tests/test_fileio.py | 12 +- 4 files changed, 254 insertions(+), 56 deletions(-) diff --git a/hexsample/bin/hxsim.py b/hexsample/bin/hxsim.py index c5a1ede..2c726f0 100755 --- a/hexsample/bin/hxsim.py +++ b/hexsample/bin/hxsim.py @@ -30,7 +30,7 @@ from hexsample import HEXSAMPLE_DATA from hexsample.app import ArgumentParser from hexsample.readout import HexagonalReadoutMode, readout_chip -from hexsample.fileio import DigiOutputFile +from hexsample.fileio import _digioutput_class from hexsample.hexagon import HexagonalLayout from hexsample.mc import PhotonList from hexsample.roi import Padding @@ -79,7 +79,7 @@ def hxsim(**kwargs): readout = readout_chip(readout_mode, *args) logger.info(f'Readout chip: {readout}') output_file_path = kwargs.get('outfile') - output_file = DigiOutputFile(output_file_path) + output_file = _digioutput_class(readout_mode)(output_file_path) output_file.update_header(**kwargs) logger.info('Starting the event loop...') for mc_event in tqdm(photon_list): diff --git a/hexsample/fileio.py b/hexsample/fileio.py index 6846800..2fd6bae 100644 --- a/hexsample/fileio.py +++ b/hexsample/fileio.py @@ -32,8 +32,8 @@ from hexsample import __version__, __tagdate__ from hexsample.mc import MonteCarloEvent -from hexsample.digi import DigiEventBase, DigiEventRectangular -from hexsample.readout import HexagonalReadoutCircular +from hexsample.digi import DigiEventBase, DigiEventSparse, DigiEventRectangular, DigiEventCircular +from hexsample.readout import HexagonalReadoutMode, HexagonalReadoutCircular from hexsample.recon import ReconEvent @@ -71,7 +71,9 @@ def _fill_mc_row(row: tables.tableextension.Row, event: MonteCarloEvent) -> None class DigiDescriptionBase(tables.IsDescription): - """ + """Base class for the description of the (flat) digi part of the file format. + It contains the trigger_id and time coordinates of the event, common to all + readout types. """ trigger_id = tables.Int32Col(pos=0) @@ -79,41 +81,53 @@ class DigiDescriptionBase(tables.IsDescription): microseconds = tables.Int32Col(pos=2) livetime = tables.Int32Col(pos=3) - def _fill_digi_row(row: tables.tableextension.Row, event: DigiEventBase) -> None: - """Helper function to fill an output table row, given a DigiEvent object. +def _fill_digi_row_base(row: tables.tableextension.Row, event: DigiEventBase) -> None: + """Helper function to fill an output table row, given a DigiEventBase object. - Note that this method of the base class is not calling the row.appen() hook, - which is delegated to the implementations in derived classes. + Note that this method of the base class is not calling the row.append() hook, + which is delegated to the implementations in derived classes. - .. note:: - This would have naturally belonged to the DigiDescription class as - a @staticmethod, but doing so is apparently breaking something into the - tables internals, and all of a sudden you get an exception due to the - fact that a staticmethod cannot be pickled. - """ - row['trigger_id'] = event.trigger_id - row['seconds'] = event.seconds - row['microseconds'] = event.microseconds - row['livetime'] = event.livetime + .. note:: + This would have naturally belonged to the DigiDescriptionBase class as + a @staticmethod, but doing so is apparently breaking something into the + tables internals, and all of a sudden you get an exception due to the + fact that a staticmethod cannot be pickled. + """ + row['trigger_id'] = event.trigger_id + row['seconds'] = event.seconds + row['microseconds'] = event.microseconds + row['livetime'] = event.livetime class DigiDescriptionSparse(DigiDescriptionBase): + """Description of the (flat) digi part of the file format for a sparse readout + DigiEvent. """ - """ + #columns = tables.Int16Col(pos=4) + #rows = tables.Int16Col(pos=5) - def _fill_digi_row(row: tables.tableextension.Row, event: DigiEventBase) -> None: - """Overloaded method. - """ - DigiDescriptionBase._fill_digi_row(row, event) - row.append() +def _fill_digi_row_sparse(row: tables.tableextension.Row, event: DigiEventBase) -> None: + """Overloaded method. + It uses the _fill_digi_row_base() function for filling the trigger_id and time + coordinates of the event. + + .. note:: + This would have naturally belonged to the DigiDescriptionSparse class as + a @staticmethod, but doing so is apparently breaking something into the + tables internals, and all of a sudden you get an exception due to the + fact that a staticmethod cannot be pickled. + """ + _fill_digi_row_base(row, event) + row.append() class DigiDescriptionRectangular(DigiDescriptionBase): - """ + """Description of the (flat) digi part of the file format for a rectangular readout + DigiEvent. """ min_col = tables.Int16Col(pos=4) @@ -125,41 +139,56 @@ class DigiDescriptionRectangular(DigiDescriptionBase): padding_bottom = tables.Int8Col(pos=10) padding_left = tables.Int8Col(pos=11) - def _fill_digi_row(row: tables.tableextension.Row, event: DigiEventBase) -> None: - """Overloaded method. - """ - DigiDescriptionBase._fill_digi_row(row, event) - row['min_col'] = event.roi.min_col - row['max_col'] = event.roi.max_col - row['min_row'] = event.roi.min_row - row['max_row'] = event.roi.max_row - row['padding_top'] = event.roi.padding.top - row['padding_right'] = event.roi.padding.right - row['padding_bottom'] = event.roi.padding.bottom - row['padding_left'] = event.roi.padding.left - row.append() +def _fill_digi_row_rectangular(row: tables.tableextension.Row, event: DigiEventBase) -> None: + """Overloaded method. + It uses the _fill_digi_row_base() function for filling the trigger_id and time + coordinates of the event. + + .. note:: + This would have naturally belonged to the DigiDescriptionRectangular class as + a @staticmethod, but doing so is apparently breaking something into the + tables internals, and all of a sudden you get an exception due to the + fact that a staticmethod cannot be pickled. + """ + _fill_digi_row_base(row, event) + row['min_col'] = event.roi.min_col + row['max_col'] = event.roi.max_col + row['min_row'] = event.roi.min_row + row['max_row'] = event.roi.max_row + row['padding_top'] = event.roi.padding.top + row['padding_right'] = event.roi.padding.right + row['padding_bottom'] = event.roi.padding.bottom + row['padding_left'] = event.roi.padding.left + row.append() class DigiDescriptionCircular(DigiDescriptionBase): - """ + """Description of the (flat) digi part of the file format for a rectangular readout + DigiEvent. """ column = tables.Int16Col(pos=4) row = tables.Int16Col(pos=5) pha = tables.Int16Col(shape=HexagonalReadoutCircular.NUM_PIXELS, pos=6) - def _fill_digi_row(row: tables.tableextension.Row, event: DigiEventBase) -> None: - """Overloaded method. - """ - DigiDescriptionBase._fill_digi_row(row, event) - row['column'] = event.column - row['row'] = event.row - row['pha'] = event.pha - row.append() - +def _fill_digi_row_circular(row: tables.tableextension.Row, event: DigiEventBase) -> None: + """Overloaded method. + It uses the _fill_digi_row_base() function for filling the trigger_id and time + coordinates of the event. + .. note:: + This would have naturally belonged to the DigiDescriptionCircular class as + a @staticmethod, but doing so is apparently breaking something into the + tables internals, and all of a sudden you get an exception due to the + fact that a staticmethod cannot be pickled. + """ + _fill_digi_row_base(row, event) + row['column'] = event.column + row['row'] = event.row + row['pha'] = event.pha + row.append() class DigiDescription(tables.IsDescription): @@ -335,9 +364,177 @@ def flush(self) -> None: This needs to be reimplemented in derived classes. """ raise NotImplementedError + +class DigiOutputFileSparse(OutputFileBase): + + """Description of a sparse digitized output file. + + This can represent either a digitized files written by the DAQ, or the + output of a simulation, in which case it contains an additional group and + table for the Monte Carlo ground truth. + + Arguments + --------- + file_path : str + The path to the file on disk. + """ + + _FILE_TYPE = FileType.DIGI + DIGI_TABLE_SPECS = ('digi_table', DigiDescriptionSparse, 'Digi data') + PHA_ARRAY_SPECS = ('pha', tables.Int32Atom(shape=())) + MC_TABLE_SPECS = ('mc_table', MonteCarloDescription, 'Monte Carlo data') + + def __init__(self, file_path: str): + """Constructor. + """ + super().__init__(file_path) + self.digi_group = self.create_group(self.root, 'digi', 'Digi') + self.digi_table = self.create_table(self.digi_group, *self.DIGI_TABLE_SPECS) + self.pha_array = self.create_vlarray(self.digi_group, *self.PHA_ARRAY_SPECS) + self.mc_group = self.create_group(self.root, 'mc', 'Monte Carlo') + self.mc_table = self.create_table(self.mc_group, *self.MC_TABLE_SPECS) + + def add_row(self, digi_event: DigiEventSparse, mc_event: MonteCarloEvent) -> None: + """Add one row to the file. + + Arguments + --------- + digi : DigiEventSparse + The digitized event contribution. + + mc : MonteCarloEvent + The Monte Carlo event contribution. + """ + # pylint: disable=arguments-differ + _fill_digi_row_sparse(self.digi_table.row, digi_event) + self.pha_array.append(digi_event.pha.flatten()) + _fill_mc_row(self.mc_table.row, mc_event) + + def flush(self) -> None: + """Flush the basic file components. + """ + self.digi_table.flush() + self.pha_array.flush() + self.mc_table.flush() + +class DigiOutputFileRectangular(OutputFileBase): + """Description of a rectangular digitized output file. + This can represent either a digitized files written by the DAQ, or the + output of a simulation, in which case it contains an additional group and + table for the Monte Carlo ground truth. + Arguments + --------- + file_path : str + The path to the file on disk. + """ + + _FILE_TYPE = FileType.DIGI + DIGI_TABLE_SPECS = ('digi_table', DigiDescriptionRectangular, 'Digi data') + PHA_ARRAY_SPECS = ('pha', tables.Int32Atom(shape=())) + MC_TABLE_SPECS = ('mc_table', MonteCarloDescription, 'Monte Carlo data') + + def __init__(self, file_path: str): + """Constructor. + """ + super().__init__(file_path) + self.digi_group = self.create_group(self.root, 'digi', 'Digi') + self.digi_table = self.create_table(self.digi_group, *self.DIGI_TABLE_SPECS) + self.pha_array = self.create_vlarray(self.digi_group, *self.PHA_ARRAY_SPECS) + self.mc_group = self.create_group(self.root, 'mc', 'Monte Carlo') + self.mc_table = self.create_table(self.mc_group, *self.MC_TABLE_SPECS) + + def add_row(self, digi_event: DigiEventRectangular, mc_event: MonteCarloEvent) -> None: + """Add one row to the file. + + Arguments + --------- + digi : DigiEventRectangular + The digitized event contribution. + + mc : MonteCarloEvent + The Monte Carlo event contribution. + """ + # pylint: disable=arguments-differ + _fill_digi_row_rectangular(self.digi_table.row, digi_event) + self.pha_array.append(digi_event.pha.flatten()) + _fill_mc_row(self.mc_table.row, mc_event) + + def flush(self) -> None: + """Flush the basic file components. + """ + self.digi_table.flush() + self.pha_array.flush() + self.mc_table.flush() + +class DigiOutputFileCircular(OutputFileBase): + + """Description of a circular digitized output file. + + This can represent either a digitized files written by the DAQ, or the + output of a simulation, in which case it contains an additional group and + table for the Monte Carlo ground truth. + + Arguments + --------- + file_path : str + The path to the file on disk. + """ + + _FILE_TYPE = FileType.DIGI + DIGI_TABLE_SPECS = ('digi_table', DigiDescriptionCircular, 'Digi data') + PHA_ARRAY_SPECS = ('pha', tables.Int32Atom(shape=())) + MC_TABLE_SPECS = ('mc_table', MonteCarloDescription, 'Monte Carlo data') + + def __init__(self, file_path: str): + """Constructor. + """ + super().__init__(file_path) + self.digi_group = self.create_group(self.root, 'digi', 'Digi') + self.digi_table = self.create_table(self.digi_group, *self.DIGI_TABLE_SPECS) + self.pha_array = self.create_vlarray(self.digi_group, *self.PHA_ARRAY_SPECS) + self.mc_group = self.create_group(self.root, 'mc', 'Monte Carlo') + self.mc_table = self.create_table(self.mc_group, *self.MC_TABLE_SPECS) + + def add_row(self, digi_event: DigiEventCircular, mc_event: MonteCarloEvent) -> None: + """Add one row to the file. + + Arguments + --------- + digi : DigiEventCircular + The digitized event contribution. + + mc : MonteCarloEvent + The Monte Carlo event contribution. + """ + # pylint: disable=arguments-differ + _fill_digi_row_circular(self.digi_table.row, digi_event) + self.pha_array.append(digi_event.pha.flatten()) + _fill_mc_row(self.mc_table.row, mc_event) + + def flush(self) -> None: + """Flush the basic file components. + """ + self.digi_table.flush() + self.pha_array.flush() + self.mc_table.flush() + +# Mapping for the digi description classes for each readout mode. +_FILEIO_CLASS_DICT = { + HexagonalReadoutMode.SPARSE: DigiOutputFileSparse, + HexagonalReadoutMode.RECTANGULAR: DigiOutputFileRectangular, + HexagonalReadoutMode.CIRCULAR: DigiOutputFileCircular +} + +def _digioutput_class(mode: HexagonalReadoutMode) -> type: + """Return the proper class to be used as DigiOutputFile, depending on the + readout mode of the chip. + """ + return _FILEIO_CLASS_DICT[mode] + +''' class DigiOutputFile(OutputFileBase): """Description of a digitized output file. @@ -390,7 +587,7 @@ def flush(self) -> None: self.pha_array.flush() self.mc_table.flush() - +''' class ReconOutputFile(OutputFileBase): diff --git a/hexsample/readout.py b/hexsample/readout.py index 2c2385d..471f5c6 100644 --- a/hexsample/readout.py +++ b/hexsample/readout.py @@ -469,7 +469,6 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl # See: https://stackoverflow.com/questions/70094914/max-on-collections-counter coord_max = max(sparse_signal, key=sparse_signal.get) col_max, row_max = coord_max - print(coord_max) #... and converting it in ADC channel coordinates (value from 0 to 6)... adc_max = self.adc_channel(*coord_max) # ... creating a 7-elements array containing the PHA of the ADC channels from 0 to 6 diff --git a/tests/test_fileio.py b/tests/test_fileio.py index 33881f2..52042f9 100644 --- a/tests/test_fileio.py +++ b/tests/test_fileio.py @@ -21,8 +21,9 @@ from hexsample import HEXSAMPLE_DATA from hexsample.digi import DigiEventRectangular -from hexsample.fileio import DigiInputFile, DigiOutputFile, ReconInputFile, ReconOutputFile,\ - FileType, peek_file_type, open_input_file +from hexsample.readout import HexagonalReadoutMode +from hexsample.fileio import DigiInputFile, ReconInputFile, ReconOutputFile,\ + FileType, peek_file_type, open_input_file, _FILEIO_CLASS_DICT, _digioutput_class from hexsample.mc import MonteCarloEvent from hexsample.roi import RegionOfInterest, Padding @@ -42,8 +43,9 @@ def _digi_event(index : int) -> DigiEventRectangular: def _test_write(file_path, num_events : int = 10): """Small test writing a bunch of toy event strcutures to file. + 10 events for every readout methods are created. """ - output_file = DigiOutputFile(file_path) + output_file = _digioutput_class(HexagonalReadoutMode.RECTANGULAR)(file_path) for i in range(num_events): output_file.add_row(_digi_event(i), _mc_event(i)) output_file.close() @@ -78,8 +80,8 @@ def test_file_type(): """ # Test for the digi files. file_path = HEXSAMPLE_DATA / 'test_digi_filetype.h5' - digi_file = DigiOutputFile(file_path) - digi_file.close() + #digi_file = DigiOutputFile(file_path) + #digi_file.close() digi_file = DigiInputFile(file_path) assert digi_file.file_type == FileType.DIGI digi_file.close() From a63b4599d2a864e1ce368314a56b2fb6126328d7 Mon Sep 17 00:00:00 2001 From: chiara Date: Tue, 23 Apr 2024 10:52:29 +0200 Subject: [PATCH 27/52] continued implementing I/O for the different modes of readout --- hexsample/digi.py | 25 +++++- hexsample/fileio.py | 205 +++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 226 insertions(+), 4 deletions(-) diff --git a/hexsample/digi.py b/hexsample/digi.py index ea0353e..d1002fd 100644 --- a/hexsample/digi.py +++ b/hexsample/digi.py @@ -121,6 +121,17 @@ def __post_init__(self) -> None: if not len(self.rows) == len(self.columns) == len(self.pha): raise RuntimeError(f'{self.__class__.__name__} has {len(self.rows)} rows' f', {len(self.columns)} columns, and {len(self.pha)} PHA values') + + @classmethod + def from_digi(cls, file_row: np.ndarray, pha: np.ndarray): + """Alternative constructor rebuilding an object from a row on a digi file. + + This is used internally when we access event data in a digi file, and + we need to reassemble a DigiEvent object from a given row of a digi table. + """ + # pylint: disable=too-many-locals + trigger_id, seconds, microseconds, livetime, pha, columns, rows = file_row + return cls(trigger_id, seconds, microseconds, livetime, pha, columns, rows) def as_dict(self) -> dict: """Return the pixel content of the event in the form of a {(col, row): pha} @@ -190,7 +201,7 @@ def __eq__(self, other) -> bool: return DigiEventBase.__eq__(self, other) and self.roi == other.roi @classmethod - def from_digi(cls, row: np.ndarray, pha: np.ndarray): + def from_digi(cls, file_row: np.ndarray, pha: np.ndarray): """Alternative constructor rebuilding an object from a row on a digi file. This is used internally when we access event data in a digi file, and @@ -198,7 +209,7 @@ def from_digi(cls, row: np.ndarray, pha: np.ndarray): """ # pylint: disable=too-many-locals trigger_id, seconds, microseconds, livetime, min_col, max_col, min_row, max_row,\ - pad_top, pad_right, pad_bottom, pad_left = row + pad_top, pad_right, pad_bottom, pad_left = file_row padding = Padding(pad_top, pad_right, pad_bottom, pad_left) roi = RegionOfInterest(min_col, max_col, min_row, max_row, padding) return cls(trigger_id, seconds, microseconds, livetime, pha, roi) @@ -296,6 +307,16 @@ def __post_init__(self) -> None: raise RuntimeError(f'{self.__class__.__name__} has {len(self.pha)} PHA values' f'instead of {7}.') + @classmethod + def from_digi(cls, file_row: np.ndarray, pha: np.ndarray): + """Alternative constructor rebuilding an object from a row on a digi file. + + This is used internally when we access event data in a digi file, and + we need to reassemble a DigiEvent object from a given row of a digi table. + """ + # pylint: disable=too-many-locals + trigger_id, seconds, microseconds, livetime, pha, column, row = file_row + return cls(trigger_id, seconds, microseconds, livetime, pha, column, row) def ascii(self, pha_width: int = 5) -> str: """Ascii representation. diff --git a/hexsample/fileio.py b/hexsample/fileio.py index 2fd6bae..c46dd9a 100644 --- a/hexsample/fileio.py +++ b/hexsample/fileio.py @@ -237,7 +237,8 @@ def _fill_digi_row(row: tables.tableextension.Row, event: DigiEventBase) -> None class ReconDescription(tables.IsDescription): - """Description of the recon file format. + """Description of the recon file format. This should be common to all the + modes of readout, so it is the same aside from the DigiDescription type. """ # pylint: disable=too-few-public-methods @@ -591,7 +592,8 @@ def flush(self) -> None: class ReconOutputFile(OutputFileBase): - """Description of a reconstructed output file. + """Description of a reconstructed output file. This should be the same for + all types of DigiEvent. Arguments --------- @@ -676,6 +678,205 @@ def header_value(self, key: str, default: Any = None) -> Any: return self.header.get(key, default) +class DigiInputFileSparse(InputFileBase): + + """Description of a sparse digitized input file. + + This has a very simple interface: we cache references to the relevant tables + when we open the file and we provide methods to reassemble a specific table + row into the corresponding DigiEvent or MonteCarloEvent objects, along with + an implementation of the iterator protocol to make event loops easier. + """ + + def __init__(self, file_path: str): + """Constructor. + """ + super().__init__(file_path) + self.digi_table = self.root.digi.digi_table + self.pha_array = self.root.digi.pha + self.mc_table = self.root.mc.mc_table + self.__index = -1 + + def column(self, name: str) -> np.ndarray: + """Return a given column in the digi table. + """ + return self.digi_table.col(name) + + def mc_column(self, name: str) -> np.ndarray: + """Return a given column in the Monte Carlo table. + """ + return self.mc_table.col(name) + + def digi_event(self, row_index: int) -> DigiEventSparse: + """Random access to the DigiEventSparse part of the event contribution. + + Arguments + --------- + row_index : int + The index of the target row in the event file. + """ + row = self.digi_table[row_index] + pha = self.pha_array[row_index] + return DigiEventSparse.from_digi(row, pha) + + def mc_event(self, row_index: int) -> MonteCarloEvent: + """Random access to the MonteCarloEvent part of the event contribution. + + Arguments + --------- + row_index : int + The index of the target row in the event file. + """ + row = self.mc_table[row_index] + return MonteCarloEvent(*row) + + def __iter__(self): + """Overloaded method for the implementation of the iterator protocol. + """ + self.__index = -1 + return self + + def __next__(self) -> DigiEventSparse: + """Overloaded method for the implementation of the iterator protocol. + """ + self.__index += 1 + if self.__index == len(self.digi_table): + raise StopIteration + return self.digi_event(self.__index) + +class DigiInputFileRectangular(InputFileBase): + + """Description of a rectangular digitized input file. + + This has a very simple interface: we cache references to the relevant tables + when we open the file and we provide methods to reassemble a specific table + row into the corresponding DigiEvent or MonteCarloEvent objects, along with + an implementation of the iterator protocol to make event loops easier. + """ + + def __init__(self, file_path: str): + """Constructor. + """ + super().__init__(file_path) + self.digi_table = self.root.digi.digi_table + self.pha_array = self.root.digi.pha + self.mc_table = self.root.mc.mc_table + self.__index = -1 + + def column(self, name: str) -> np.ndarray: + """Return a given column in the digi table. + """ + return self.digi_table.col(name) + + def mc_column(self, name: str) -> np.ndarray: + """Return a given column in the Monte Carlo table. + """ + return self.mc_table.col(name) + + def digi_event(self, row_index: int) -> DigiEventRectangular: + """Random access to the DigiEventSparse part of the event contribution. + + Arguments + --------- + row_index : int + The index of the target row in the event file. + """ + row = self.digi_table[row_index] + pha = self.pha_array[row_index] + return DigiEventRectangular.from_digi(row, pha) + + def mc_event(self, row_index: int) -> MonteCarloEvent: + """Random access to the MonteCarloEvent part of the event contribution. + + Arguments + --------- + row_index : int + The index of the target row in the event file. + """ + row = self.mc_table[row_index] + return MonteCarloEvent(*row) + + def __iter__(self): + """Overloaded method for the implementation of the iterator protocol. + """ + self.__index = -1 + return self + + def __next__(self) -> DigiEventRectangular: + """Overloaded method for the implementation of the iterator protocol. + """ + self.__index += 1 + if self.__index == len(self.digi_table): + raise StopIteration + return self.digi_event(self.__index) + +class DigiInputFileCircular(InputFileBase): + + """Description of a circular digitized input file. + + This has a very simple interface: we cache references to the relevant tables + when we open the file and we provide methods to reassemble a specific table + row into the corresponding DigiEvent or MonteCarloEvent objects, along with + an implementation of the iterator protocol to make event loops easier. + """ + + def __init__(self, file_path: str): + """Constructor. + """ + super().__init__(file_path) + self.digi_table = self.root.digi.digi_table + self.pha_array = self.root.digi.pha + self.mc_table = self.root.mc.mc_table + self.__index = -1 + + def column(self, name: str) -> np.ndarray: + """Return a given column in the digi table. + """ + return self.digi_table.col(name) + + def mc_column(self, name: str) -> np.ndarray: + """Return a given column in the Monte Carlo table. + """ + return self.mc_table.col(name) + + def digi_event(self, row_index: int) -> DigiEventCircular: + """Random access to the DigiEventSparse part of the event contribution. + + Arguments + --------- + row_index : int + The index of the target row in the event file. + """ + row = self.digi_table[row_index] + pha = self.pha_array[row_index] + return DigiEventCircular.from_digi(row, pha) + + def mc_event(self, row_index: int) -> MonteCarloEvent: + """Random access to the MonteCarloEvent part of the event contribution. + + Arguments + --------- + row_index : int + The index of the target row in the event file. + """ + row = self.mc_table[row_index] + return MonteCarloEvent(*row) + + def __iter__(self): + """Overloaded method for the implementation of the iterator protocol. + """ + self.__index = -1 + return self + + def __next__(self) -> DigiEventCircular: + """Overloaded method for the implementation of the iterator protocol. + """ + self.__index += 1 + if self.__index == len(self.digi_table): + raise StopIteration + return self.digi_event(self.__index) + + class DigiInputFile(InputFileBase): From fff434c493bc87f75a3ce752df36e01e02c29c23 Mon Sep 17 00:00:00 2001 From: chiara Date: Tue, 23 Apr 2024 11:42:57 +0200 Subject: [PATCH 28/52] continued implementing I/O for different readout types --- hexsample/fileio.py | 4 +++- tests/test_fileio.py | 23 +++++++++++++++++++---- 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/hexsample/fileio.py b/hexsample/fileio.py index c46dd9a..f247c31 100644 --- a/hexsample/fileio.py +++ b/hexsample/fileio.py @@ -391,6 +391,8 @@ def __init__(self, file_path: str): super().__init__(file_path) self.digi_group = self.create_group(self.root, 'digi', 'Digi') self.digi_table = self.create_table(self.digi_group, *self.DIGI_TABLE_SPECS) + self.columns_array = self.create_vlarray(self.digi_group, *self.PHA_ARRAY_SPECS) + self.rows_array = self.create_vlarray(self.digi_group, *self.PHA_ARRAY_SPECS) self.pha_array = self.create_vlarray(self.digi_group, *self.PHA_ARRAY_SPECS) self.mc_group = self.create_group(self.root, 'mc', 'Monte Carlo') self.mc_table = self.create_table(self.mc_group, *self.MC_TABLE_SPECS) @@ -495,7 +497,7 @@ def __init__(self, file_path: str): super().__init__(file_path) self.digi_group = self.create_group(self.root, 'digi', 'Digi') self.digi_table = self.create_table(self.digi_group, *self.DIGI_TABLE_SPECS) - self.pha_array = self.create_vlarray(self.digi_group, *self.PHA_ARRAY_SPECS) + self.pha_array = self.create_array(self.digi_group, *self.PHA_ARRAY_SPECS) self.mc_group = self.create_group(self.root, 'mc', 'Monte Carlo') self.mc_table = self.create_table(self.mc_group, *self.MC_TABLE_SPECS) diff --git a/tests/test_fileio.py b/tests/test_fileio.py index 52042f9..a7bdde9 100644 --- a/tests/test_fileio.py +++ b/tests/test_fileio.py @@ -20,7 +20,7 @@ import numpy as np from hexsample import HEXSAMPLE_DATA -from hexsample.digi import DigiEventRectangular +from hexsample.digi import DigiEventSparse, DigiEventRectangular, DigiEventCircular from hexsample.readout import HexagonalReadoutMode from hexsample.fileio import DigiInputFile, ReconInputFile, ReconOutputFile,\ FileType, peek_file_type, open_input_file, _FILEIO_CLASS_DICT, _digioutput_class @@ -34,20 +34,35 @@ def _mc_event(index : int) -> MonteCarloEvent: """ return MonteCarloEvent(0.1 * index, 5.9, 0., 0., 0.02, 1000 + index) -def _digi_event(index : int) -> DigiEventRectangular: +def _digi_event_rectangular(index : int) -> DigiEventRectangular: """Create a bogus DigiEvent object with index-dependent properties. """ roi = RegionOfInterest(100, 107, 150, 155 + index * 2, Padding(2)) pha = np.full(roi.size, index) return DigiEventRectangular(index, index, index, 0, pha, roi) +def _digi_event_sparse(index : int) -> DigiEventSparse: + columns = np.randint(0, high=10, size=10) + rows = np.randint(0, high=10, size=10) + pha = np.full(10, index) + return DigiEventSparse(columns, rows, pha) + +def _test_write_sparse(file_path, num_events : int = 10): + """Small test writing a bunch of toy event strcutures to file. + 10 events for every readout methods are created. + """ + output_file = _digioutput_class(HexagonalReadoutMode.RECTANGULAR)(file_path) + for i in range(num_events): + output_file.add_row(_digi_event_rectangular(i), _mc_event(i)) + output_file.close() + def _test_write(file_path, num_events : int = 10): """Small test writing a bunch of toy event strcutures to file. 10 events for every readout methods are created. """ output_file = _digioutput_class(HexagonalReadoutMode.RECTANGULAR)(file_path) for i in range(num_events): - output_file.add_row(_digi_event(i), _mc_event(i)) + output_file.add_row(_digi_event_rectangular(i), _mc_event(i)) output_file.close() def _test_read(file_path): @@ -58,7 +73,7 @@ def _test_read(file_path): for i, event in enumerate(input_file): print(event.ascii()) print(input_file.mc_event(i)) - target = _digi_event(i) + target = _digi_event_rectangular(i) assert event.trigger_id == target.trigger_id assert event.seconds == target.seconds assert event.microseconds == target.microseconds From eaf27930ec5d46b0bced7460998c577e21e1ca40 Mon Sep 17 00:00:00 2001 From: chiara Date: Tue, 23 Apr 2024 12:08:28 +0200 Subject: [PATCH 29/52] Continued implementing classes for file I/O --- hexsample/fileio.py | 27 ++++++++++++++++----------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/hexsample/fileio.py b/hexsample/fileio.py index f247c31..143400a 100644 --- a/hexsample/fileio.py +++ b/hexsample/fileio.py @@ -105,13 +105,11 @@ class DigiDescriptionSparse(DigiDescriptionBase): """Description of the (flat) digi part of the file format for a sparse readout DigiEvent. """ - #columns = tables.Int16Col(pos=4) - #rows = tables.Int16Col(pos=5) def _fill_digi_row_sparse(row: tables.tableextension.Row, event: DigiEventBase) -> None: """Overloaded method. It uses the _fill_digi_row_base() function for filling the trigger_id and time - coordinates of the event. + coordinates of the event. .. note:: This would have naturally belonged to the DigiDescriptionSparse class as @@ -162,7 +160,6 @@ def _fill_digi_row_rectangular(row: tables.tableextension.Row, event: DigiEventB row.append() - class DigiDescriptionCircular(DigiDescriptionBase): """Description of the (flat) digi part of the file format for a rectangular readout @@ -193,6 +190,8 @@ def _fill_digi_row_circular(row: tables.tableextension.Row, event: DigiEventBase class DigiDescription(tables.IsDescription): """Description of the (flat) digi part of the file format. + NOTE: This should be eliminated when the above three classes will be fully + implemented and tested. """ # pylint: disable=too-few-public-methods @@ -275,6 +274,7 @@ def _fill_recon_row(row: tables.tableextension.Row, event: ReconEvent) -> None: class FileType(Enum): """Enum class for the different file types. + ****** IS IT POSSIBLE TO DEFINE SUBCLASSES FOR DIGI? ****** """ DIGI = 'Digi' @@ -381,8 +381,11 @@ class DigiOutputFileSparse(OutputFileBase): """ _FILE_TYPE = FileType.DIGI + _READOUT_TYPE = HexagonalReadoutMode.SPARSE #not sure if useful DIGI_TABLE_SPECS = ('digi_table', DigiDescriptionSparse, 'Digi data') - PHA_ARRAY_SPECS = ('pha', tables.Int32Atom(shape=())) + COLUMNS_ARRAY_SPECS = ('columns', tables.Int16Atom(shape=())) + ROWS_ARRAY_SPECS = ('rows', tables.Int16Atom(shape=())) + PHA_ARRAY_SPECS = ('pha', tables.Int16Atom(shape=())) MC_TABLE_SPECS = ('mc_table', MonteCarloDescription, 'Monte Carlo data') def __init__(self, file_path: str): @@ -391,8 +394,8 @@ def __init__(self, file_path: str): super().__init__(file_path) self.digi_group = self.create_group(self.root, 'digi', 'Digi') self.digi_table = self.create_table(self.digi_group, *self.DIGI_TABLE_SPECS) - self.columns_array = self.create_vlarray(self.digi_group, *self.PHA_ARRAY_SPECS) - self.rows_array = self.create_vlarray(self.digi_group, *self.PHA_ARRAY_SPECS) + self.columns_array = self.create_vlarray(self.digi_group, *self.COLUMNS_ARRAY_SPECS) + self.rows_array = self.create_vlarray(self.digi_group, *self.ROWS_ARRAY_SPECS) self.pha_array = self.create_vlarray(self.digi_group, *self.PHA_ARRAY_SPECS) self.mc_group = self.create_group(self.root, 'mc', 'Monte Carlo') self.mc_table = self.create_table(self.mc_group, *self.MC_TABLE_SPECS) @@ -410,6 +413,8 @@ def add_row(self, digi_event: DigiEventSparse, mc_event: MonteCarloEvent) -> Non """ # pylint: disable=arguments-differ _fill_digi_row_sparse(self.digi_table.row, digi_event) + self.columns_array.append(digi_event.columns.flatten()) + self.rows_array.append(digi_event.rows.flatten()) self.pha_array.append(digi_event.pha.flatten()) _fill_mc_row(self.mc_table.row, mc_event) @@ -417,6 +422,8 @@ def flush(self) -> None: """Flush the basic file components. """ self.digi_table.flush() + self.columns_array.flush() + self.rows_array.flush() self.pha_array.flush() self.mc_table.flush() @@ -435,6 +442,7 @@ class DigiOutputFileRectangular(OutputFileBase): """ _FILE_TYPE = FileType.DIGI + _READOUT_TYPE = HexagonalReadoutMode.RECTANGULAR #not sure if useful DIGI_TABLE_SPECS = ('digi_table', DigiDescriptionRectangular, 'Digi data') PHA_ARRAY_SPECS = ('pha', tables.Int32Atom(shape=())) MC_TABLE_SPECS = ('mc_table', MonteCarloDescription, 'Monte Carlo data') @@ -487,8 +495,8 @@ class DigiOutputFileCircular(OutputFileBase): """ _FILE_TYPE = FileType.DIGI + _READOUT_TYPE = HexagonalReadoutMode.CIRCULAR #not sure if useful DIGI_TABLE_SPECS = ('digi_table', DigiDescriptionCircular, 'Digi data') - PHA_ARRAY_SPECS = ('pha', tables.Int32Atom(shape=())) MC_TABLE_SPECS = ('mc_table', MonteCarloDescription, 'Monte Carlo data') def __init__(self, file_path: str): @@ -497,7 +505,6 @@ def __init__(self, file_path: str): super().__init__(file_path) self.digi_group = self.create_group(self.root, 'digi', 'Digi') self.digi_table = self.create_table(self.digi_group, *self.DIGI_TABLE_SPECS) - self.pha_array = self.create_array(self.digi_group, *self.PHA_ARRAY_SPECS) self.mc_group = self.create_group(self.root, 'mc', 'Monte Carlo') self.mc_table = self.create_table(self.mc_group, *self.MC_TABLE_SPECS) @@ -514,14 +521,12 @@ def add_row(self, digi_event: DigiEventCircular, mc_event: MonteCarloEvent) -> N """ # pylint: disable=arguments-differ _fill_digi_row_circular(self.digi_table.row, digi_event) - self.pha_array.append(digi_event.pha.flatten()) _fill_mc_row(self.mc_table.row, mc_event) def flush(self) -> None: """Flush the basic file components. """ self.digi_table.flush() - self.pha_array.flush() self.mc_table.flush() # Mapping for the digi description classes for each readout mode. From 7f6d16cdd34574964f703c80d039d3d1d0653047 Mon Sep 17 00:00:00 2001 From: chiara Date: Tue, 23 Apr 2024 15:58:34 +0200 Subject: [PATCH 30/52] Continued implementing I/O facilities --- hexsample/digi.py | 7 ++- hexsample/fileio.py | 76 ++++++++++++++++++-------- tests/test_fileio.py | 124 ++++++++++++++++++++++++++++++++++++------- 3 files changed, 163 insertions(+), 44 deletions(-) diff --git a/hexsample/digi.py b/hexsample/digi.py index d1002fd..25328b0 100644 --- a/hexsample/digi.py +++ b/hexsample/digi.py @@ -26,7 +26,6 @@ from loguru import logger import numpy as np -#from hexsample.readout import HexagonalReadoutCircular from hexsample.pprint import AnsiFontEffect, ansi_format, space, line from hexsample.roi import Padding, RegionOfInterest @@ -123,14 +122,14 @@ def __post_init__(self) -> None: f', {len(self.columns)} columns, and {len(self.pha)} PHA values') @classmethod - def from_digi(cls, file_row: np.ndarray, pha: np.ndarray): + def from_digi(cls, file_row: np.ndarray, pha: np.ndarray, columns: np.ndarray, rows: np.ndarray): """Alternative constructor rebuilding an object from a row on a digi file. This is used internally when we access event data in a digi file, and we need to reassemble a DigiEvent object from a given row of a digi table. """ # pylint: disable=too-many-locals - trigger_id, seconds, microseconds, livetime, pha, columns, rows = file_row + trigger_id, seconds, microseconds, livetime = file_row return cls(trigger_id, seconds, microseconds, livetime, pha, columns, rows) def as_dict(self) -> dict: @@ -308,7 +307,7 @@ def __post_init__(self) -> None: f'instead of {7}.') @classmethod - def from_digi(cls, file_row: np.ndarray, pha: np.ndarray): + def from_digi(cls, file_row: np.ndarray): """Alternative constructor rebuilding an object from a row on a digi file. This is used internally when we access event data in a digi file, and diff --git a/hexsample/fileio.py b/hexsample/fileio.py index 143400a..f995876 100644 --- a/hexsample/fileio.py +++ b/hexsample/fileio.py @@ -165,10 +165,9 @@ class DigiDescriptionCircular(DigiDescriptionBase): """Description of the (flat) digi part of the file format for a rectangular readout DigiEvent. """ - - column = tables.Int16Col(pos=4) - row = tables.Int16Col(pos=5) - pha = tables.Int16Col(shape=HexagonalReadoutCircular.NUM_PIXELS, pos=6) + pha = tables.Int16Col(shape=HexagonalReadoutCircular.NUM_PIXELS, pos=4) + column = tables.Int16Col(pos=5) + row = tables.Int16Col(pos=6) def _fill_digi_row_circular(row: tables.tableextension.Row, event: DigiEventBase) -> None: """Overloaded method. @@ -182,9 +181,9 @@ def _fill_digi_row_circular(row: tables.tableextension.Row, event: DigiEventBase fact that a staticmethod cannot be pickled. """ _fill_digi_row_base(row, event) + row['pha'] = event.pha row['column'] = event.column row['row'] = event.row - row['pha'] = event.pha row.append() class DigiDescription(tables.IsDescription): @@ -301,6 +300,7 @@ class OutputFileBase(tables.File): _DATE_FORMAT = '%a, %d %b %Y %H:%M:%S %z' _FILE_TYPE = None + _READOUT_TYPE = None def __init__(self, file_path: str) -> None: """Constructor. @@ -310,8 +310,14 @@ def __init__(self, file_path: str) -> None: self.header_group = self.create_group(self.root, 'header', 'File header') date = time.strftime(self._DATE_FORMAT) creator = pathlib.Path(inspect.stack()[-1].filename).name - self.update_header(filetype=self._FILE_TYPE.value, date=date, creator=creator, - version=__version__, tagdate=__tagdate__) + self.update_header(filetype=self._FILE_TYPE.value, readouttype=self._READOUT_TYPE.value,\ + date=date, creator=creator, version=__version__, tagdate=__tagdate__) + #if self._FILE_TYPE.value == FileType.DIGI: + # self.update_header(filetype=self._FILE_TYPE.value, readouttype=self._READOUT_TYPE.value,\ + # date=date, creator=creator, version=__version__, tagdate=__tagdate__) + # elif self._FILE_TYPE.value == FileType.RECON: + # self.update_header(filetype=self._FILE_TYPE.value,\ + # date=date, creator=creator, version=__version__, tagdate=__tagdate__) def update_header(self, **kwargs) -> None: """Update the user attributes in the header group. @@ -381,11 +387,11 @@ class DigiOutputFileSparse(OutputFileBase): """ _FILE_TYPE = FileType.DIGI - _READOUT_TYPE = HexagonalReadoutMode.SPARSE #not sure if useful + _READOUT_TYPE = HexagonalReadoutMode.SPARSE DIGI_TABLE_SPECS = ('digi_table', DigiDescriptionSparse, 'Digi data') - COLUMNS_ARRAY_SPECS = ('columns', tables.Int16Atom(shape=())) - ROWS_ARRAY_SPECS = ('rows', tables.Int16Atom(shape=())) - PHA_ARRAY_SPECS = ('pha', tables.Int16Atom(shape=())) + COLUMNS_ARRAY_SPECS = ('columns', tables.Int32Atom(shape=())) + ROWS_ARRAY_SPECS = ('rows', tables.Int32Atom(shape=())) + PHA_ARRAY_SPECS = ('pha', tables.Int32Atom(shape=())) MC_TABLE_SPECS = ('mc_table', MonteCarloDescription, 'Monte Carlo data') def __init__(self, file_path: str): @@ -542,7 +548,7 @@ def _digioutput_class(mode: HexagonalReadoutMode) -> type: """ return _FILEIO_CLASS_DICT[mode] -''' + class DigiOutputFile(OutputFileBase): """Description of a digitized output file. @@ -558,6 +564,7 @@ class DigiOutputFile(OutputFileBase): """ _FILE_TYPE = FileType.DIGI + _READOUT_TYPE = HexagonalReadoutMode.RECTANGULAR DIGI_TABLE_SPECS = ('digi_table', DigiDescription, 'Digi data') PHA_ARRAY_SPECS = ('pha', tables.Int32Atom(shape=())) MC_TABLE_SPECS = ('mc_table', MonteCarloDescription, 'Monte Carlo data') @@ -595,7 +602,7 @@ def flush(self) -> None: self.pha_array.flush() self.mc_table.flush() -''' + class ReconOutputFile(OutputFileBase): @@ -609,6 +616,7 @@ class ReconOutputFile(OutputFileBase): """ _FILE_TYPE = FileType.RECON + #_READOUT_TYPE = HexagonalReadoutMode.RECTANGULAR RECON_TABLE_SPECS = ('recon_table', ReconDescription, 'Recon data') MC_TABLE_SPECS = ('mc_table', MonteCarloDescription, 'Monte Carlo data') @@ -700,6 +708,8 @@ def __init__(self, file_path: str): """ super().__init__(file_path) self.digi_table = self.root.digi.digi_table + self.columns_array = self.root.digi.columns + self.rows_array = self.root.digi.rows self.pha_array = self.root.digi.pha self.mc_table = self.root.mc.mc_table self.__index = -1 @@ -723,8 +733,10 @@ def digi_event(self, row_index: int) -> DigiEventSparse: The index of the target row in the event file. """ row = self.digi_table[row_index] + columns = self.columns_array[row_index] + rows = self.rows_array[row_index] pha = self.pha_array[row_index] - return DigiEventSparse.from_digi(row, pha) + return DigiEventSparse.from_digi(row, pha, columns, rows) def mc_event(self, row_index: int) -> MonteCarloEvent: """Random access to the MonteCarloEvent part of the event contribution. @@ -781,7 +793,7 @@ def mc_column(self, name: str) -> np.ndarray: return self.mc_table.col(name) def digi_event(self, row_index: int) -> DigiEventRectangular: - """Random access to the DigiEventSparse part of the event contribution. + """Random access to the DigiEvent part of the event contribution. Arguments --------- @@ -832,7 +844,6 @@ def __init__(self, file_path: str): """ super().__init__(file_path) self.digi_table = self.root.digi.digi_table - self.pha_array = self.root.digi.pha self.mc_table = self.root.mc.mc_table self.__index = -1 @@ -855,8 +866,7 @@ def digi_event(self, row_index: int) -> DigiEventCircular: The index of the target row in the event file. """ row = self.digi_table[row_index] - pha = self.pha_array[row_index] - return DigiEventCircular.from_digi(row, pha) + return DigiEventCircular.from_digi(row) def mc_event(self, row_index: int) -> MonteCarloEvent: """Random access to the MonteCarloEvent part of the event contribution. @@ -888,6 +898,8 @@ def __next__(self) -> DigiEventCircular: class DigiInputFile(InputFileBase): """Description of a digitized input file. + NOTE: this class should be eliminated when the above three classes have been + fully implemented and tested. This has a very simple interface: we cache references to the relevant tables when we open the file and we provide methods to reassemble a specific table @@ -992,9 +1004,24 @@ def peek_file_type(file_path: str) -> FileType: return FileType(input_file.root.header._v_attrs['filetype']) except KeyError as exception: raise RuntimeError(f'File {file_path} has no type information.') from exception + +def peek_readout_type(file_path: str) -> HexagonalReadoutMode: + """Peek into the header of a HDF5 Digi file and determing the readout type. + + Arguments + --------- + file_path : str + The path to the input file. + """ + with tables.open_file(file_path, 'r') as input_file: + try: + return HexagonalReadoutMode(input_file.root.header._v_attrs['readouttype']) + except KeyError as exception: + raise RuntimeError(f'File {file_path} has no readout information.') from exception + def open_input_file(file_path: str) -> InputFileBase: - """Open an input file automatically determining the file type. + """Open an input file automatically determining the file type and readout type. Arguments --------- @@ -1002,8 +1029,15 @@ def open_input_file(file_path: str) -> InputFileBase: The path to the output file. """ file_type = peek_file_type(file_path) + readout_type = peek_readout_type(file_path) if file_type == FileType.DIGI: - return DigiInputFile(file_path) + if readout_type == HexagonalReadoutMode.SPARSE: + return DigiInputFileSparse(file_path) + elif readout_type == HexagonalReadoutMode.RECTANGULAR: + return DigiInputFileRectangular(file_path) + elif readout_type == HexagonalReadoutMode.CIRCULAR: + return DigiInputFileCircular(file_path) + #return DigiInputFile(file_path) if file_type == FileType.RECON: return ReconInputFile(file_path) - raise RuntimeError(f'Invalid input file type {file_type}') + raise RuntimeError(f'Invalid input file type {file_type} or invalid readout type for file type {readout_type}') diff --git a/tests/test_fileio.py b/tests/test_fileio.py index a7bdde9..ee61fd0 100644 --- a/tests/test_fileio.py +++ b/tests/test_fileio.py @@ -21,9 +21,11 @@ from hexsample import HEXSAMPLE_DATA from hexsample.digi import DigiEventSparse, DigiEventRectangular, DigiEventCircular -from hexsample.readout import HexagonalReadoutMode -from hexsample.fileio import DigiInputFile, ReconInputFile, ReconOutputFile,\ - FileType, peek_file_type, open_input_file, _FILEIO_CLASS_DICT, _digioutput_class +from hexsample.readout import HexagonalReadoutMode, HexagonalReadoutCircular +from hexsample.fileio import DigiInputFile, DigiInputFileSparse, DigiInputFileRectangular,\ + DigiInputFileCircular, ReconInputFile, ReconOutputFile,FileType,\ + DigiOutputFile, DigiOutputFileSparse, DigiOutputFileRectangular, DigiOutputFileCircular,\ + peek_file_type, peek_readout_type, open_input_file, _digioutput_class from hexsample.mc import MonteCarloEvent from hexsample.roi import RegionOfInterest, Padding @@ -42,31 +44,50 @@ def _digi_event_rectangular(index : int) -> DigiEventRectangular: return DigiEventRectangular(index, index, index, 0, pha, roi) def _digi_event_sparse(index : int) -> DigiEventSparse: - columns = np.randint(0, high=10, size=10) - rows = np.randint(0, high=10, size=10) - pha = np.full(10, index) - return DigiEventSparse(columns, rows, pha) + columns = np.array([0, 3, 6, 7], dtype=int) + rows = np.array([0, 1, 2, 4], dtype=int) + assert len(columns) == len(rows) + pha = np.full(len(columns), index) + return DigiEventSparse(index, index, index, 0, pha, columns, rows) -def _test_write_sparse(file_path, num_events : int = 10): - """Small test writing a bunch of toy event strcutures to file. - 10 events for every readout methods are created. +def _digi_event_circular(index : int) -> DigiEventCircular: + column = index + row = index + pha = np.array([100, 200, 0, 75, 43, 31, 98]) + assert len(pha) == HexagonalReadoutCircular.NUM_PIXELS + return DigiEventCircular(index, index, index, 0, pha, column, row) + +def _test_write(file_path, num_events : int = 10): + """Small test writing a bunch of toy event structures to file. + This test considers a rectangular type of readout. """ output_file = _digioutput_class(HexagonalReadoutMode.RECTANGULAR)(file_path) for i in range(num_events): output_file.add_row(_digi_event_rectangular(i), _mc_event(i)) output_file.close() -def _test_write(file_path, num_events : int = 10): - """Small test writing a bunch of toy event strcutures to file. - 10 events for every readout methods are created. +def _test_write_sparse(file_path, num_events : int = 10): + """Small test writing a bunch of toy event structures to file. + This test considers a sparse type of readout. """ - output_file = _digioutput_class(HexagonalReadoutMode.RECTANGULAR)(file_path) + output_file = _digioutput_class(HexagonalReadoutMode.SPARSE)(file_path) + #logger.info(f'Output readout type: {output_file.root.header._v_attrs['readouttype']}') for i in range(num_events): - output_file.add_row(_digi_event_rectangular(i), _mc_event(i)) + output_file.add_row(_digi_event_sparse(i), _mc_event(i)) + output_file.close() + +def _test_write_circular(file_path, num_events : int = 10): + """Small test writing a bunch of toy event structures to file. + This test considers a sparse type of readout. + """ + output_file = _digioutput_class(HexagonalReadoutMode.CIRCULAR)(file_path) + #logger.info(f'Output readout type: {output_file.root.header._v_attrs['readouttype']}') + for i in range(num_events): + output_file.add_row(_digi_event_circular(i), _mc_event(i)) output_file.close() def _test_read(file_path): - """Small test interating over an input file. + """Small test iterating over an input file. """ input_file = DigiInputFile(file_path) print(input_file) @@ -81,6 +102,40 @@ def _test_read(file_path): assert (event.pha == target.pha).all() input_file.close() +def _test_read_sparse(file_path): + """Small test iterating over a sparse input file. + """ + input_file = DigiInputFileSparse(file_path) + print(input_file) + for i, event in enumerate(input_file): + print(event.ascii()) + print(input_file.mc_event(i)) + target = _digi_event_sparse(i) + assert event.trigger_id == target.trigger_id + assert event.seconds == target.seconds + assert event.microseconds == target.microseconds + assert (event.columns == target.columns).all() + assert (event.rows == target.rows).all() + assert (event.pha == target.pha).all() + input_file.close() + +def _test_read_circular(file_path): + """Small test iterating over a sparse input file. + """ + input_file = DigiInputFileCircular(file_path) + print(input_file) + for i, event in enumerate(input_file): + print(event.ascii()) + print(input_file.mc_event(i)) + target = _digi_event_circular(i) + assert event.trigger_id == target.trigger_id + assert event.seconds == target.seconds + assert event.microseconds == target.microseconds + assert event.column == target.column + assert event.row == target.row + assert (event.pha == target.pha).all() + input_file.close() + def test(): """Write and read back a simple digi file. """ @@ -90,19 +145,50 @@ def test(): logger.info(f'Testing input file {file_path}...') _test_read(file_path) +def test_sparse(): + """Write and read back a sparse digi file. + """ + file_path = HEXSAMPLE_DATA / 'test_io.h5' + logger.info(f'Testing output file {file_path}...') + _test_write_sparse(file_path) + logger.info(f'Testing input file {file_path}...') + _test_read_sparse(file_path) + +def test_circular(): + """Write and read back a sparse digi file. + """ + file_path = HEXSAMPLE_DATA / 'test_io.h5' + logger.info(f'Testing output file {file_path}...') + _test_write_circular(file_path) + logger.info(f'Testing input file {file_path}...') + _test_read_circular(file_path) + def test_file_type(): """Test the auto-recognition machinery for input file types. """ # Test for the digi files. file_path = HEXSAMPLE_DATA / 'test_digi_filetype.h5' - #digi_file = DigiOutputFile(file_path) - #digi_file.close() + digi_file = DigiOutputFile(file_path) + digi_file.close() digi_file = DigiInputFile(file_path) assert digi_file.file_type == FileType.DIGI digi_file.close() assert peek_file_type(file_path) == FileType.DIGI digi_file = open_input_file(file_path) - assert isinstance(digi_file, DigiInputFile) + #assert isinstance(digi_file, DigiInputFile) + assert digi_file.file_type == FileType.DIGI + digi_file.close() + # Test for the digi sparse files. + file_path = HEXSAMPLE_DATA / 'test_digi_filetype_sparse_readouttype.h5' + digi_file = DigiOutputFileSparse(file_path) + digi_file.close() + digi_file = DigiInputFileSparse(file_path) + assert digi_file.file_type == FileType.DIGI + digi_file.close() + assert peek_file_type(file_path) == FileType.DIGI + assert peek_readout_type(file_path) == HexagonalReadoutMode.SPARSE + digi_file = open_input_file(file_path) + assert isinstance(digi_file, DigiInputFileSparse) assert digi_file.file_type == FileType.DIGI digi_file.close() # Test for the recon files. From ed9518c0ec076ff8f95a447f92c7a05c9938e31a Mon Sep 17 00:00:00 2001 From: chiara Date: Wed, 24 Apr 2024 09:52:34 +0200 Subject: [PATCH 31/52] Finished implementing I/O for every radout mode and tested the reading and writing --- hexsample/fileio.py | 29 +++++++++++++---------------- tests/test_fileio.py | 15 ++++++--------- 2 files changed, 19 insertions(+), 25 deletions(-) diff --git a/hexsample/fileio.py b/hexsample/fileio.py index f995876..e49cf40 100644 --- a/hexsample/fileio.py +++ b/hexsample/fileio.py @@ -300,7 +300,6 @@ class OutputFileBase(tables.File): _DATE_FORMAT = '%a, %d %b %Y %H:%M:%S %z' _FILE_TYPE = None - _READOUT_TYPE = None def __init__(self, file_path: str) -> None: """Constructor. @@ -310,14 +309,9 @@ def __init__(self, file_path: str) -> None: self.header_group = self.create_group(self.root, 'header', 'File header') date = time.strftime(self._DATE_FORMAT) creator = pathlib.Path(inspect.stack()[-1].filename).name - self.update_header(filetype=self._FILE_TYPE.value, readouttype=self._READOUT_TYPE.value,\ - date=date, creator=creator, version=__version__, tagdate=__tagdate__) - #if self._FILE_TYPE.value == FileType.DIGI: - # self.update_header(filetype=self._FILE_TYPE.value, readouttype=self._READOUT_TYPE.value,\ - # date=date, creator=creator, version=__version__, tagdate=__tagdate__) - # elif self._FILE_TYPE.value == FileType.RECON: - # self.update_header(filetype=self._FILE_TYPE.value,\ - # date=date, creator=creator, version=__version__, tagdate=__tagdate__) + self.update_header(filetype=self._FILE_TYPE.value, date=date,\ + creator=creator, version=__version__, tagdate=__tagdate__) + def update_header(self, **kwargs) -> None: """Update the user attributes in the header group. @@ -387,7 +381,7 @@ class DigiOutputFileSparse(OutputFileBase): """ _FILE_TYPE = FileType.DIGI - _READOUT_TYPE = HexagonalReadoutMode.SPARSE + _READOUT_MODE = HexagonalReadoutMode.SPARSE DIGI_TABLE_SPECS = ('digi_table', DigiDescriptionSparse, 'Digi data') COLUMNS_ARRAY_SPECS = ('columns', tables.Int32Atom(shape=())) ROWS_ARRAY_SPECS = ('rows', tables.Int32Atom(shape=())) @@ -398,6 +392,7 @@ def __init__(self, file_path: str): """Constructor. """ super().__init__(file_path) + self.update_header(readoutmode=self._READOUT_MODE.value) self.digi_group = self.create_group(self.root, 'digi', 'Digi') self.digi_table = self.create_table(self.digi_group, *self.DIGI_TABLE_SPECS) self.columns_array = self.create_vlarray(self.digi_group, *self.COLUMNS_ARRAY_SPECS) @@ -448,7 +443,7 @@ class DigiOutputFileRectangular(OutputFileBase): """ _FILE_TYPE = FileType.DIGI - _READOUT_TYPE = HexagonalReadoutMode.RECTANGULAR #not sure if useful + _READOUT_MODE = HexagonalReadoutMode.RECTANGULAR #not sure if useful DIGI_TABLE_SPECS = ('digi_table', DigiDescriptionRectangular, 'Digi data') PHA_ARRAY_SPECS = ('pha', tables.Int32Atom(shape=())) MC_TABLE_SPECS = ('mc_table', MonteCarloDescription, 'Monte Carlo data') @@ -457,6 +452,7 @@ def __init__(self, file_path: str): """Constructor. """ super().__init__(file_path) + self.update_header(readoutmode=self._READOUT_MODE.value) self.digi_group = self.create_group(self.root, 'digi', 'Digi') self.digi_table = self.create_table(self.digi_group, *self.DIGI_TABLE_SPECS) self.pha_array = self.create_vlarray(self.digi_group, *self.PHA_ARRAY_SPECS) @@ -501,7 +497,7 @@ class DigiOutputFileCircular(OutputFileBase): """ _FILE_TYPE = FileType.DIGI - _READOUT_TYPE = HexagonalReadoutMode.CIRCULAR #not sure if useful + _READOUT_MODE = HexagonalReadoutMode.CIRCULAR #not sure if useful DIGI_TABLE_SPECS = ('digi_table', DigiDescriptionCircular, 'Digi data') MC_TABLE_SPECS = ('mc_table', MonteCarloDescription, 'Monte Carlo data') @@ -509,6 +505,7 @@ def __init__(self, file_path: str): """Constructor. """ super().__init__(file_path) + self.update_header(readoutmode=self._READOUT_MODE.value) self.digi_group = self.create_group(self.root, 'digi', 'Digi') self.digi_table = self.create_table(self.digi_group, *self.DIGI_TABLE_SPECS) self.mc_group = self.create_group(self.root, 'mc', 'Monte Carlo') @@ -564,7 +561,7 @@ class DigiOutputFile(OutputFileBase): """ _FILE_TYPE = FileType.DIGI - _READOUT_TYPE = HexagonalReadoutMode.RECTANGULAR + _READOUT_MODE = HexagonalReadoutMode.RECTANGULAR DIGI_TABLE_SPECS = ('digi_table', DigiDescription, 'Digi data') PHA_ARRAY_SPECS = ('pha', tables.Int32Atom(shape=())) MC_TABLE_SPECS = ('mc_table', MonteCarloDescription, 'Monte Carlo data') @@ -573,6 +570,7 @@ def __init__(self, file_path: str): """Constructor. """ super().__init__(file_path) + self.update_header(readoutmode=self._READOUT_MODE.value) self.digi_group = self.create_group(self.root, 'digi', 'Digi') self.digi_table = self.create_table(self.digi_group, *self.DIGI_TABLE_SPECS) self.pha_array = self.create_vlarray(self.digi_group, *self.PHA_ARRAY_SPECS) @@ -616,7 +614,6 @@ class ReconOutputFile(OutputFileBase): """ _FILE_TYPE = FileType.RECON - #_READOUT_TYPE = HexagonalReadoutMode.RECTANGULAR RECON_TABLE_SPECS = ('recon_table', ReconDescription, 'Recon data') MC_TABLE_SPECS = ('mc_table', MonteCarloDescription, 'Monte Carlo data') @@ -1015,7 +1012,7 @@ def peek_readout_type(file_path: str) -> HexagonalReadoutMode: """ with tables.open_file(file_path, 'r') as input_file: try: - return HexagonalReadoutMode(input_file.root.header._v_attrs['readouttype']) + return HexagonalReadoutMode(input_file.root.header._v_attrs['readoutmode']) except KeyError as exception: raise RuntimeError(f'File {file_path} has no readout information.') from exception @@ -1029,8 +1026,8 @@ def open_input_file(file_path: str) -> InputFileBase: The path to the output file. """ file_type = peek_file_type(file_path) - readout_type = peek_readout_type(file_path) if file_type == FileType.DIGI: + readout_type = peek_readout_type(file_path) if readout_type == HexagonalReadoutMode.SPARSE: return DigiInputFileSparse(file_path) elif readout_type == HexagonalReadoutMode.RECTANGULAR: diff --git a/tests/test_fileio.py b/tests/test_fileio.py index ee61fd0..e653e57 100644 --- a/tests/test_fileio.py +++ b/tests/test_fileio.py @@ -71,7 +71,6 @@ def _test_write_sparse(file_path, num_events : int = 10): This test considers a sparse type of readout. """ output_file = _digioutput_class(HexagonalReadoutMode.SPARSE)(file_path) - #logger.info(f'Output readout type: {output_file.root.header._v_attrs['readouttype']}') for i in range(num_events): output_file.add_row(_digi_event_sparse(i), _mc_event(i)) output_file.close() @@ -81,7 +80,6 @@ def _test_write_circular(file_path, num_events : int = 10): This test considers a sparse type of readout. """ output_file = _digioutput_class(HexagonalReadoutMode.CIRCULAR)(file_path) - #logger.info(f'Output readout type: {output_file.root.header._v_attrs['readouttype']}') for i in range(num_events): output_file.add_row(_digi_event_circular(i), _mc_event(i)) output_file.close() @@ -139,7 +137,7 @@ def _test_read_circular(file_path): def test(): """Write and read back a simple digi file. """ - file_path = HEXSAMPLE_DATA / 'test_io.h5' + file_path = HEXSAMPLE_DATA / 'test_io_rectangular.h5' logger.info(f'Testing output file {file_path}...') _test_write(file_path) logger.info(f'Testing input file {file_path}...') @@ -148,7 +146,7 @@ def test(): def test_sparse(): """Write and read back a sparse digi file. """ - file_path = HEXSAMPLE_DATA / 'test_io.h5' + file_path = HEXSAMPLE_DATA / 'test_io_sparse.h5' logger.info(f'Testing output file {file_path}...') _test_write_sparse(file_path) logger.info(f'Testing input file {file_path}...') @@ -157,7 +155,7 @@ def test_sparse(): def test_circular(): """Write and read back a sparse digi file. """ - file_path = HEXSAMPLE_DATA / 'test_io.h5' + file_path = HEXSAMPLE_DATA / 'test_io_circular.h5' logger.info(f'Testing output file {file_path}...') _test_write_circular(file_path) logger.info(f'Testing input file {file_path}...') @@ -167,15 +165,14 @@ def test_file_type(): """Test the auto-recognition machinery for input file types. """ # Test for the digi files. - file_path = HEXSAMPLE_DATA / 'test_digi_filetype.h5' - digi_file = DigiOutputFile(file_path) + file_path = HEXSAMPLE_DATA / 'test_digi_filetype_rectangular.h5' + digi_file = DigiOutputFileRectangular(file_path) digi_file.close() - digi_file = DigiInputFile(file_path) + digi_file = DigiInputFileRectangular(file_path) assert digi_file.file_type == FileType.DIGI digi_file.close() assert peek_file_type(file_path) == FileType.DIGI digi_file = open_input_file(file_path) - #assert isinstance(digi_file, DigiInputFile) assert digi_file.file_type == FileType.DIGI digi_file.close() # Test for the digi sparse files. From 4c533fff547019f260243168da8b4fa0abf3902b Mon Sep 17 00:00:00 2001 From: chiara Date: Wed, 24 Apr 2024 09:58:29 +0200 Subject: [PATCH 32/52] Minor --- hexsample/fileio.py | 1 - 1 file changed, 1 deletion(-) diff --git a/hexsample/fileio.py b/hexsample/fileio.py index e49cf40..4f424b4 100644 --- a/hexsample/fileio.py +++ b/hexsample/fileio.py @@ -1034,7 +1034,6 @@ def open_input_file(file_path: str) -> InputFileBase: return DigiInputFileRectangular(file_path) elif readout_type == HexagonalReadoutMode.CIRCULAR: return DigiInputFileCircular(file_path) - #return DigiInputFile(file_path) if file_type == FileType.RECON: return ReconInputFile(file_path) raise RuntimeError(f'Invalid input file type {file_type} or invalid readout type for file type {readout_type}') From 9359be42614994efdc68239d7b33b2ae88667887 Mon Sep 17 00:00:00 2001 From: chiara Date: Wed, 24 Apr 2024 16:12:42 +0200 Subject: [PATCH 33/52] Created DigiInputFileBase class --- hexsample/fileio.py | 139 ++++++++++++++++---------------------------- 1 file changed, 49 insertions(+), 90 deletions(-) diff --git a/hexsample/fileio.py b/hexsample/fileio.py index 4f424b4..814e4cc 100644 --- a/hexsample/fileio.py +++ b/hexsample/fileio.py @@ -688,9 +688,55 @@ def header_value(self, key: str, default: Any = None) -> Any: """Return the header value corresponding to a given key. """ return self.header.get(key, default) + +class DigiInputFileBase(InputFileBase): + def __init__(self, file_path: str): + """Constructor. + """ + super().__init__(file_path) + self.digi_table = self.root.digi.digi_table + self.mc_table = self.root.mc.mc_table + self.__index = -1 + + def column(self, name: str) -> np.ndarray: + """Return a given column in the digi table. + """ + return self.digi_table.col(name) + def mc_column(self, name: str) -> np.ndarray: + """Return a given column in the Monte Carlo table. + """ + return self.mc_table.col(name) + + def mc_event(self, row_index: int) -> MonteCarloEvent: + """Random access to the MonteCarloEvent part of the event contribution. -class DigiInputFileSparse(InputFileBase): + Arguments + --------- + row_index : int + The index of the target row in the event file. + """ + row = self.mc_table[row_index] + return MonteCarloEvent(*row) + + def __iter__(self): + """Overloaded method for the implementation of the iterator protocol. + """ + self.__index = -1 + return self + + def __next__(self) -> DigiEventBase: + """Overloaded method for the implementation of the iterator protocol. + """ + self.__index += 1 + if self.__index == len(self.digi_table): + raise StopIteration + return self.digi_event(self.__index) + + pass + + +class DigiInputFileSparse(DigiInputFileBase): """Description of a sparse digitized input file. @@ -704,23 +750,11 @@ def __init__(self, file_path: str): """Constructor. """ super().__init__(file_path) - self.digi_table = self.root.digi.digi_table self.columns_array = self.root.digi.columns self.rows_array = self.root.digi.rows self.pha_array = self.root.digi.pha - self.mc_table = self.root.mc.mc_table self.__index = -1 - def column(self, name: str) -> np.ndarray: - """Return a given column in the digi table. - """ - return self.digi_table.col(name) - - def mc_column(self, name: str) -> np.ndarray: - """Return a given column in the Monte Carlo table. - """ - return self.mc_table.col(name) - def digi_event(self, row_index: int) -> DigiEventSparse: """Random access to the DigiEventSparse part of the event contribution. @@ -735,23 +769,6 @@ def digi_event(self, row_index: int) -> DigiEventSparse: pha = self.pha_array[row_index] return DigiEventSparse.from_digi(row, pha, columns, rows) - def mc_event(self, row_index: int) -> MonteCarloEvent: - """Random access to the MonteCarloEvent part of the event contribution. - - Arguments - --------- - row_index : int - The index of the target row in the event file. - """ - row = self.mc_table[row_index] - return MonteCarloEvent(*row) - - def __iter__(self): - """Overloaded method for the implementation of the iterator protocol. - """ - self.__index = -1 - return self - def __next__(self) -> DigiEventSparse: """Overloaded method for the implementation of the iterator protocol. """ @@ -760,7 +777,7 @@ def __next__(self) -> DigiEventSparse: raise StopIteration return self.digi_event(self.__index) -class DigiInputFileRectangular(InputFileBase): +class DigiInputFileRectangular(DigiInputFileBase): """Description of a rectangular digitized input file. @@ -774,21 +791,9 @@ def __init__(self, file_path: str): """Constructor. """ super().__init__(file_path) - self.digi_table = self.root.digi.digi_table self.pha_array = self.root.digi.pha - self.mc_table = self.root.mc.mc_table self.__index = -1 - def column(self, name: str) -> np.ndarray: - """Return a given column in the digi table. - """ - return self.digi_table.col(name) - - def mc_column(self, name: str) -> np.ndarray: - """Return a given column in the Monte Carlo table. - """ - return self.mc_table.col(name) - def digi_event(self, row_index: int) -> DigiEventRectangular: """Random access to the DigiEvent part of the event contribution. @@ -801,23 +806,6 @@ def digi_event(self, row_index: int) -> DigiEventRectangular: pha = self.pha_array[row_index] return DigiEventRectangular.from_digi(row, pha) - def mc_event(self, row_index: int) -> MonteCarloEvent: - """Random access to the MonteCarloEvent part of the event contribution. - - Arguments - --------- - row_index : int - The index of the target row in the event file. - """ - row = self.mc_table[row_index] - return MonteCarloEvent(*row) - - def __iter__(self): - """Overloaded method for the implementation of the iterator protocol. - """ - self.__index = -1 - return self - def __next__(self) -> DigiEventRectangular: """Overloaded method for the implementation of the iterator protocol. """ @@ -826,7 +814,7 @@ def __next__(self) -> DigiEventRectangular: raise StopIteration return self.digi_event(self.__index) -class DigiInputFileCircular(InputFileBase): +class DigiInputFileCircular(DigiInputFileBase): """Description of a circular digitized input file. @@ -840,20 +828,8 @@ def __init__(self, file_path: str): """Constructor. """ super().__init__(file_path) - self.digi_table = self.root.digi.digi_table - self.mc_table = self.root.mc.mc_table self.__index = -1 - def column(self, name: str) -> np.ndarray: - """Return a given column in the digi table. - """ - return self.digi_table.col(name) - - def mc_column(self, name: str) -> np.ndarray: - """Return a given column in the Monte Carlo table. - """ - return self.mc_table.col(name) - def digi_event(self, row_index: int) -> DigiEventCircular: """Random access to the DigiEventSparse part of the event contribution. @@ -865,23 +841,6 @@ def digi_event(self, row_index: int) -> DigiEventCircular: row = self.digi_table[row_index] return DigiEventCircular.from_digi(row) - def mc_event(self, row_index: int) -> MonteCarloEvent: - """Random access to the MonteCarloEvent part of the event contribution. - - Arguments - --------- - row_index : int - The index of the target row in the event file. - """ - row = self.mc_table[row_index] - return MonteCarloEvent(*row) - - def __iter__(self): - """Overloaded method for the implementation of the iterator protocol. - """ - self.__index = -1 - return self - def __next__(self) -> DigiEventCircular: """Overloaded method for the implementation of the iterator protocol. """ From 7559e3edf9b368e3413912cbebede8e359d2503d Mon Sep 17 00:00:00 2001 From: chiara Date: Wed, 24 Apr 2024 16:13:12 +0200 Subject: [PATCH 34/52] Modified run() function for different modes of readout --- hexsample/clustering.py | 40 +++++++++++++++++++++++++++++----------- 1 file changed, 29 insertions(+), 11 deletions(-) diff --git a/hexsample/clustering.py b/hexsample/clustering.py index c1f154c..60560b0 100644 --- a/hexsample/clustering.py +++ b/hexsample/clustering.py @@ -25,8 +25,9 @@ import numpy as np -from hexsample.digi import DigiEventRectangular +from hexsample.digi import DigiEventSparse, DigiEventRectangular, DigiEventCircular from hexsample.hexagon import HexagonalGrid +from hexsample.readout import HexagonalReadoutCircular @dataclass @@ -104,23 +105,40 @@ class ClusteringNN(ClusteringBase): num_neighbors: int - def run(self, event: DigiEventRectangular) -> Cluster: + def run(self, event: DigiEventSparse | DigiEventRectangular | DigiEventCircular) -> Cluster: """Overladed method. .. warning:: The loop ever the neighbors might likely be vectorized and streamlined for speed using proper numpy array for the offset indexes. """ + if isinstance(event, DigiEventSparse): + pass + if isinstance(event, DigiEventCircular): + #If the readout is circular, we want to take all the neirest neighbors. + self.num_neighbors = HexagonalReadoutCircular.NUM_PIXELS - 1 # -1 is bc the central px is already considered + col = [event.column] + row = [event.row] + adc_channel_order = [] + # Taking the NN in logical coordinates ... + for _col, _row in self.grid.neighbors(event.column, event.row): + col.append(_col) + row.append(_row) + # ... transforming the coordinates of the NN in its corresponding ADC channel ... + adc_channel_order.append(self.grid.adc_channel(_col, row)) + # ... reordering the pha array for the correspondance (col[i], row[i]) with pha[i]. + pha = event.pha[adc_channel_order] # pylint: disable = invalid-name - seed_col, seed_row = event.highest_pixel() - col = [seed_col] - row = [seed_row] - for _col, _row in self.grid.neighbors(seed_col, seed_row): - col.append(_col) - row.append(_row) - col = np.array(col) - row = np.array(row) - pha = np.array([event(_col, _row) for _col, _row in zip(col, row)]) + if isinstance(event, DigiEventRectangular): + seed_col, seed_row = event.highest_pixel() + col = [seed_col] + row = [seed_row] + for _col, _row in self.grid.neighbors(seed_col, seed_row): + col.append(_col) + row.append(_row) + col = np.array(col) + row = np.array(row) + pha = np.array([event(_col, _row) for _col, _row in zip(col, row)]) pha = self.zero_suppress(pha) # Array indexes in order of decreasing pha---note that we use -pha to # trick argsort into sorting values in decreasing order. From 1ddfe5d031939b622802e66f521537bf97b7c513 Mon Sep 17 00:00:00 2001 From: chiara Date: Wed, 24 Apr 2024 16:13:25 +0200 Subject: [PATCH 35/52] minor --- hexsample/bin/hxdisplay.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hexsample/bin/hxdisplay.py b/hexsample/bin/hxdisplay.py index 980f77f..53adefb 100755 --- a/hexsample/bin/hxdisplay.py +++ b/hexsample/bin/hxdisplay.py @@ -24,7 +24,7 @@ from hexsample import logger from hexsample.app import ArgumentParser -from hexsample.digi import HexagonalReadoutRectangular +from hexsample.readout import HexagonalReadoutRectangular from hexsample.display import HexagonalGridDisplay from hexsample.fileio import DigiInputFile from hexsample.hexagon import HexagonalLayout From 13b592103bd1ea1fa0566d61a5e81b1d24fac238 Mon Sep 17 00:00:00 2001 From: chiara Date: Wed, 24 Apr 2024 17:37:15 +0200 Subject: [PATCH 36/52] continued updating header in I/O files for every readout mode --- hexsample/bin/hxrecon.py | 33 +++++++++++++++++++++++++++--- hexsample/fileio.py | 2 +- hexsample/recon.py | 43 +++++++++++++++++++++++++++++++++++++++- 3 files changed, 73 insertions(+), 5 deletions(-) diff --git a/hexsample/bin/hxrecon.py b/hexsample/bin/hxrecon.py index 67a7d90..dc6e9d9 100755 --- a/hexsample/bin/hxrecon.py +++ b/hexsample/bin/hxrecon.py @@ -27,8 +27,8 @@ from hexsample import logger from hexsample.app import ArgumentParser, check_required_args from hexsample.clustering import ClusteringNN -from hexsample.readout import HexagonalReadoutRectangular -from hexsample.fileio import DigiInputFile, ReconOutputFile +from hexsample.readout import HexagonalReadoutMode, HexagonalReadoutSparse, HexagonalReadoutRectangular, HexagonalReadoutCircular +from hexsample.fileio import DigiInputFileBase, DigiInputFileSparse, DigiInputFileRectangular, DigiInputFileCircular, DigiInputFile, ReconOutputFile, peek_readout_type from hexsample.hexagon import HexagonalLayout from hexsample.recon import ReconEvent @@ -52,10 +52,37 @@ def hxrecon(**kwargs): input_file_path = str(kwargs['infile']) if not input_file_path.endswith('.h5'): raise RuntimeError('Input file {input_file_path} does not look like a HDF5 file') + ''' + readout_mode = peek_readout_type(input_file_path) + # Now we can construct a set of isinstance() + if readout_mode is HexagonalReadoutMode.SPARSE: + input_file = DigiInputFileSparse(input_file_path) + header = input_file.header + args = HexagonalLayout(header['layout']), header['numcolumns'], header['numrows'],\ + header['pitch'], header['noise'], header['gain'] + readout = HexagonalReadoutSparse(*args) + logger.info(f'Readout chip: {readout}') + #input_file = DigiInputFile(input_file_path) + elif readout_mode is HexagonalReadoutMode.RECTANGULAR: + input_file = DigiInputFile(input_file_path) + header = input_file.header + args = HexagonalLayout(header['layout']), header['numcolumns'], header['numrows'],\ + header['pitch'], header['noise'], header['gain'] + readout = HexagonalReadoutRectangular(*args) + logger.info(f'Readout chip: {readout}') + elif readout_mode is HexagonalReadoutMode.CIRCULAR: + input_file = DigiInputFileCircular(input_file_path) + header = input_file.header + args = HexagonalLayout(header['layout']), header['numcolumns'], header['numrows'],\ + header['pitch'], header['noise'], header['gain'] + readout = HexagonalReadoutCircular(*args) + logger.info(f'Readout chip: {readout}') + ''' input_file = DigiInputFile(input_file_path) header = input_file.header args = HexagonalLayout(header['layout']), header['numcolumns'], header['numrows'],\ header['pitch'], header['noise'], header['gain'] + #readout = HexagonalReadoutCircular(*args) readout = HexagonalReadoutRectangular(*args) logger.info(f'Readout chip: {readout}') clustering = ClusteringNN(readout, kwargs['zsupthreshold'], kwargs['nneighbors']) @@ -66,7 +93,7 @@ def hxrecon(**kwargs): output_file.update_digi_header(**input_file.header) for i, event in tqdm(enumerate(input_file)): cluster = clustering.run(event) - args = event.trigger_id, event.timestamp(), event.livetime, event.roi.size, cluster + args = event.trigger_id, event.timestamp(), event.livetime, cluster recon_event = ReconEvent(*args) mc_event = input_file.mc_event(i) output_file.add_row(recon_event, mc_event) diff --git a/hexsample/fileio.py b/hexsample/fileio.py index 814e4cc..8f4290e 100644 --- a/hexsample/fileio.py +++ b/hexsample/fileio.py @@ -262,7 +262,7 @@ def _fill_recon_row(row: tables.tableextension.Row, event: ReconEvent) -> None: row['trigger_id'] = event.trigger_id row['timestamp'] = event.timestamp row['livetime'] = event.livetime - row['roi_size'] = event.roi_size + #row['roi_size'] = event.roi_size row['cluster_size'] = event.cluster.size() row['energy'] = event.energy() row['posx'], row['posy'] = event.position() diff --git a/hexsample/recon.py b/hexsample/recon.py index c88de30..972407c 100644 --- a/hexsample/recon.py +++ b/hexsample/recon.py @@ -30,6 +30,47 @@ DEFAULT_IONIZATION_POTENTIAL = xraydb.ionization_potential('Si') +@dataclass +class ReconEventBase: + + """Descriptor for a reconstructed event. + + Arguments + --------- + trigger_id : int + The trigger identifier. + + timestamp : float + The timestamp (in s) of the event. + + livetime : int + The livetime (in us) since the last event. + + cluster : Cluster + The reconstructed cluster for the event. + """ + + trigger_id: int + timestamp: float + livetime: int + cluster: Cluster + + def energy(self, ionization_potential: float = DEFAULT_IONIZATION_POTENTIAL) -> float: + """Return the energy of the event in eV. + + .. warning:: + This is currently using the ionization energy of Silicon to do the + conversion, assuming a detector gain of 1. We will need to do some + bookkeeping, here, to make this work reliably. + """ + return ionization_potential * self.cluster.pulse_height() + + def position(self) -> Tuple[float, float]: + """Return the reconstructed position of the event. + """ + return self.cluster.centroid() + + @dataclass class ReconEvent: @@ -57,7 +98,7 @@ class ReconEvent: trigger_id: int timestamp: float livetime: int - roi_size: int + #roi_size: int cluster: Cluster def energy(self, ionization_potential: float = DEFAULT_IONIZATION_POTENTIAL) -> float: From 8dedfed5763365c120f1a9d163539a672442ca42 Mon Sep 17 00:00:00 2001 From: chiara Date: Wed, 24 Apr 2024 17:46:02 +0200 Subject: [PATCH 37/52] minor --- hexsample/bin/hxsim.py | 8 ++++---- hexsample/fileio.py | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/hexsample/bin/hxsim.py b/hexsample/bin/hxsim.py index 2c726f0..b44920d 100755 --- a/hexsample/bin/hxsim.py +++ b/hexsample/bin/hxsim.py @@ -65,12 +65,12 @@ def hxsim(**kwargs): photon_list = PhotonList(source, sensor, kwargs['numevents']) readout_mode = HexagonalReadoutMode(kwargs['mode']) # Is there any nicer way to do this? See https://github.com/lucabaldini/hexsample/issues/51 - if readout_mode == HexagonalReadoutMode.SPARSE: + if readout_mode is HexagonalReadoutMode.SPARSE: readout_args = kwargs['trgthreshold'], kwargs['zsupthreshold'], kwargs['offset'] - elif readout_mode == HexagonalReadoutMode.RECTANGULAR: + elif readout_mode is HexagonalReadoutMode.RECTANGULAR: padding = Padding(*kwargs['padding']) readout_args = kwargs['trgthreshold'], padding, kwargs['zsupthreshold'], kwargs['offset'] - elif readout_mode == HexagonalReadoutMode.CIRCULAR: + elif readout_mode is HexagonalReadoutMode.CIRCULAR: readout_args = kwargs['trgthreshold'], kwargs['zsupthreshold'], kwargs['offset'] else: raise RuntimeError @@ -80,7 +80,7 @@ def hxsim(**kwargs): logger.info(f'Readout chip: {readout}') output_file_path = kwargs.get('outfile') output_file = _digioutput_class(readout_mode)(output_file_path) - output_file.update_header(**kwargs) + output_file.update_header(*args) logger.info('Starting the event loop...') for mc_event in tqdm(photon_list): x, y = mc_event.propagate(sensor.trans_diffusion_sigma) diff --git a/hexsample/fileio.py b/hexsample/fileio.py index 8f4290e..f381263 100644 --- a/hexsample/fileio.py +++ b/hexsample/fileio.py @@ -539,7 +539,7 @@ def flush(self) -> None: HexagonalReadoutMode.CIRCULAR: DigiOutputFileCircular } -def _digioutput_class(mode: HexagonalReadoutMode) -> type: +def _digioutput_class(mode: HexagonalReadoutMode) -> DigiOutputFileSparse | DigiOutputFileRectangular | DigiDescriptionCircular: """Return the proper class to be used as DigiOutputFile, depending on the readout mode of the chip. """ From e0b0f884e02a10628e169c0c80b8f45a2231af6d Mon Sep 17 00:00:00 2001 From: chiara Date: Wed, 24 Apr 2024 18:40:43 +0200 Subject: [PATCH 38/52] minor --- hexsample/bin/hxsim.py | 6 +++--- hexsample/fileio.py | 2 +- tests/test_fileio.py | 8 ++++---- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/hexsample/bin/hxsim.py b/hexsample/bin/hxsim.py index b44920d..9c97461 100755 --- a/hexsample/bin/hxsim.py +++ b/hexsample/bin/hxsim.py @@ -30,7 +30,7 @@ from hexsample import HEXSAMPLE_DATA from hexsample.app import ArgumentParser from hexsample.readout import HexagonalReadoutMode, readout_chip -from hexsample.fileio import _digioutput_class +from hexsample.fileio import DigiDescriptionSparse, DigiDescriptionRectangular, DigiDescriptionCircular, digioutput_class from hexsample.hexagon import HexagonalLayout from hexsample.mc import PhotonList from hexsample.roi import Padding @@ -79,8 +79,8 @@ def hxsim(**kwargs): readout = readout_chip(readout_mode, *args) logger.info(f'Readout chip: {readout}') output_file_path = kwargs.get('outfile') - output_file = _digioutput_class(readout_mode)(output_file_path) - output_file.update_header(*args) + output_file = digioutput_class(readout_mode)(output_file_path) + output_file.update_header(**kwargs) logger.info('Starting the event loop...') for mc_event in tqdm(photon_list): x, y = mc_event.propagate(sensor.trans_diffusion_sigma) diff --git a/hexsample/fileio.py b/hexsample/fileio.py index f381263..0bf04a9 100644 --- a/hexsample/fileio.py +++ b/hexsample/fileio.py @@ -539,7 +539,7 @@ def flush(self) -> None: HexagonalReadoutMode.CIRCULAR: DigiOutputFileCircular } -def _digioutput_class(mode: HexagonalReadoutMode) -> DigiOutputFileSparse | DigiOutputFileRectangular | DigiDescriptionCircular: +def digioutput_class(mode: HexagonalReadoutMode) -> DigiOutputFileSparse | DigiOutputFileRectangular | DigiDescriptionCircular: """Return the proper class to be used as DigiOutputFile, depending on the readout mode of the chip. """ diff --git a/tests/test_fileio.py b/tests/test_fileio.py index e653e57..4103b79 100644 --- a/tests/test_fileio.py +++ b/tests/test_fileio.py @@ -25,7 +25,7 @@ from hexsample.fileio import DigiInputFile, DigiInputFileSparse, DigiInputFileRectangular,\ DigiInputFileCircular, ReconInputFile, ReconOutputFile,FileType,\ DigiOutputFile, DigiOutputFileSparse, DigiOutputFileRectangular, DigiOutputFileCircular,\ - peek_file_type, peek_readout_type, open_input_file, _digioutput_class + peek_file_type, peek_readout_type, open_input_file, digioutput_class from hexsample.mc import MonteCarloEvent from hexsample.roi import RegionOfInterest, Padding @@ -61,7 +61,7 @@ def _test_write(file_path, num_events : int = 10): """Small test writing a bunch of toy event structures to file. This test considers a rectangular type of readout. """ - output_file = _digioutput_class(HexagonalReadoutMode.RECTANGULAR)(file_path) + output_file = digioutput_class(HexagonalReadoutMode.RECTANGULAR)(file_path) for i in range(num_events): output_file.add_row(_digi_event_rectangular(i), _mc_event(i)) output_file.close() @@ -70,7 +70,7 @@ def _test_write_sparse(file_path, num_events : int = 10): """Small test writing a bunch of toy event structures to file. This test considers a sparse type of readout. """ - output_file = _digioutput_class(HexagonalReadoutMode.SPARSE)(file_path) + output_file = digioutput_class(HexagonalReadoutMode.SPARSE)(file_path) for i in range(num_events): output_file.add_row(_digi_event_sparse(i), _mc_event(i)) output_file.close() @@ -79,7 +79,7 @@ def _test_write_circular(file_path, num_events : int = 10): """Small test writing a bunch of toy event structures to file. This test considers a sparse type of readout. """ - output_file = _digioutput_class(HexagonalReadoutMode.CIRCULAR)(file_path) + output_file = digioutput_class(HexagonalReadoutMode.CIRCULAR)(file_path) for i in range(num_events): output_file.add_row(_digi_event_circular(i), _mc_event(i)) output_file.close() From 32b55c431674a0edee2f649fe24f9ce9bae2143b Mon Sep 17 00:00:00 2001 From: chiara Date: Mon, 29 Apr 2024 11:46:58 +0200 Subject: [PATCH 39/52] minor --- hexsample/clustering.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/hexsample/clustering.py b/hexsample/clustering.py index 60560b0..f1e4f66 100644 --- a/hexsample/clustering.py +++ b/hexsample/clustering.py @@ -115,19 +115,23 @@ def run(self, event: DigiEventSparse | DigiEventRectangular | DigiEventCircular) if isinstance(event, DigiEventSparse): pass if isinstance(event, DigiEventCircular): - #If the readout is circular, we want to take all the neirest neighbors. + # If the readout is circular, we want to take all the neirest neighbors. self.num_neighbors = HexagonalReadoutCircular.NUM_PIXELS - 1 # -1 is bc the central px is already considered col = [event.column] row = [event.row] - adc_channel_order = [] + adc_channel_order = [self.grid.adc_channel(event.column, event.row)] # Taking the NN in logical coordinates ... for _col, _row in self.grid.neighbors(event.column, event.row): col.append(_col) row.append(_row) # ... transforming the coordinates of the NN in its corresponding ADC channel ... - adc_channel_order.append(self.grid.adc_channel(_col, row)) + adc_channel_order.append(self.grid.adc_channel(_col, _row)) # ... reordering the pha array for the correspondance (col[i], row[i]) with pha[i]. pha = event.pha[adc_channel_order] + # Converting lists into numpy arrays + col = np.array(col) + row = np.array(row) + pha = np.array(pha) # pylint: disable = invalid-name if isinstance(event, DigiEventRectangular): seed_col, seed_row = event.highest_pixel() @@ -139,6 +143,7 @@ def run(self, event: DigiEventSparse | DigiEventRectangular | DigiEventCircular) col = np.array(col) row = np.array(row) pha = np.array([event(_col, _row) for _col, _row in zip(col, row)]) + # Zero suppressing the event (whatever the readout type)... pha = self.zero_suppress(pha) # Array indexes in order of decreasing pha---note that we use -pha to # trick argsort into sorting values in decreasing order. From c46036f9297b7ff33a9071f52599b5a890a3cea0 Mon Sep 17 00:00:00 2001 From: chiara Date: Mon, 29 Apr 2024 12:15:16 +0200 Subject: [PATCH 40/52] /minor/ --- hexsample/bin/hxrecon.py | 61 ++++++++++++++++++++++------------------ 1 file changed, 34 insertions(+), 27 deletions(-) diff --git a/hexsample/bin/hxrecon.py b/hexsample/bin/hxrecon.py index dc6e9d9..67763c1 100755 --- a/hexsample/bin/hxrecon.py +++ b/hexsample/bin/hxrecon.py @@ -52,42 +52,49 @@ def hxrecon(**kwargs): input_file_path = str(kwargs['infile']) if not input_file_path.endswith('.h5'): raise RuntimeError('Input file {input_file_path} does not look like a HDF5 file') - ''' - readout_mode = peek_readout_type(input_file_path) - # Now we can construct a set of isinstance() - if readout_mode is HexagonalReadoutMode.SPARSE: - input_file = DigiInputFileSparse(input_file_path) - header = input_file.header - args = HexagonalLayout(header['layout']), header['numcolumns'], header['numrows'],\ - header['pitch'], header['noise'], header['gain'] - readout = HexagonalReadoutSparse(*args) - logger.info(f'Readout chip: {readout}') - #input_file = DigiInputFile(input_file_path) - elif readout_mode is HexagonalReadoutMode.RECTANGULAR: + + # It is necessary to extract the reaodut type because every readout type + # corresponds to a different DigiEvent type. + # If there is no info about the readout, we try to reconstruct with a Rectangular + # readout mode, that is the mode for all the old reconstruction before the + # distinction between different readouts was implemented. + try: + readout_mode = peek_readout_type(input_file_path) + # Now we can construct a set of if/else. + if readout_mode is HexagonalReadoutMode.SPARSE: + input_file = DigiInputFileSparse(input_file_path) + header = input_file.header + args = HexagonalLayout(header['layout']), header['numcolumns'], header['numrows'],\ + header['pitch'], header['noise'], header['gain'] + readout = HexagonalReadoutSparse(*args) + logger.info(f'Readout chip: {readout}') + elif readout_mode is HexagonalReadoutMode.RECTANGULAR: + input_file = DigiInputFile(input_file_path) + header = input_file.header + args = HexagonalLayout(header['layout']), header['numcolumns'], header['numrows'],\ + header['pitch'], header['noise'], header['gain'] + readout = HexagonalReadoutRectangular(*args) + logger.info(f'Readout chip: {readout}') + elif readout_mode is HexagonalReadoutMode.CIRCULAR: + input_file = DigiInputFileCircular(input_file_path) + header = input_file.header + args = HexagonalLayout(header['layout']), header['numcolumns'], header['numrows'],\ + header['pitch'], header['noise'], header['gain'] + readout = HexagonalReadoutCircular(*args) + logger.info(f'Readout chip: {readout}') + except RuntimeError: input_file = DigiInputFile(input_file_path) header = input_file.header args = HexagonalLayout(header['layout']), header['numcolumns'], header['numrows'],\ header['pitch'], header['noise'], header['gain'] + #readout = HexagonalReadoutCircular(*args) readout = HexagonalReadoutRectangular(*args) logger.info(f'Readout chip: {readout}') - elif readout_mode is HexagonalReadoutMode.CIRCULAR: - input_file = DigiInputFileCircular(input_file_path) - header = input_file.header - args = HexagonalLayout(header['layout']), header['numcolumns'], header['numrows'],\ - header['pitch'], header['noise'], header['gain'] - readout = HexagonalReadoutCircular(*args) - logger.info(f'Readout chip: {readout}') - ''' - input_file = DigiInputFile(input_file_path) - header = input_file.header - args = HexagonalLayout(header['layout']), header['numcolumns'], header['numrows'],\ - header['pitch'], header['noise'], header['gain'] - #readout = HexagonalReadoutCircular(*args) - readout = HexagonalReadoutRectangular(*args) - logger.info(f'Readout chip: {readout}') + # When the readout tipology is determined, the event is clustered ... clustering = ClusteringNN(readout, kwargs['zsupthreshold'], kwargs['nneighbors']) suffix = kwargs['suffix'] output_file_path = input_file_path.replace('.h5', f'_{suffix}.h5') + # ... and saved into an output file. output_file = ReconOutputFile(output_file_path) output_file.update_header(**kwargs) output_file.update_digi_header(**input_file.header) From 460f2412dfff7cca6ff333abc87775d937558706 Mon Sep 17 00:00:00 2001 From: chiara Date: Mon, 29 Apr 2024 13:47:49 +0200 Subject: [PATCH 41/52] minor --- hexsample/bin/hxsim.py | 26 ++++++++++++++++++++++++++ hexsample/clustering.py | 4 ++-- 2 files changed, 28 insertions(+), 2 deletions(-) diff --git a/hexsample/bin/hxsim.py b/hexsample/bin/hxsim.py index 9c97461..9603468 100755 --- a/hexsample/bin/hxsim.py +++ b/hexsample/bin/hxsim.py @@ -82,6 +82,32 @@ def hxsim(**kwargs): output_file = digioutput_class(readout_mode)(output_file_path) output_file.update_header(**kwargs) logger.info('Starting the event loop...') + """ + while trigger_id <= kwargs['numevents']: + photon_event = PhotonList(source, sensor, 1) + readout_mode = HexagonalReadoutMode(kwargs['mode']) + # Is there any nicer way to do this? See https://github.com/lucabaldini/hexsample/issues/51 + if readout_mode is HexagonalReadoutMode.SPARSE: + readout_args = kwargs['trgthreshold'], kwargs['zsupthreshold'], kwargs['offset'] + elif readout_mode is HexagonalReadoutMode.RECTANGULAR: + padding = Padding(*kwargs['padding']) + readout_args = kwargs['trgthreshold'], padding, kwargs['zsupthreshold'], kwargs['offset'] + elif readout_mode is HexagonalReadoutMode.CIRCULAR: + readout_args = kwargs['trgthreshold'], kwargs['zsupthreshold'], kwargs['offset'] + else: + raise RuntimeError + args = HexagonalLayout(kwargs['layout']), kwargs['numcolumns'], kwargs['numrows'],\ + kwargs['pitch'], kwargs['noise'], kwargs['gain'] + readout = readout_chip(readout_mode, *args) + logger.info(f'Readout chip: {readout}') + output_file_path = kwargs.get('outfile') + output_file = digioutput_class(readout_mode)(output_file_path) + output_file.update_header(**kwargs) + logger.info('Starting the event loop...') + x, y = mc_event.propagate(sensor.trans_diffusion_sigma) + digi_event = readout.read(mc_event.timestamp, x, y, *readout_args) + output_file.add_row(digi_event, mc_event) + """ for mc_event in tqdm(photon_list): x, y = mc_event.propagate(sensor.trans_diffusion_sigma) digi_event = readout.read(mc_event.timestamp, x, y, *readout_args) diff --git a/hexsample/clustering.py b/hexsample/clustering.py index f1e4f66..f84f974 100644 --- a/hexsample/clustering.py +++ b/hexsample/clustering.py @@ -114,7 +114,7 @@ def run(self, event: DigiEventSparse | DigiEventRectangular | DigiEventCircular) """ if isinstance(event, DigiEventSparse): pass - if isinstance(event, DigiEventCircular): + elif isinstance(event, DigiEventCircular): # If the readout is circular, we want to take all the neirest neighbors. self.num_neighbors = HexagonalReadoutCircular.NUM_PIXELS - 1 # -1 is bc the central px is already considered col = [event.column] @@ -133,7 +133,7 @@ def run(self, event: DigiEventSparse | DigiEventRectangular | DigiEventCircular) row = np.array(row) pha = np.array(pha) # pylint: disable = invalid-name - if isinstance(event, DigiEventRectangular): + elif isinstance(event, DigiEventRectangular): seed_col, seed_row = event.highest_pixel() col = [seed_col] row = [seed_row] From c4ba7ead048dd41447d4f460416c24be82627e61 Mon Sep 17 00:00:00 2001 From: chiara Date: Mon, 29 Apr 2024 13:51:46 +0200 Subject: [PATCH 42/52] minor --- hexsample/bin/hxsim.py | 26 -------------------------- 1 file changed, 26 deletions(-) diff --git a/hexsample/bin/hxsim.py b/hexsample/bin/hxsim.py index 9603468..9c97461 100755 --- a/hexsample/bin/hxsim.py +++ b/hexsample/bin/hxsim.py @@ -82,32 +82,6 @@ def hxsim(**kwargs): output_file = digioutput_class(readout_mode)(output_file_path) output_file.update_header(**kwargs) logger.info('Starting the event loop...') - """ - while trigger_id <= kwargs['numevents']: - photon_event = PhotonList(source, sensor, 1) - readout_mode = HexagonalReadoutMode(kwargs['mode']) - # Is there any nicer way to do this? See https://github.com/lucabaldini/hexsample/issues/51 - if readout_mode is HexagonalReadoutMode.SPARSE: - readout_args = kwargs['trgthreshold'], kwargs['zsupthreshold'], kwargs['offset'] - elif readout_mode is HexagonalReadoutMode.RECTANGULAR: - padding = Padding(*kwargs['padding']) - readout_args = kwargs['trgthreshold'], padding, kwargs['zsupthreshold'], kwargs['offset'] - elif readout_mode is HexagonalReadoutMode.CIRCULAR: - readout_args = kwargs['trgthreshold'], kwargs['zsupthreshold'], kwargs['offset'] - else: - raise RuntimeError - args = HexagonalLayout(kwargs['layout']), kwargs['numcolumns'], kwargs['numrows'],\ - kwargs['pitch'], kwargs['noise'], kwargs['gain'] - readout = readout_chip(readout_mode, *args) - logger.info(f'Readout chip: {readout}') - output_file_path = kwargs.get('outfile') - output_file = digioutput_class(readout_mode)(output_file_path) - output_file.update_header(**kwargs) - logger.info('Starting the event loop...') - x, y = mc_event.propagate(sensor.trans_diffusion_sigma) - digi_event = readout.read(mc_event.timestamp, x, y, *readout_args) - output_file.add_row(digi_event, mc_event) - """ for mc_event in tqdm(photon_list): x, y = mc_event.propagate(sensor.trans_diffusion_sigma) digi_event = readout.read(mc_event.timestamp, x, y, *readout_args) From 9e05588e85f6e363eb6ec9f8dad3f050a368cadc Mon Sep 17 00:00:00 2001 From: chiara Date: Mon, 29 Apr 2024 17:00:31 +0200 Subject: [PATCH 43/52] Implemented some plots for evaluating reconstruction --- hexsample/bin/hxview.py | 50 ++++++++++++++++++++++++++++------------- 1 file changed, 35 insertions(+), 15 deletions(-) diff --git a/hexsample/bin/hxview.py b/hexsample/bin/hxview.py index 8003c3e..dcf41d7 100755 --- a/hexsample/bin/hxview.py +++ b/hexsample/bin/hxview.py @@ -27,7 +27,8 @@ from hexsample.app import ArgumentParser from hexsample.fileio import ReconInputFile -from hexsample.plot import plt +from hexsample.hist import Histogram1d, Histogram2d +from hexsample.plot import plt, setup_gca from hexsample.analysis import create_histogram @@ -38,26 +39,45 @@ # Parser object. HXVIEW_ARGPARSER = ArgumentParser(description=__description__) HXVIEW_ARGPARSER.add_infile() -HXVIEW_ARGPARSER.add_argument("attribute", type=str, help='Attribute to be viewed.\ - To be taken from\ the following list:\ - trigger_id (int), timestamp (float), livetime (int) \ - roi_size (int), energy (float), position (Tuple[float,float])\ - cluster_size (int), roi_size (int)') -HXVIEW_ARGPARSER.add_argument("mc_table", type=str, help='Tells if the quantities are in mc table\ - accepts True or False.') def hxview(**kwargs): """View the file content. - Shows histograms of energy and cluster_size. + Shows histograms of energy and cluster_size of recon events vs their MC truth. """ input_file = ReconInputFile(kwargs['infile']) - attribute = kwargs['attribute'] - is_mc = literal_eval(kwargs['mc_table']) - print(is_mc) - histo = create_histogram(input_file, attribute, mc = is_mc) - plt.figure() + # Plotting the reconstructed energy and the true energy + histo = create_histogram(input_file, 'energy', mc = False) + mc_histo = create_histogram(input_file, 'energy', mc = True) + plt.figure('Photons energy') + histo.plot(label='Reconstructed') + mc_histo.plot(label='MonteCarlo') + plt.xlabel('Energy [eV]') + plt.legend() + + # Plotting the reconstructed x and y position and the true position. + plt.figure('Reconstructed photons position') + binning = np.linspace(-5. * 0.1, 5. * 0.1, 100) + x = input_file.column('posx') + y = input_file.column('posy') + histo = Histogram2d(binning, binning).fill(x, y) histo.plot() - plt.xlabel(attribute) + setup_gca(xlabel='x [cm]', ylabel='y [cm]') + plt.figure('True photons position') + x_mc = input_file.mc_column('absx') + y_mc = input_file.mc_column('absy') + histo_mc = Histogram2d(binning, binning).fill(x_mc, y_mc) + histo_mc.plot() + setup_gca(xlabel='x [cm]', ylabel='y [cm]') + #Closing the file and showing the figures. + plt.figure('x-direction resolution') + binning = np.linspace((x-x_mc).min(), (x-x_mc).max(), 100) + histx = Histogram1d(binning, xlabel=r'$x - x_{MC}$ [cm]').fill(x-x_mc) + histx.plot() + plt.figure('y-direction resolution') + binning = np.linspace((y-y_mc).min(), (y-y_mc).max(), 100) + histy = Histogram1d(binning, xlabel=r'$y - y_{MC}$ [cm]').fill(y-y_mc) + histy.plot() + input_file.close() plt.show() From 2b9eeef7e9a83018981142acc0dc9c21e4d15e73 Mon Sep 17 00:00:00 2001 From: chiara Date: Mon, 29 Apr 2024 17:04:32 +0200 Subject: [PATCH 44/52] minor --- hexsample/readout.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/hexsample/readout.py b/hexsample/readout.py index 471f5c6..2744d96 100644 --- a/hexsample/readout.py +++ b/hexsample/readout.py @@ -17,7 +17,7 @@ # with this program; if not, write to the Free Software Foundation Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. -"""Readout facilities facilities. +"""Readout facilities. """ from collections import Counter @@ -231,6 +231,7 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl # .. and digitize the pha values. pha = self.digitize(pha, zero_sup_threshold, offset) seconds, microseconds, livetime = self.latch_timestamp(timestamp) + # And do not forget to increment the trigger identifier! self.trigger_id += 1 return DigiEventSparse(self.trigger_id, seconds, microseconds, livetime, pha, columns, rows) @@ -486,7 +487,7 @@ def read(self, timestamp: float, x: np.ndarray, y: np.ndarray, trg_threshold: fl # .. and digitize the pha values. pha = self.digitize(pha, zero_sup_threshold, offset) seconds, microseconds, livetime = self.latch_timestamp(timestamp) - # Do not forget to update the trigger_id! + # And do not forget to increment the trigger identifier! self.trigger_id += 1 #The pha array is always in the order [pha(adc0), pha(adc1), pha(adc2), pha(adc3), pha(adc4), pha(adc5), pha(adc6)] return DigiEventCircular(self.trigger_id, seconds, microseconds, livetime, pha, *coord_max) From c0be1ef6d9e576d5927bce7ce39350d5cd4712f0 Mon Sep 17 00:00:00 2001 From: chiara Date: Mon, 29 Apr 2024 17:08:06 +0200 Subject: [PATCH 45/52] minor --- hexsample/bin/hxview.py | 52 +++++++++++++---------------------------- 1 file changed, 16 insertions(+), 36 deletions(-) diff --git a/hexsample/bin/hxview.py b/hexsample/bin/hxview.py index dcf41d7..d01819c 100755 --- a/hexsample/bin/hxview.py +++ b/hexsample/bin/hxview.py @@ -27,8 +27,7 @@ from hexsample.app import ArgumentParser from hexsample.fileio import ReconInputFile -from hexsample.hist import Histogram1d, Histogram2d -from hexsample.plot import plt, setup_gca +from hexsample.plot import plt from hexsample.analysis import create_histogram @@ -39,47 +38,28 @@ # Parser object. HXVIEW_ARGPARSER = ArgumentParser(description=__description__) HXVIEW_ARGPARSER.add_infile() +HXVIEW_ARGPARSER.add_argument("attribute", type=str, help='Attribute to be viewed.\ + To be taken from\ the following list:\ + trigger_id (int), timestamp (float), livetime (int) \ + roi_size (int), energy (float), position (Tuple[float,float])\ + cluster_size (int), roi_size (int)') +HXVIEW_ARGPARSER.add_argument("mc_table", type=str, help='Tells if the quantities are in mc table\ + accepts True or False.') def hxview(**kwargs): """View the file content. - Shows histograms of energy and cluster_size of recon events vs their MC truth. + Shows histograms of energy and cluster_size. """ input_file = ReconInputFile(kwargs['infile']) - # Plotting the reconstructed energy and the true energy - histo = create_histogram(input_file, 'energy', mc = False) - mc_histo = create_histogram(input_file, 'energy', mc = True) - plt.figure('Photons energy') - histo.plot(label='Reconstructed') - mc_histo.plot(label='MonteCarlo') - plt.xlabel('Energy [eV]') - plt.legend() - - # Plotting the reconstructed x and y position and the true position. - plt.figure('Reconstructed photons position') - binning = np.linspace(-5. * 0.1, 5. * 0.1, 100) - x = input_file.column('posx') - y = input_file.column('posy') - histo = Histogram2d(binning, binning).fill(x, y) + attribute = kwargs['attribute'] + is_mc = literal_eval(kwargs['mc_table']) + print(is_mc) + histo = create_histogram(input_file, attribute, mc = is_mc) + plt.figure() histo.plot() - setup_gca(xlabel='x [cm]', ylabel='y [cm]') - plt.figure('True photons position') - x_mc = input_file.mc_column('absx') - y_mc = input_file.mc_column('absy') - histo_mc = Histogram2d(binning, binning).fill(x_mc, y_mc) - histo_mc.plot() - setup_gca(xlabel='x [cm]', ylabel='y [cm]') - #Closing the file and showing the figures. - plt.figure('x-direction resolution') - binning = np.linspace((x-x_mc).min(), (x-x_mc).max(), 100) - histx = Histogram1d(binning, xlabel=r'$x - x_{MC}$ [cm]').fill(x-x_mc) - histx.plot() - plt.figure('y-direction resolution') - binning = np.linspace((y-y_mc).min(), (y-y_mc).max(), 100) - histy = Histogram1d(binning, xlabel=r'$y - y_{MC}$ [cm]').fill(y-y_mc) - histy.plot() - + plt.xlabel(attribute) input_file.close() plt.show() if __name__ == '__main__': - hxview(**vars(HXVIEW_ARGPARSER.parse_args())) + hxview(**vars(HXVIEW_ARGPARSER.parse_args())) \ No newline at end of file From 255fa85b9bb3a444058adb6900e1cb70cde6d9fe Mon Sep 17 00:00:00 2001 From: chiara Date: Mon, 29 Apr 2024 17:08:36 +0200 Subject: [PATCH 46/52] created file for recon vs MC quantities evaluation --- hexsample/bin/hxplots.py | 85 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 85 insertions(+) create mode 100644 hexsample/bin/hxplots.py diff --git a/hexsample/bin/hxplots.py b/hexsample/bin/hxplots.py new file mode 100644 index 0000000..9aee675 --- /dev/null +++ b/hexsample/bin/hxplots.py @@ -0,0 +1,85 @@ +#!/usr/bin/env python +# +# Copyright (C) 2023 luca.baldini@pi.infn.it +# +# For the license terms see the file LICENSE, distributed along with this +# software. +# +# This program is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 2 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# General Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +"""Event file viewer comparing reconstructed quantities with MC truth. +""" +from ast import literal_eval + +import numpy as np + +from hexsample.app import ArgumentParser +from hexsample.fileio import ReconInputFile +from hexsample.hist import Histogram1d, Histogram2d +from hexsample.plot import plt, setup_gca +from hexsample.analysis import create_histogram + + +__description__ = \ +"""Simple viewer for reconstructed event lists. +""" + +# Parser object. +HXVIEW_ARGPARSER = ArgumentParser(description=__description__) +HXVIEW_ARGPARSER.add_infile() + +def hxview(**kwargs): + """View the file content. + Shows histograms of energy and cluster_size of recon events vs their MC truth. + """ + input_file = ReconInputFile(kwargs['infile']) + # Plotting the reconstructed energy and the true energy + histo = create_histogram(input_file, 'energy', mc = False) + mc_histo = create_histogram(input_file, 'energy', mc = True) + plt.figure('Photons energy') + histo.plot(label='Reconstructed') + mc_histo.plot(label='MonteCarlo') + plt.xlabel('Energy [eV]') + plt.legend() + + # Plotting the reconstructed x and y position and the true position. + plt.figure('Reconstructed photons position') + binning = np.linspace(-5. * 0.1, 5. * 0.1, 100) + x = input_file.column('posx') + y = input_file.column('posy') + histo = Histogram2d(binning, binning).fill(x, y) + histo.plot() + setup_gca(xlabel='x [cm]', ylabel='y [cm]') + plt.figure('True photons position') + x_mc = input_file.mc_column('absx') + y_mc = input_file.mc_column('absy') + histo_mc = Histogram2d(binning, binning).fill(x_mc, y_mc) + histo_mc.plot() + setup_gca(xlabel='x [cm]', ylabel='y [cm]') + #Closing the file and showing the figures. + plt.figure('x-direction resolution') + binning = np.linspace((x-x_mc).min(), (x-x_mc).max(), 100) + histx = Histogram1d(binning, xlabel=r'$x - x_{MC}$ [cm]').fill(x-x_mc) + histx.plot() + plt.figure('y-direction resolution') + binning = np.linspace((y-y_mc).min(), (y-y_mc).max(), 100) + histy = Histogram1d(binning, xlabel=r'$y - y_{MC}$ [cm]').fill(y-y_mc) + histy.plot() + + input_file.close() + plt.show() + +if __name__ == '__main__': + hxview(**vars(HXVIEW_ARGPARSER.parse_args())) From f0d4f2668099f869685f2641e97c56ed8e3ca9f3 Mon Sep 17 00:00:00 2001 From: chiara Date: Mon, 29 Apr 2024 17:09:37 +0200 Subject: [PATCH 47/52] minor --- hexsample/bin/hxplots.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/hexsample/bin/hxplots.py b/hexsample/bin/hxplots.py index 9aee675..356324e 100644 --- a/hexsample/bin/hxplots.py +++ b/hexsample/bin/hxplots.py @@ -33,7 +33,8 @@ __description__ = \ -"""Simple viewer for reconstructed event lists. +"""Simple viewer for comparing reconstructed energy and position with the MC +true values. """ # Parser object. From 2854e3941fe51e25ae34e6cc0aa841d3c983c6e3 Mon Sep 17 00:00:00 2001 From: Luca Baldini Date: Tue, 30 Apr 2024 08:22:40 +0200 Subject: [PATCH 48/52] Minor. --- hexsample/app.py | 2 +- hexsample/bin/hxsim.py | 2 +- hexsample/fileio.py | 42 +++++++++++++++++++++--------------------- 3 files changed, 23 insertions(+), 23 deletions(-) diff --git a/hexsample/app.py b/hexsample/app.py index a40f2d8..c2a4c41 100644 --- a/hexsample/app.py +++ b/hexsample/app.py @@ -156,7 +156,7 @@ def add_readout_options(self) -> None: group.add_argument('--pitch', type=float, default=0.005, help='pitch of the readout chip in cm') modes = [item.value for item in HexagonalReadoutMode] - group.add_argument('--mode', type=str, choices=modes, default='RECTANGULAR', + group.add_argument('--readoutmode', type=str, choices=modes, default='RECTANGULAR', help='readout mode') group.add_argument('--padding', type=int, nargs=4, default=(2, 2, 2, 2), help='padding on the four sides of the ROT') diff --git a/hexsample/bin/hxsim.py b/hexsample/bin/hxsim.py index 9c97461..f37ac50 100755 --- a/hexsample/bin/hxsim.py +++ b/hexsample/bin/hxsim.py @@ -63,7 +63,7 @@ def hxsim(**kwargs): material = Material(kwargs['actmedium'], kwargs['fano']) sensor = Sensor(material, kwargs['thickness'], kwargs['transdiffsigma']) photon_list = PhotonList(source, sensor, kwargs['numevents']) - readout_mode = HexagonalReadoutMode(kwargs['mode']) + readout_mode = HexagonalReadoutMode(kwargs['readoutmode']) # Is there any nicer way to do this? See https://github.com/lucabaldini/hexsample/issues/51 if readout_mode is HexagonalReadoutMode.SPARSE: readout_args = kwargs['trgthreshold'], kwargs['zsupthreshold'], kwargs['offset'] diff --git a/hexsample/fileio.py b/hexsample/fileio.py index 0bf04a9..c8ccb91 100644 --- a/hexsample/fileio.py +++ b/hexsample/fileio.py @@ -73,7 +73,7 @@ class DigiDescriptionBase(tables.IsDescription): """Base class for the description of the (flat) digi part of the file format. It contains the trigger_id and time coordinates of the event, common to all - readout types. + readout types. """ trigger_id = tables.Int32Col(pos=0) @@ -140,7 +140,7 @@ class DigiDescriptionRectangular(DigiDescriptionBase): def _fill_digi_row_rectangular(row: tables.tableextension.Row, event: DigiEventBase) -> None: """Overloaded method. It uses the _fill_digi_row_base() function for filling the trigger_id and time - coordinates of the event. + coordinates of the event. .. note:: This would have naturally belonged to the DigiDescriptionRectangular class as @@ -172,7 +172,7 @@ class DigiDescriptionCircular(DigiDescriptionBase): def _fill_digi_row_circular(row: tables.tableextension.Row, event: DigiEventBase) -> None: """Overloaded method. It uses the _fill_digi_row_base() function for filling the trigger_id and time - coordinates of the event. + coordinates of the event. .. note:: This would have naturally belonged to the DigiDescriptionCircular class as @@ -189,7 +189,7 @@ def _fill_digi_row_circular(row: tables.tableextension.Row, event: DigiEventBase class DigiDescription(tables.IsDescription): """Description of the (flat) digi part of the file format. - NOTE: This should be eliminated when the above three classes will be fully + NOTE: This should be eliminated when the above three classes will be fully implemented and tested. """ @@ -365,7 +365,7 @@ def flush(self) -> None: This needs to be reimplemented in derived classes. """ raise NotImplementedError - + class DigiOutputFileSparse(OutputFileBase): """Description of a sparse digitized output file. @@ -381,7 +381,7 @@ class DigiOutputFileSparse(OutputFileBase): """ _FILE_TYPE = FileType.DIGI - _READOUT_MODE = HexagonalReadoutMode.SPARSE + #_READOUT_MODE = HexagonalReadoutMode.SPARSE DIGI_TABLE_SPECS = ('digi_table', DigiDescriptionSparse, 'Digi data') COLUMNS_ARRAY_SPECS = ('columns', tables.Int32Atom(shape=())) ROWS_ARRAY_SPECS = ('rows', tables.Int32Atom(shape=())) @@ -392,7 +392,7 @@ def __init__(self, file_path: str): """Constructor. """ super().__init__(file_path) - self.update_header(readoutmode=self._READOUT_MODE.value) + #self.update_header(readoutmode=self._READOUT_MODE.value) self.digi_group = self.create_group(self.root, 'digi', 'Digi') self.digi_table = self.create_table(self.digi_group, *self.DIGI_TABLE_SPECS) self.columns_array = self.create_vlarray(self.digi_group, *self.COLUMNS_ARRAY_SPECS) @@ -443,7 +443,7 @@ class DigiOutputFileRectangular(OutputFileBase): """ _FILE_TYPE = FileType.DIGI - _READOUT_MODE = HexagonalReadoutMode.RECTANGULAR #not sure if useful + #_READOUT_MODE = HexagonalReadoutMode.RECTANGULAR #not sure if useful DIGI_TABLE_SPECS = ('digi_table', DigiDescriptionRectangular, 'Digi data') PHA_ARRAY_SPECS = ('pha', tables.Int32Atom(shape=())) MC_TABLE_SPECS = ('mc_table', MonteCarloDescription, 'Monte Carlo data') @@ -452,7 +452,7 @@ def __init__(self, file_path: str): """Constructor. """ super().__init__(file_path) - self.update_header(readoutmode=self._READOUT_MODE.value) + #self.update_header(readoutmode=self._READOUT_MODE.value) self.digi_group = self.create_group(self.root, 'digi', 'Digi') self.digi_table = self.create_table(self.digi_group, *self.DIGI_TABLE_SPECS) self.pha_array = self.create_vlarray(self.digi_group, *self.PHA_ARRAY_SPECS) @@ -497,7 +497,7 @@ class DigiOutputFileCircular(OutputFileBase): """ _FILE_TYPE = FileType.DIGI - _READOUT_MODE = HexagonalReadoutMode.CIRCULAR #not sure if useful + #_READOUT_MODE = HexagonalReadoutMode.CIRCULAR #not sure if useful DIGI_TABLE_SPECS = ('digi_table', DigiDescriptionCircular, 'Digi data') MC_TABLE_SPECS = ('mc_table', MonteCarloDescription, 'Monte Carlo data') @@ -505,7 +505,7 @@ def __init__(self, file_path: str): """Constructor. """ super().__init__(file_path) - self.update_header(readoutmode=self._READOUT_MODE.value) + #self.update_header(readoutmode=self._READOUT_MODE.value) self.digi_group = self.create_group(self.root, 'digi', 'Digi') self.digi_table = self.create_table(self.digi_group, *self.DIGI_TABLE_SPECS) self.mc_group = self.create_group(self.root, 'mc', 'Monte Carlo') @@ -561,7 +561,7 @@ class DigiOutputFile(OutputFileBase): """ _FILE_TYPE = FileType.DIGI - _READOUT_MODE = HexagonalReadoutMode.RECTANGULAR + #_READOUT_MODE = HexagonalReadoutMode.RECTANGULAR DIGI_TABLE_SPECS = ('digi_table', DigiDescription, 'Digi data') PHA_ARRAY_SPECS = ('pha', tables.Int32Atom(shape=())) MC_TABLE_SPECS = ('mc_table', MonteCarloDescription, 'Monte Carlo data') @@ -570,7 +570,7 @@ def __init__(self, file_path: str): """Constructor. """ super().__init__(file_path) - self.update_header(readoutmode=self._READOUT_MODE.value) + #self.update_header(readoutmode=self._READOUT_MODE.value) self.digi_group = self.create_group(self.root, 'digi', 'Digi') self.digi_table = self.create_table(self.digi_group, *self.DIGI_TABLE_SPECS) self.pha_array = self.create_vlarray(self.digi_group, *self.PHA_ARRAY_SPECS) @@ -604,7 +604,7 @@ def flush(self) -> None: class ReconOutputFile(OutputFileBase): - """Description of a reconstructed output file. This should be the same for + """Description of a reconstructed output file. This should be the same for all types of DigiEvent. Arguments @@ -688,7 +688,7 @@ def header_value(self, key: str, default: Any = None) -> Any: """Return the header value corresponding to a given key. """ return self.header.get(key, default) - + class DigiInputFileBase(InputFileBase): def __init__(self, file_path: str): """Constructor. @@ -697,7 +697,7 @@ def __init__(self, file_path: str): self.digi_table = self.root.digi.digi_table self.mc_table = self.root.mc.mc_table self.__index = -1 - + def column(self, name: str) -> np.ndarray: """Return a given column in the digi table. """ @@ -707,7 +707,7 @@ def mc_column(self, name: str) -> np.ndarray: """Return a given column in the Monte Carlo table. """ return self.mc_table.col(name) - + def mc_event(self, row_index: int) -> MonteCarloEvent: """Random access to the MonteCarloEvent part of the event contribution. @@ -724,7 +724,7 @@ def __iter__(self): """ self.__index = -1 return self - + def __next__(self) -> DigiEventBase: """Overloaded method for the implementation of the iterator protocol. """ @@ -776,7 +776,7 @@ def __next__(self) -> DigiEventSparse: if self.__index == len(self.digi_table): raise StopIteration return self.digi_event(self.__index) - + class DigiInputFileRectangular(DigiInputFileBase): """Description of a rectangular digitized input file. @@ -855,7 +855,7 @@ class DigiInputFile(InputFileBase): """Description of a digitized input file. NOTE: this class should be eliminated when the above three classes have been - fully implemented and tested. + fully implemented and tested. This has a very simple interface: we cache references to the relevant tables when we open the file and we provide methods to reassemble a specific table @@ -960,7 +960,7 @@ def peek_file_type(file_path: str) -> FileType: return FileType(input_file.root.header._v_attrs['filetype']) except KeyError as exception: raise RuntimeError(f'File {file_path} has no type information.') from exception - + def peek_readout_type(file_path: str) -> HexagonalReadoutMode: """Peek into the header of a HDF5 Digi file and determing the readout type. From a35f823b0bfbdce8e67ee05defc23b9f7d58bf36 Mon Sep 17 00:00:00 2001 From: chiara Date: Tue, 30 Apr 2024 13:58:04 +0200 Subject: [PATCH 49/52] minor --- hexsample/clustering.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/hexsample/clustering.py b/hexsample/clustering.py index f84f974..892286e 100644 --- a/hexsample/clustering.py +++ b/hexsample/clustering.py @@ -126,7 +126,7 @@ def run(self, event: DigiEventSparse | DigiEventRectangular | DigiEventCircular) row.append(_row) # ... transforming the coordinates of the NN in its corresponding ADC channel ... adc_channel_order.append(self.grid.adc_channel(_col, _row)) - # ... reordering the pha array for the correspondance (col[i], row[i]) with pha[i]. + # ... reordering the pha array for the correspondance (col[i], row[i]) with pha[i]. pha = event.pha[adc_channel_order] # Converting lists into numpy arrays col = np.array(col) @@ -149,6 +149,8 @@ def run(self, event: DigiEventSparse | DigiEventRectangular | DigiEventCircular) # trick argsort into sorting values in decreasing order. idx = np.argsort(-pha) # Only pick the seed and the N highest pixels. + # This is useless for the circular readout because in that case all + # neighbors are used for track reconstruction. mask = idx[:self.num_neighbors + 1] # If there's any zero left in the target pixels, get rid of it. mask = mask[pha[mask] > 0] From 6e21816222020903df552362e8f9986a84a51b29 Mon Sep 17 00:00:00 2001 From: chiara Date: Tue, 30 Apr 2024 14:06:12 +0200 Subject: [PATCH 50/52] hexsample/bin/hxplots.py --- hexsample/bin/hxplots.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100644 => 100755 hexsample/bin/hxplots.py diff --git a/hexsample/bin/hxplots.py b/hexsample/bin/hxplots.py old mode 100644 new mode 100755 From 34827bc1458524de3a6d9f90c73e771b95bbf2d7 Mon Sep 17 00:00:00 2001 From: Luca Baldini Date: Tue, 30 Apr 2024 14:19:13 +0200 Subject: [PATCH 51/52] Fixed unit tests. --- hexsample/bin/hxsim.py | 3 ++- tests/test_fileio.py | 2 ++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/hexsample/bin/hxsim.py b/hexsample/bin/hxsim.py index f37ac50..8326af7 100755 --- a/hexsample/bin/hxsim.py +++ b/hexsample/bin/hxsim.py @@ -30,7 +30,8 @@ from hexsample import HEXSAMPLE_DATA from hexsample.app import ArgumentParser from hexsample.readout import HexagonalReadoutMode, readout_chip -from hexsample.fileio import DigiDescriptionSparse, DigiDescriptionRectangular, DigiDescriptionCircular, digioutput_class +from hexsample.fileio import DigiDescriptionSparse, DigiDescriptionRectangular,\ + DigiDescriptionCircular, digioutput_class from hexsample.hexagon import HexagonalLayout from hexsample.mc import PhotonList from hexsample.roi import Padding diff --git a/tests/test_fileio.py b/tests/test_fileio.py index 4103b79..69ee195 100644 --- a/tests/test_fileio.py +++ b/tests/test_fileio.py @@ -167,6 +167,7 @@ def test_file_type(): # Test for the digi files. file_path = HEXSAMPLE_DATA / 'test_digi_filetype_rectangular.h5' digi_file = DigiOutputFileRectangular(file_path) + digi_file.update_header(readoutmode='RECTANGULAR') digi_file.close() digi_file = DigiInputFileRectangular(file_path) assert digi_file.file_type == FileType.DIGI @@ -178,6 +179,7 @@ def test_file_type(): # Test for the digi sparse files. file_path = HEXSAMPLE_DATA / 'test_digi_filetype_sparse_readouttype.h5' digi_file = DigiOutputFileSparse(file_path) + digi_file.update_header(readoutmode='SPARSE') digi_file.close() digi_file = DigiInputFileSparse(file_path) assert digi_file.file_type == FileType.DIGI From 18dd12d25bdf5db701b77ef4ad6072f7af15e3e4 Mon Sep 17 00:00:00 2001 From: Luca Baldini Date: Tue, 30 Apr 2024 14:21:26 +0200 Subject: [PATCH 52/52] Updated release notes. --- docs/release.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/release.rst b/docs/release.rst index 808632d..af717e2 100644 --- a/docs/release.rst +++ b/docs/release.rst @@ -4,6 +4,8 @@ Release notes ============= +* Merging https://github.com/lucabaldini/hexsample/pull/46 +* Implemented circular and sparse readout modes. * Minor tweaks to the docs to make clear that all the lengths are in cm. * Minor formatting changes.