Compare commits

1 Commits

Author SHA1 Message Date
Erik Fröjdh
9bf21244b0 added custom io for my302 2025-11-27 17:17:43 +01:00
60 changed files with 1157 additions and 2441 deletions

View File

@@ -387,7 +387,6 @@ set(SourceFiles
${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/to_string.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/utils/task.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/utils/ifstream_helpers.cpp
)
@@ -461,7 +460,6 @@ if(AARE_TESTS)
${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/utils/task.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/to_string.test.cpp
)
target_sources(tests PRIVATE ${TestSources} )

View File

@@ -1,17 +1,5 @@
# Release notes
## head
### New Features:
- Expanding 24 to 32 bit data
- Decoding digital data from Mythen 302
- added ``transform_eta_values``. Function transforms :math:`\eta` to uniform spatial coordinates. Should only be used for easier debugging.
- New to_string, string_to for aare
- Added exptime and period members to RawMasterFile including decoding
- Removed redundant arr.value(ix,iy...) on NDArray use arr(ix,iy...)
- Removed Print/Print_some/Print_all form NDArray (operator << still works)
- Added const* version of .data()
### 2025.11.21

View File

@@ -12,19 +12,15 @@ set(SPHINX_SOURCE ${CMAKE_CURRENT_SOURCE_DIR}/src)
set(SPHINX_BUILD ${CMAKE_CURRENT_BINARY_DIR})
file(GLOB_RECURSE SPHINX_SOURCE_FILES
CONFIGURE_DEPENDS
RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}/src"
"${CMAKE_CURRENT_SOURCE_DIR}/src/*.rst"
)
file(GLOB SPHINX_SOURCE_FILES CONFIGURE_DEPENDS "src/*.rst")
foreach(relpath IN LISTS SPHINX_SOURCE_FILES)
set(src "${CMAKE_CURRENT_SOURCE_DIR}/src/${relpath}")
set(dst "${SPHINX_BUILD}/src/${relpath}")
message(STATUS "Copying ${src} to ${dst}")
configure_file("${src}" "${dst}" COPYONLY)
endforeach()
foreach(filename ${SPHINX_SOURCE_FILES})
get_filename_component(fname ${filename} NAME)
message(STATUS "Copying ${filename} to ${SPHINX_BUILD}/src/${fname}")
configure_file(${filename} "${SPHINX_BUILD}/src/${fname}")
endforeach(filename ${SPHINX_SOURCE_FILES})
configure_file(
"${CMAKE_CURRENT_SOURCE_DIR}/conf.py.in"

Binary file not shown.

Before

Width:  |  Height:  |  Size: 9.3 KiB

After

Width:  |  Height:  |  Size: 6.7 KiB

Binary file not shown.

Binary file not shown.

Before

Width:  |  Height:  |  Size: 14 KiB

After

Width:  |  Height:  |  Size: 10 KiB

Binary file not shown.

Binary file not shown.

Before

Width:  |  Height:  |  Size: 20 KiB

After

Width:  |  Height:  |  Size: 13 KiB

View File

@@ -1,43 +1,12 @@
.. _Interpolation_C++API:
Interpolation
==============
The Interpolation class implements the :math:`\eta`-interpolation method.
This interpolation technique is based on charge sharing: for detected photon hits (e.g. clusters), it refines the estimated photon hit using information from neighboring pixels.
Interpolation class for :math:`\eta` Interpolation.
The method relies on the so-called :math:`\eta`-functions, which describe the relationship between the energy measured in the central cluster pixel (the initially estimated photon hit) and the energies measured in its neighboring pixels.
Depending on how much energy each neighboring pixel receives relative to the central pixel, the estimated photon hit is shifted toward that neighbor by a certain offset to the actual photon hit position in the pixel :math:`(x, y)`.
The mapping between the :math:`\eta` values and the corresponding spatial photon position :math:`(x,y)` can be viewed as an optimal transport problem.
One can readily compute the probability distribution :math:`P_{\eta}` of the :math:`\eta` values by forming a 2D histogram.
However, the probability distribution :math:`P_{x,y}` of the true photon positions is generally unknown unless the detector is illuminated uniformly (i.e. under flat-field conditions).
In a flat-field, the photon positions are uniformly distributed.
With this assumption, the problem reduces to determining a transport map :math:`T:(\eta_x,\eta_y) \rightarrow (x,y)`, that pushes forward the distribution of :math:`(\eta_x, \eta_y)` to the known uniform distribution of photon positions of a flatfield.
The map :math:`T` is given by:
.. math::
\begin{align*}
T_1: & F_{x}^{-1} F_{\eta_x|\eta_y} \\
T_2: & F_{y}^{-1} F_{\eta_y|\eta_x},
\end{align*}
where :math:`F_{\eta_x|\eta_y}` and :math:`F_{\eta_y|\eta_x}` are the conditional cumulative distribution functions e.g. :math:`F_{\eta_x|\eta_y}(\eta_x', \eta_y') = P_{\eta_x, \eta_y}(\eta_x \leq \eta_x' | \eta_y = \eta_y')`.
And :math:`F_{x}` and :math:`F_{y}` are the cumulative distribution functions of :math:`x` and :math:`y`. Note as :math:`x` and :math:`y` are uniformly distributed :math:`F_{x}` and :math:`F_{y}` are the identity functions. The map :math:`T` thus simplifies to
.. math::
\begin{align*}
T_1: & F_{\eta_x|\eta_y} \\
T_2: & F_{\eta_y|\eta_x}.
\end{align*}
Note that for the implementation :math:`P_{\eta}` is not only a distribution of :math:`\eta_x`, :math:`\eta_y` but also of the estimated photon energy :math:`e`.
The energy level correlates slightly with the z-depth. Higher z-depth leads to more charge sharing and a different :math:`\eta` distribution. Thus we create a mapping :math:`T` for each energy level.
The Interpolator class provides methods to interpolate the positions of photons based on their :math:`\eta` values.
.. warning::
The interpolation might lead to erroneous photon positions for clusters at the boarders of a frame. Make sure to filter out such cases.
:math:`\eta`-Functions:
---------------------------
@@ -52,8 +21,6 @@ The energy level correlates slightly with the z-depth. Higher z-depth leads to m
Supported are the following :math:`\eta`-functions:
:math:`\eta`-Function on 2x2 Clusters:
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. image:: ../figures/Eta2x2.png
:target: ../figures/Eta2x2.png
@@ -67,18 +34,11 @@ Supported are the following :math:`\eta`-functions:
{\color{green}{\eta_y}} = \frac{Q_{1,1}}{Q_{0,1} + Q_{1,1}}
\end{equation*}
The :math:`\eta` values can range between 0,1. Note they only range between 0,1 because the position of the center pixel (red) can change.
If the center pixel is in the bottom left pixel :math:`\eta_x` will be close to zero. If the center pixel is in the bottom right pixel :math:`\eta_y` will be close to 1.
One can apply this :math:`\eta` not only on 2x2 clusters but on clusters with any size. Then the 2x2 subcluster with maximum energy is choosen and the :math:`\eta` function applied on the subcluster.
.. doxygenfunction:: aare::calculate_eta2(const ClusterVector<ClusterType>&)
.. doxygenfunction:: aare::calculate_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>&)
Full :math:`\eta`-Function on 2x2 Clusters:
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. image:: ../figures/Eta2x2Full.png
:target: ../figures/Eta2x2Full.png
:width: 650px
@@ -87,20 +47,15 @@ Full :math:`\eta`-Function on 2x2 Clusters:
.. math::
\begin{equation*}
{\color{blue}{\eta_x}} = \frac{Q_{0,1} + Q_{1,1}}{\sum_i^{1}\sum_j^{1}Q_{i,j}} \quad \quad
{\textcolor{green}{\eta_y}} = \frac{Q_{1,0} + Q_{1,1}}{\sum_i^{1}\sum_j^{1}Q_{i,j}}
{\color{blue}{\eta_x}} = \frac{Q_{0,1} + Q_{1,1}}{\sum_i^{2}\sum_j^{2}Q_{i,j}} \quad \quad
{\textcolor{green}{\eta_y}} = \frac{Q_{1,0} + Q_{1,1}}{\sum_i^{2}\sum_j^{2}Q_{i,j}}
\end{equation*}
The :math:`\eta` values can range between 0,1. Note they only range between 0,1 because the position of the center pixel (red) can change.
If the center pixel is in the bottom left pixel :math:`\eta_x` will be close to zero. If the center pixel is in the bottom right pixel :math:`\eta_y` will be close to 1.
.. doxygenfunction:: aare::calculate_full_eta2(const ClusterVector<ClusterType>&)
.. doxygenfunction:: aare::calculate_full_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>&)
Full :math:`\eta`-Function on 3x3 Clusters:
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. image:: ../figures/Eta3x3.png
:target: ../figures/Eta3x3.png
:width: 650px
@@ -109,18 +64,13 @@ Full :math:`\eta`-Function on 3x3 Clusters:
.. math::
\begin{equation*}
{\color{blue}{\eta_x}} = \frac{\sum_{i=0}^{2} Q_{i,2} - \sum_{i=0}^{2} Q_{i,0}}{\sum_{i=0}^{2}\sum_{j=0}^{2} Q_{i,j}} \quad \quad
{\color{green}{\eta_y}} = \frac{\sum_{j=0}^{2} Q_{2,j} - \sum_{j=0}^{2} Q_{0,j}}{\sum_{i=0}^{2}\sum_{j=0}^{2} Q_{i,j}}
{\color{blue}{\eta_x}} = \frac{\sum_{i}^{3} Q_{i,2} - \sum_{i}^{3} Q_{i,0}}{\sum_{i}^{3}\sum_{j}^{3} Q_{i,j}} \quad \quad
{\color{green}{\eta_y}} = \frac{\sum_{j}^{3} Q_{2,j} - \sum_{j}^{3} Q_{0,j}}{\sum_{i}^{3}\sum_{j}^{3} Q_{i,j}}
\end{equation*}
The :math:`\eta` values can range between -0.5,0.5.
.. doxygenfunction:: aare::calculate_eta3(const ClusterVector<Cluster<T, 3,3, CoordType>>&)
.. doxygenfunction:: aare::calculate_eta3(const ClusterVector<ClusterType>&)
.. doxygenfunction:: aare::calculate_eta3(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>&)
Cross :math:`\eta`-Function on 3x3 Clusters:
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. doxygenfunction:: aare::calculate_eta3(const Cluster<T, 3, 3, CoordType>&)
.. image:: ../figures/Eta3x3Cross.png
:target: ../figures/Eta3x3Cross.png
@@ -130,28 +80,20 @@ Cross :math:`\eta`-Function on 3x3 Clusters:
.. math::
\begin{equation*}
{\color{blue}{\eta_x}} = \frac{Q_{1,2} - Q_{1,0}}{Q_{1,0} + Q_{1,1} + Q_{1,2}} \quad \quad
{\color{green}{\eta_y}} = \frac{Q_{0,2} - Q_{0,1}}{Q_{0,1} + Q_{1,1} + Q_{2,1}}
{\color{blue}{\eta_x}} = \frac{Q_{1,2} - Q_{1,0}}{Q_{1,0} + Q_{1,1} + Q_{1,0}} \quad \quad
{\color{green}{\eta_y}} = \frac{Q_{0,2} - Q_{0,1}}{Q_{0,1} + Q_{1,1} + Q_{1,2}}
\end{equation*}
The :math:`\eta` values can range between -0.5,0.5.
.. doxygenfunction:: aare::calculate_cross_eta3(const ClusterVector<Cluster<T, 3,3, CoordType>>&)
.. doxygenfunction:: aare::calculate_cross_eta3(const ClusterVector<ClusterType>&)
.. doxygenfunction:: aare::calculate_cross_eta3(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>&)
.. doxygenfunction:: aare::calculate_cross_eta3(const Cluster<T, 3, 3, CoordType>&)
Interpolation class:
---------------------
.. warning::
The interpolation might lead to erroneous photon positions for clusters at the borders of a frame. Make sure to filter out such cases.
.. Warning::
Make sure to use the same :math:`\eta`-function during interpolation as given by the joint :math:`\eta`-distribution passed to the constructor.
.. Note::
Make sure to use resonable energy bins, when constructing the joint distribution. If data is too sparse for a given energy the interpolation will lead to erreneous results.
.. doxygenclass:: aare::Interpolator
:members:
:undoc-members:

View File

@@ -22,14 +22,21 @@ AARE
.. toctree::
:caption: Python API
:maxdepth: 3
:hidden:
:maxdepth: 1
pyFile
pycalibration
python/cluster/index
python/file/index
pyFit
pyCtbRawFile
pyClusterFile
pyClusterVector
pyCluster
pyInterpolation
pyJungfrauDataFile
pyRawFile
pyRawMasterFile
pyVarClusterFinder
pyFit
.. toctree::

11
docs/src/pyCtbRawFile.rst Normal file
View File

@@ -0,0 +1,11 @@
CtbRawFile
============
.. py:currentmodule:: aare
.. autoclass:: CtbRawFile
:members:
:undoc-members:
:show-inheritance:
:inherited-members:

View File

@@ -0,0 +1,94 @@
Interpolation
==============
Interpolation class for :math:`\eta` Interpolation.
The Interpolator class provides methods to interpolate the positions of photons based on their :math:`\eta` values.
.. warning::
The interpolation might lead to erroneous photon positions for clusters at the boarders of a frame. Make sure to filter out such cases.
Below is an example of the Eta class of type ``double``. Supported are ``Etaf`` of type ``float`` and ``Etai`` of type ``int``.
.. autoclass:: aare._aare.Etad
:members:
:private-members:
.. note::
The corner value ``c`` is only relevant when one uses ``calculate_eta_2`` or ``calculate_full_eta2``. Otherwise its default value is ``cTopLeft``.
Supported are the following :math:`\eta`-functions:
.. py:currentmodule:: aare
.. image:: ../figures/Eta2x2.png
:target: ../figures/Eta2x2.png
:width: 650px
:align: center
:alt: Eta2x2
.. math::
\begin{equation*}
{\color{blue}{\eta_x}} = \frac{Q_{1,1}}{Q_{1,0} + Q_{1,1}} \quad \quad
{\color{green}{\eta_y}} = \frac{Q_{1,1}}{Q_{0,1} + Q_{1,1}}
\end{equation*}
.. autofunction:: calculate_eta2
.. image:: ../figures/Eta2x2Full.png
:target: ../figures/Eta2x2Full.png
:width: 650px
:align: center
:alt: Eta2x2 Full
.. math::
\begin{equation*}
{\color{blue}{\eta_x}} = \frac{Q_{0,1} + Q_{1,1}}{\sum_i^{2}\sum_j^{2}Q_{i,j}} \quad \quad
{\textcolor{green}{\eta_y}} = \frac{Q_{1,0} + Q_{1,1}}{\sum_i^{2}\sum_j^{2}Q_{i,j}}
\end{equation*}
.. autofunction:: calculate_full_eta2
.. image:: ../figures/Eta3x3.png
:target: ../figures/Eta3x3.png
:width: 650px
:align: center
:alt: Eta3x3
.. math::
\begin{equation*}
{\color{blue}{\eta_x}} = \frac{\sum_{i}^{3} Q_{i,2} - \sum_{i}^{3} Q_{i,0}}{\sum_{i}^{3}\sum_{j}^{3} Q_{i,j}} \quad \quad
{\color{green}{\eta_y}} = \frac{\sum_{j}^{3} Q_{2,j} - \sum_{j}^{3} Q_{0,j}}{\sum_{i}^{3}\sum_{j}^{3} Q_{i,j}}
\end{equation*}
.. autofunction:: calculate_eta3
.. image:: ../figures/Eta3x3Cross.png
:target: ../figures/Eta3x3Cross.png
:width: 650px
:align: center
:alt: Cross Eta3x3
.. math::
\begin{equation*}
{\color{blue}{\eta_x}} = \frac{Q_{1,2} - Q_{1,0}}{Q_{1,0} + Q_{1,1} + Q_{1,0}} \quad \quad
{\color{green}{\eta_y}} = \frac{Q_{0,2} - Q_{0,1}}{Q_{0,1} + Q_{1,1} + Q_{1,2}}
\end{equation*}
.. autofunction:: calculate_cross_eta3
Interpolation class for :math:`\eta`-Interpolation
----------------------------------------------------
.. Warning::
Make sure to use the same :math:`\eta`-function during interpolation as given by the joint :math:`\eta`-distribution passed to the constructor.
.. py:currentmodule:: aare
.. autoclass:: Interpolator
:special-members: __init__
:members:
:undoc-members:
:inherited-members:

View File

@@ -1,11 +0,0 @@
Cluster & Interpolation
==========================
.. toctree::
:caption: Cluster & Interpolation
:maxdepth: 1
pyCluster
pyClusterVector
pyInterpolation
pyVarClusterFinder

View File

@@ -1,124 +0,0 @@
Interpolation
==============
The Interpolation class implements the :math:`\eta`-interpolation method.
This interpolation technique is based on charge sharing: for detected photon hits (e.g. clusters), it refines the estimated photon hit using information from neighboring pixels.
See :ref:`Interpolation_C++API` for a more elaborate documentation and explanation of the method.
:math:`\eta`-Functions:
--------------------------
Below is an example of the Eta class of type ``double``. Supported are ``Etaf`` of type ``float`` and ``Etai`` of type ``int``.
.. autoclass:: aare._aare.Etad
:members:
:private-members:
.. note::
The corner value ``c`` is only relevant when one uses ``calculate_eta_2`` or ``calculate_full_eta2``. Otherwise its default value is ``cTopLeft``.
Supported are the following :math:`\eta`-functions:
:math:`\eta`-Function on 2x2 Clusters:
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. py:currentmodule:: aare
.. image:: ../../../figures/Eta2x2.png
:target: ../../../figures/Eta2x2.png
:width: 650px
:align: center
:alt: Eta2x2
.. math::
\begin{equation*}
{\color{blue}{\eta_x}} = \frac{Q_{1,1}}{Q_{1,0} + Q_{1,1}} \quad \quad
{\color{green}{\eta_y}} = \frac{Q_{1,1}}{Q_{0,1} + Q_{1,1}}
\end{equation*}
The :math:`\eta` values can range between 0,1. Note they only range between 0,1 because the position of the center pixel (red) can change.
If the center pixel is in the bottom left pixel :math:`\eta_x` will be close to zero. If the center pixel is in the bottom right pixel :math:`\eta_y` will be close to 1.
.. autofunction:: calculate_eta2
Full :math:`\eta`-Function on 2x2 Clusters:
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. image:: ../../../figures/Eta2x2Full.png
:target: ../../../figures/Eta2x2Full.png
:width: 650px
:align: center
:alt: Eta2x2 Full
.. math::
\begin{equation*}
{\color{blue}{\eta_x}} = \frac{Q_{0,1} + Q_{1,1}}{\sum_{i=0}^{1}\sum_{j=0}^{1}Q_{i,j}} \quad \quad
{\textcolor{green}{\eta_y}} = \frac{Q_{1,0} + Q_{1,1}}{\sum_{i=0}^{1}\sum_{j=0}^{1}Q_{i,j}}
\end{equation*}
The :math:`\eta` values can range between 0,1. Note they only range between 0,1 because the position of the center pixel (red) can change.
If the center pixel is in the bottom left pixel :math:`\eta_x` will be close to zero. If the center pixel is in the bottom right pixel :math:`\eta_y` will be close to 1.
.. autofunction:: calculate_full_eta2
Full :math:`\eta`-Function on 3x3 Clusters:
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. image:: ../../../figures/Eta3x3.png
:target: ../../../figures/Eta3x3.png
:width: 650px
:align: center
:alt: Eta3x3
.. math::
\begin{equation*}
{\color{blue}{\eta_x}} = \frac{\sum_{i=0}^{2} Q_{i,2} - \sum_{i=0}^{2} Q_{i,0}}{\sum_{i=0}^{2}\sum_{j}^{3} Q_{i,j}} \quad \quad
{\color{green}{\eta_y}} = \frac{\sum_{j=0}^{2} Q_{2,j} - \sum_{j=0}^{2} Q_{0,j}}{\sum_{i=0}^{2}\sum_{j}^{3} Q_{i,j}}
\end{equation*}
The :math:`\eta` values can range between -0.5,0.5.
.. autofunction:: calculate_eta3
Cross :math:`\eta`-Function on 3x3 Clusters:
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. image:: ../../../figures/Eta3x3Cross.png
:target: ../../../figures/Eta3x3Cross.png
:width: 650px
:align: center
:alt: Cross Eta3x3
.. math::
\begin{equation*}
{\color{blue}{\eta_x}} = \frac{Q_{1,2} - Q_{1,0}}{Q_{1,0} + Q_{1,1} + Q_{1,2}} \quad \quad
{\color{green}{\eta_y}} = \frac{Q_{0,2} - Q_{0,1}}{Q_{0,1} + Q_{1,1} + Q_{2,1}}
\end{equation*}
The :math:`\eta` values can range between -0.5,0.5.
.. autofunction:: calculate_cross_eta3
Interpolation class for :math:`\eta`-Interpolation
----------------------------------------------------
.. Warning::
Make sure to use the same :math:`\eta`-function during interpolation as given by the joint :math:`\eta`-distribution passed to the constructor.
.. Warning::
The interpolation might lead to erroneous photon positions for clusters at the boarders of a frame. Make sure to filter out such cases.
.. Note::
Make sure to use resonable energy bins, when constructing the joint distribution. If data is too sparse for a given energy the interpolation will lead to erreneous results.
.. py:currentmodule:: aare
.. autoclass:: Interpolator
:special-members: __init__
:members:
:undoc-members:
:inherited-members:

View File

@@ -1,14 +0,0 @@
File I/O
===================
.. toctree::
:caption: File I/O
:maxdepth: 1
pyClusterFile
pyCtbRawFile
pyFile
pyJungfrauDataFile
pyRawFile
pyRawMasterFile
pyTransform

View File

@@ -1,25 +0,0 @@
CtbRawFile
============
Read analog, digital and transceiver samples from a raw file containing
data from the Chip Test Board. Uses :mod:`aare.transform` to decode the
data into a format that the user can work with.
.. code:: python
import aare
from aare.transform import Mythen302Transform
my302 = Mythen302Transform(offset = 4)
with aare.CtbRawFile(fname, transform = my302) as f:
for header, data in f:
#do something with the data
.. py:currentmodule:: aare
.. autoclass:: CtbRawFile
:members:
:undoc-members:
:show-inheritance:
:inherited-members:

View File

@@ -1,27 +0,0 @@
Transform
===================
The transform module takes data read by :class:`aare.CtbRawFile` and decodes it
to a useful image format. Depending on detector it supports both analog
and digital samples.
For convenience the following transform objects are defined with a short name
.. code:: python
moench05 = Moench05Transform()
moench05_1g = Moench05Transform1g()
moench05_old = Moench05TransformOld()
matterhorn02 = Matterhorn02Transform()
adc_sar_04_64to16 = AdcSar04Transform64to16()
adc_sar_05_64to16 = AdcSar05Transform64to16()
.. py:currentmodule:: aare
.. automodule:: aare.transform
:members:
:undoc-members:
:private-members:
:special-members: __call__
:show-inheritance:
:inherited-members:

View File

@@ -29,8 +29,7 @@ template <typename T> struct Eta2 {
double x{};
/// @brief eta in y direction
double y{};
/// @brief index of subcluster with highest energy value (given as corner
/// relative to cluster center)
/// @brief index of subcluster given as corner relative to cluster center
corner c{0};
/// @brief photon energy (cluster sum)
T sum{};

View File

@@ -34,15 +34,15 @@ struct FileConfig {
DetectorType detector_type{DetectorType::Unknown};
int max_frames_per_file{};
size_t total_frames{};
// std::string to_string() const {
// return "{ dtype: " + dtype.to_string() +
// ", rows: " + std::to_string(rows) +
// ", cols: " + std::to_string(cols) +
// ", geometry: " + geometry.to_string() +
// ", detector_type: " + ToString(detector_type) +
// ", max_frames_per_file: " + std::to_string(max_frames_per_file) +
// ", total_frames: " + std::to_string(total_frames) + " }";
// }
std::string to_string() const {
return "{ dtype: " + dtype.to_string() +
", rows: " + std::to_string(rows) +
", cols: " + std::to_string(cols) +
", geometry: " + geometry.to_string() +
", detector_type: " + ToString(detector_type) +
", max_frames_per_file: " + std::to_string(max_frames_per_file) +
", total_frames: " + std::to_string(total_frames) + " }";
}
};
/**

View File

@@ -17,27 +17,11 @@ struct Photon {
double energy;
};
struct Coordinate2D {
double x{};
double y{};
};
class Interpolator {
/**
* @brief
* marginal CDF of eta_x (if rosenblatt applied), conditional
* CDF of eta_x conditioned on eta_y
* value at (i, j, e): F(eta_x[i] |
*eta_y[j], energy[e])
*/
// marginal CDF of eta_x (if rosenblatt applied), conditional
// CDF of eta_x conditioned on eta_y
NDArray<double, 3> m_ietax;
/**
* @brief
* conditional CDF of eta_y conditioned on eta_x
* value at (i,j,e): F(eta_y[j] | eta_x[i], energy[e])
*/
// conditional CDF of eta_y conditioned on eta_x
NDArray<double, 3> m_ietay;
NDArray<double, 1> m_etabinsx;
@@ -47,11 +31,11 @@ class Interpolator {
public:
/**
* @brief Constructor for the Interpolator class
* @param etacube joint distribution of etaX, etaY and photon energy (note
* first dimension is etaX, second etaY, third photon energy)
* @param etacube joint distribution of etaX, etaY and photon energy
* @param xbins bin edges for etaX
* @param ybins bin edges for etaY
* @param ebins bin edges for photon energy
* @note note first dimension is etaX, second etaY, third photon energy
*/
Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
NDView<double, 1> ybins, NDView<double, 1> ebins);
@@ -69,8 +53,8 @@ class Interpolator {
* @brief transforms the joint eta distribution of etaX and etaY to the two
* independant uniform distributions based on the Roseblatt transform for
* each energy level
* @param etacube joint distribution of etaX, etaY and photon energy (first
* dimension is etaX, second etaY, third photon energy)
* @param etacube joint distribution of etaX, etaY and photon energy
* @note note first dimension is etaX, second etaY, third photon energy
*/
void rosenblatttransform(NDView<double, 3> etacube);
@@ -85,24 +69,25 @@ class Interpolator {
* calculate_eta2
* @return interpolated photons (photon positions are given as double but
* following row column format e.g. x=0, y=0 means top row and first column
* of frame) (An interpolated photon position of (1.5, 2.5) corresponds to
* an estimated photon hit at the pixel center of pixel (1,2))
* of frame)
*/
template <auto EtaFunction = calculate_eta2, typename ClusterType,
typename Eanble = std::enable_if_t<is_cluster_v<ClusterType>>>
std::vector<Photon>
interpolate(const ClusterVector<ClusterType> &clusters) const;
/**
* @brief transforms the eta values to uniform coordinates based on the CDF
* ieta_x and ieta_y
* @tparam eta Eta to transform
* @return uniform coordinates {x,y}
*/
template <typename T>
Coordinate2D transform_eta_values(const Eta2<T> &eta) const;
std::vector<Photon> interpolate(const ClusterVector<ClusterType> &clusters);
private:
/**
* @brief implements underlying interpolation logic based on EtaFunction
* Type
* @tparam EtaFunction Function object that calculates desired eta default:
* @param u: transformed photon position in x between [0,1]
* @param v: transformed photon position in y between [0,1]
* @param c: corner of eta
*/
template <auto EtaFunction, typename ClusterType>
void interpolation_logic(Photon &photon, const double u, const double v,
const corner c = corner::cTopLeft);
/**
* @brief bilinear interpolation of the transformed eta values
* @param ix index of etaX bin
@@ -113,14 +98,13 @@ class Interpolator {
template <typename T>
std::pair<double, double>
bilinear_interpolation(const size_t ix, const size_t iy, const size_t ie,
const Eta2<T> &eta) const;
const Eta2<T> &eta);
};
template <typename T>
std::pair<double, double>
Interpolator::bilinear_interpolation(const size_t ix, const size_t iy,
const size_t ie,
const Eta2<T> &eta) const {
const size_t ie, const Eta2<T> &eta) {
auto next_index_y = static_cast<ssize_t>(iy + 1) >= m_ietax.shape(1)
? m_ietax.shape(1) - 1
: iy + 1;
@@ -160,36 +144,9 @@ Interpolator::bilinear_interpolation(const size_t ix, const size_t iy,
return {ietax_interpolated, ietay_interpolated};
}
template <typename T>
Coordinate2D Interpolator::transform_eta_values(const Eta2<T> &eta) const {
// Finding the index of the last element that is smaller
// should work fine as long as we have many bins
auto ie = last_smaller(m_energy_bins, static_cast<double>(eta.sum));
auto ix = last_smaller(m_etabinsx, eta.x);
auto iy = last_smaller(m_etabinsy, eta.y);
if (static_cast<ssize_t>(ix) >= m_etabinsx.size() - 1 ||
static_cast<ssize_t>(iy) >= m_etabinsy.size() - 1 ||
static_cast<ssize_t>(ie) >= m_energy_bins.size() - 1)
throw std::runtime_error(
fmt::format("Eta values out of bounds of eta distribution: eta.x = "
"{:.4f}, eta.y = {:.4f}, energy = {:.4f}",
eta.x, eta.y, eta.sum));
// TODO: bilinear interpolation only works if all bins have a size > 1 -
// otherwise bilinear interpolation with zero values which skew the
// results
// TODO: maybe trim the bins at the edges with zero values beforehand
// auto [ietax_interpolated, ietay_interpolated] =
// bilinear_interpolation(ix, iy, ie, eta);
return Coordinate2D{m_ietax(ix, iy, ie), m_ietay(ix, iy, ie)};
}
template <auto EtaFunction, typename ClusterType, typename Enable>
std::vector<Photon>
Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) const {
Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) {
std::vector<Photon> photons;
photons.reserve(clusters.size());
@@ -202,21 +159,55 @@ Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) const {
photon.y = cluster.y;
photon.energy = static_cast<decltype(photon.energy)>(eta.sum);
auto uniform_coordinates = transform_eta_values(eta);
// std::cout << "eta.x: " << eta.x << " eta.y: " << eta.y << std::endl;
// Finding the index of the last element that is smaller
// should work fine as long as we have many bins
auto ie = last_smaller(m_energy_bins, photon.energy);
auto ix = last_smaller(m_etabinsx, eta.x);
auto iy = last_smaller(m_etabinsy, eta.y);
// std::cout << "ix: " << ix << " iy: " << iy << std::endl;
// TODO: bilinear interpolation only works if all bins have a size > 1 -
// otherwise bilinear interpolation with zero values which skew the
// results
// TODO: maybe trim the bins at the edges with zero values beforehand
// auto [ietax_interpolated, ietay_interpolated] =
// bilinear_interpolation(ix, iy, ie, eta);
double ietax_interpolated = m_ietax(ix, iy, ie);
double ietay_interpolated = m_ietay(ix, iy, ie);
interpolation_logic<EtaFunction, ClusterType>(
photon, ietax_interpolated, ietay_interpolated, eta.c);
photons.push_back(photon);
}
return photons;
}
template <auto EtaFunction, typename ClusterType>
void Interpolator::interpolation_logic(Photon &photon, const double u,
const double v, const corner c) {
// std::cout << "u: " << u << " v: " << v << std::endl;
// TODO: try to call this with std::is_same_v and have it constexpr if
// possible
if (EtaFunction == &calculate_eta2<typename ClusterType::value_type,
ClusterType::cluster_size_x,
ClusterType::cluster_size_y,
typename ClusterType::coord_type> ||
EtaFunction ==
&calculate_full_eta2<typename ClusterType::value_type,
EtaFunction == &calculate_full_eta2<typename ClusterType::value_type,
ClusterType::cluster_size_x,
ClusterType::cluster_size_y,
typename ClusterType::coord_type>) {
double dX{}, dY{};
// TODO: could also chaneg the sign of the eta calculation
switch (eta.c) {
switch (c) {
case corner::cTopLeft:
dX = -1.0;
dY = -1.0;
@@ -234,21 +225,14 @@ Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) const {
dY = 0.0;
break;
}
photon.x = photon.x + 0.5 + uniform_coordinates.x +
dX; // use pixel center + 0.5
photon.y =
photon.y + 0.5 + uniform_coordinates.y +
photon.x = photon.x + 0.5 + u + dX; // use pixel center + 0.5
photon.y = photon.y + 0.5 + v +
dY; // eta2 calculates the ratio between bottom and sum of
// bottom and top shift by 1 add eta value correctly
} else {
photon.x += uniform_coordinates.x;
photon.y += uniform_coordinates.y;
photon.x += u;
photon.y += v;
}
photons.push_back(photon);
}
return photons;
}
} // namespace aare

View File

@@ -1,10 +1,13 @@
// SPDX-License-Identifier: MPL-2.0
//
// Container holding image data, or a time series of image data in contigious
// memory. Used for all data processing in Aare.
//
#pragma once
/*
Container holding image data, or a time series of image data in contigious
memory.
TODO! Add expression templates for operators
*/
#include "aare/ArrayExpr.hpp"
#include "aare/NDView.hpp"
@@ -23,17 +26,12 @@ template <typename T, ssize_t Ndim = 2>
class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
std::array<ssize_t, Ndim> shape_;
std::array<ssize_t, Ndim> strides_;
size_t size_{}; // TODO! do we need to store size when we have shape?
size_t size_{}; //TODO! do we need to store size when we have shape?
T *data_;
public:
///////////////////////////////////////////////////////////////////////////////
// Constructors
//
///////////////////////////////////////////////////////////////////////////////
/**
* @brief Default constructor. Constructs an empty NDArray.
* @brief Default constructor. Will construct an empty NDArray.
*
*/
NDArray() : shape_(), strides_(c_strides<Ndim>(shape_)), data_(nullptr) {};
@@ -46,7 +44,8 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
*/
explicit NDArray(std::array<ssize_t, Ndim> shape)
: shape_(shape), strides_(c_strides<Ndim>(shape_)),
size_(num_elements(shape_)), data_(new T[size_]) {}
size_(num_elements(shape_)),
data_(new T[size_]) {}
/**
* @brief Construct a new NDArray object with a shape and value.
@@ -58,10 +57,6 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
this->operator=(value);
}
// Allow NDArray of different type and dimension to be friend classes
// This is needed for the move constructor from NDArray<T,Ndim+1>
template <typename U, ssize_t Dim> friend class NDArray;
/**
* @brief Construct a new NDArray object from a NDView.
* @note The data is copied from the view to the NDArray.
@@ -72,67 +67,44 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
std::copy(v.begin(), v.end(), begin());
}
/**
* @brief Construct a new NDArray object from an std::array.
*/
template <size_t Size>
NDArray(const std::array<T, Size> &arr) : NDArray<T, 1>({Size}) {
std::copy(arr.begin(), arr.end(), begin());
}
/**
* @brief Move construct a new NDArray object. Cheap since it just
* reassigns the pointer and copy size/strides.
*
* @param other
*/
// Move constructor
NDArray(NDArray &&other) noexcept
: shape_(other.shape_), strides_(c_strides<Ndim>(shape_)),
size_(other.size_), data_(other.data_) {
other.reset(); // Needed to avoid double free
other.reset(); // TODO! is this necessary?
}
/**
* @brief Move construct a new NDArray object from an array with Ndim + 1.
* Can be used to drop a dimension cheaply.
* @param other
*/
//Move constructor from an an array with Ndim + 1
template <ssize_t M, typename = std::enable_if_t<(M == Ndim + 1)>>
NDArray(NDArray<T, M> &&other)
: shape_(drop_first_dim(other.shape())),
strides_(c_strides<Ndim>(shape_)), size_(num_elements(shape_)),
data_(other.data()) {
// For now only allow move if the size matches, to avoid unreachable
// data if the use case arises we can remove this check
if (size() != other.size()) {
data_ = nullptr; // avoid double free, other will clean up the
// memory in it's destructor
throw std::runtime_error(
LOCATION +
// For now only allow move if the size matches, to avoid unreachable data
// if the use case arises we can remove this check
if(size() != other.size()) {
data_ = nullptr; // avoid double free, other will clean up the memory in it's destructor
throw std::runtime_error(LOCATION +
"Size mismatch in move constructor of NDArray<T, Ndim-1>");
}
other.reset();
}
/**
* @brief Copy construct a new NDArray object from another NDArray.
*
* @param other
*/
// Copy constructor
NDArray(const NDArray &other)
: shape_(other.shape_), strides_(c_strides<Ndim>(shape_)),
size_(other.size_), data_(new T[size_]) {
std::copy(other.data_, other.data_ + size_, data_);
}
/**
* @brief Conversion from a ArrayExpr to an actual NDArray. Used when
* the expression is evaluated and data needed.
*
* @tparam E
* @param expr
*/
// Conversion operator from array expression to array
template <typename E>
NDArray(ArrayExpr<E, Ndim> &&expr) : NDArray(expr.shape()) {
for (size_t i = 0; i < size_; ++i) {
@@ -140,129 +112,23 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
}
}
/**
* @brief Destroy the NDArray object. Frees the allocated memory.
*
*/
~NDArray() { delete[] data_; }
///////////////////////////////////////////////////////////////////////////////
// Iterators and indexing
//
///////////////////////////////////////////////////////////////////////////////
auto begin() { return data_; }
auto end() { return data_ + size_; }
auto *begin() { return data_; }
const auto *begin() const { return data_; }
auto *end() { return data_ + size_; }
const auto *end() const { return data_ + size_; }
/*
* @brief Access element at given multi-dimensional index.
* i.e. arr(i,j,k,...)
*
* @note The fast index is the last index. Please take care when iterating
* through the array.
*/
template <typename... Ix>
std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) {
return data_[element_offset(strides_, index...)];
}
/*
* @brief Access element at given multi-dimensional index (const version).
* i.e. arr(i,j,k,...)
*
* @note The fast index is the last index. Please take care when iterating
* through the array.
*/
template <typename... Ix>
std::enable_if_t<sizeof...(Ix) == Ndim, const T &>
operator()(Ix... index) const {
return data_[element_offset(strides_, index...)];
}
/*
@brief Index the array as it would be a 1D array. To get a certain
pixel in a multidimensional array use the (i,j,k,...) operator instead.
*/
T &operator()(ssize_t i) { return data_[i]; }
/*
@brief Index the array as it would be a 1D array. To get a certain
pixel in a multidimensional array use the (i,j,k,...) operator instead.
*/
const T &operator()(ssize_t i) const { return data_[i]; }
/*
@brief Index the array as it would be a 1D array. To get a certain
pixel in a multidimensional array use the (i,j,k,...) operator instead.
*/
T &operator[](ssize_t i) { return data_[i]; }
/*
@brief Index the array as it would be a 1D array. To get a certain
pixel in a multidimensional array use the (i,j,k,...) operator instead.
*/
const T &operator[](ssize_t i) const { return data_[i]; }
/* @brief Return a raw pointer to the data */
T *data() { return data_; }
/* @brief Return a const raw pointer to the data */
const T *data() const { return data_; }
/* @brief Return a byte pointer to the data. Useful for memcpy like
* operations */
std::byte *buffer() { return reinterpret_cast<std::byte *>(data_); }
/**
* @brief Return the total number of elements in the array as a signed
* integer
*/
ssize_t size() const { return static_cast<ssize_t>(size_); }
/** @brief Return the total number of bytes in the array */
size_t total_bytes() const { return size_ * sizeof(T); }
/** @brief Return the shape of the array */
Shape<Ndim> shape() const noexcept { return shape_; }
/** @brief Return the size of dimension i */
ssize_t shape(ssize_t i) const noexcept { return shape_[i]; }
/** @brief Return the strides of the array */
std::array<ssize_t, Ndim> strides() const noexcept { return strides_; }
/**
* @brief Return the bitdepth of the array. Useful for checking that
* detector data can fit in the array type.
*/
size_t bitdepth() const noexcept { return sizeof(T) * 8; }
/**
* @brief Return the number of bytes to step in each dimension when
* traversing the array.
*/
std::array<ssize_t, Ndim> byte_strides() const noexcept {
auto byte_strides = strides_;
for (auto &val : byte_strides)
val *= sizeof(T);
return byte_strides;
}
auto begin() const { return data_; }
auto end() const { return data_ + size_; }
using value_type = T;
///////////////////////////////////////////////////////////////////////////////
// Assignments
//
///////////////////////////////////////////////////////////////////////////////
NDArray &operator=(NDArray &&other) noexcept; // Move assign
NDArray &operator=(const NDArray &other); // Copy assign
NDArray &operator+=(const NDArray &other);
NDArray &operator-=(const NDArray &other);
NDArray &operator*=(const NDArray &other);
/**
* @brief Copy to the NDArray from an std::array. If the size of the array
* is different we reallocate the data.
*
*/
// Write directly to the data array, or create a new one
template <size_t Size>
NDArray<T, 1> &operator=(const std::array<T, Size> &other) {
if (Size != size_) {
@@ -276,94 +142,12 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
return *this;
}
/**
* @brief Move assignment operator.
*/
NDArray &operator=(NDArray &&other) noexcept {
// TODO! Should we use swap?
if (this != &other) {
delete[] data_;
data_ = other.data_;
shape_ = other.shape_;
size_ = other.size_;
strides_ = other.strides_;
other.reset();
}
return *this;
}
// NDArray& operator/=(const NDArray& other);
/**
* @brief Copy assignment operator.
*/
NDArray &operator=(const NDArray &other) {
if (this != &other) {
delete[] data_;
shape_ = other.shape_;
strides_ = other.strides_;
size_ = other.size_;
data_ = new T[size_];
std::copy(other.data_, other.data_ + size_, data_);
}
return *this;
}
///////////////////////////////////////////////////////////////////////////////
// Math operators
//
///////////////////////////////////////////////////////////////////////////////
/**
* @brief Add elementwise from another NDArray.
*/
NDArray &operator+=(const NDArray &other) {
if (shape_ != other.shape_)
throw(std::runtime_error(
"Shape of NDArray must match for operator +="));
for (size_t i = 0; i < size_; ++i) {
data_[i] += other.data_[i];
}
return *this;
}
/**
* @brief Subtract elementwise with another NDArray.
*/
NDArray &operator-=(const NDArray &other) {
if (shape_ != other.shape_)
throw(std::runtime_error(
"Shape of NDArray must match for operator -="));
for (size_t i = 0; i < size_; ++i) {
data_[i] -= other.data_[i];
}
return *this;
}
/**
* @brief Multiply elementwise with another NDArray.
*/
NDArray &operator*=(const NDArray &other) {
if (shape_ != other.shape_)
throw(std::runtime_error(
"Shape of NDArray must match for operator *="));
for (size_t i = 0; i < size_; ++i) {
data_[i] *= other.data_[i];
}
return *this;
}
/**
* @brief Divide elementwise by another NDArray. Templated to allow division
* with different types.
*
* TODO! Why is this templated when the others are not?
*/
template <typename V> NDArray &operator/=(const NDArray<V, Ndim> &other) {
// check shape
if (shape_ == other.shape()) {
for (size_t i = 0; i < size_; ++i) {
for (uint32_t i = 0; i < size_; ++i) {
data_[i] /= other(i);
}
return *this;
@@ -371,139 +155,67 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
throw(std::runtime_error("Shape of NDArray must match"));
}
/**
* @brief Assign a scalar value to all elements in the NDArray.
*/
NDArray &operator=(const T &value) {
std::fill_n(data_, size_, value);
return *this;
}
NDArray<bool, Ndim> operator>(const NDArray &other);
/**
* @brief Add a scalar value to all elements in the NDArray.
*/
NDArray &operator+=(const T &value) {
for (size_t i = 0; i < size_; ++i)
data_[i] += value;
return *this;
}
bool operator==(const NDArray &other) const;
bool operator!=(const NDArray &other) const;
/**
* @brief Subtract a scalar value to all elements in the NDArray.
*/
NDArray &operator-=(const T &value) {
for (size_t i = 0; i < size_; ++i)
data_[i] -= value;
return *this;
}
NDArray &operator=(const T & /*value*/);
NDArray &operator+=(const T & /*value*/);
NDArray operator+(const T & /*value*/);
NDArray &operator-=(const T & /*value*/);
NDArray operator-(const T & /*value*/);
NDArray &operator*=(const T & /*value*/);
NDArray operator*(const T & /*value*/);
NDArray &operator/=(const T & /*value*/);
NDArray operator/(const T & /*value*/);
/**
* @brief Multiply all elements in the NDArray with a scalar value
*/
NDArray &operator*=(const T &value) {
for (size_t i = 0; i < size_; ++i)
data_[i] *= value;
return *this;
}
NDArray &operator&=(const T & /*mask*/);
/**
* @brief Divide all elements in the NDArray with a scalar value
*/
NDArray &operator/=(const T &value) {
for (size_t i = 0; i < size_; ++i)
data_[i] /= value;
return *this;
}
/**
* @brief Bitwise AND all elements in the NDArray with a scalar mask.
* Used for example to mask out gain bits for Jungfrau detectors.
*/
NDArray &operator&=(const T &mask) {
for (auto it = begin(); it != end(); ++it)
*it &= mask;
return *this;
}
/**
* @brief Operator + with a scalar value. Returns a new NDArray.
*
* TODO! Expression template version of this?
*/
NDArray operator+(const T &value) {
NDArray result = *this;
result += value;
return result;
}
/**
* @brief Operator - with a scalar value. Returns a new NDArray.
*
* TODO! Expression template version of this?
*/
NDArray operator-(const T &value) {
NDArray result = *this;
result -= value;
return result;
}
/**
* @brief Operator * with a scalar value. Returns a new NDArray.
*
* TODO! Expression template version of this?
*/
NDArray operator*(const T &value) {
NDArray result = *this;
result *= value;
return result;
}
/**
* @brief Operator / with a scalar value. Returns a new NDArray.
*
* TODO! Expression template version of this?
*/
NDArray operator/(const T &value) {
NDArray result = *this;
result /= value;
return result;
}
/**
* @brief Compare two NDArrays elementwise for equality.
*/
bool operator==(const NDArray &other) const {
if (shape_ != other.shape_)
return false;
for (size_t i = 0; i != size_; ++i)
if (data_[i] != other.data_[i])
return false;
return true;
}
/**
* @brief Compare two NDArrays elementwise for non-equality.
*/
bool operator!=(const NDArray &other) const { return !((*this) == other); }
/**
* @brief Compute the square root of all elements in the NDArray.
*/
void sqrt() {
for (size_t i = 0; i < size_; ++i) {
for (int i = 0; i < size_; ++i) {
data_[i] = std::sqrt(data_[i]);
}
}
/*
* @brief Prefix increment operator. Increments all elements by 1.
*/
NDArray &operator++() {
for (size_t i = 0; i < size_; ++i)
data_[i] += T{1};
return *this;
NDArray &operator++(); // pre inc
template <typename... Ix>
std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) {
return data_[element_offset(strides_, index...)];
}
template <typename... Ix>
std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) const {
return data_[element_offset(strides_, index...)];
}
template <typename... Ix>
std::enable_if_t<sizeof...(Ix) == Ndim, T> value(Ix... index) {
return data_[element_offset(strides_, index...)];
}
// TODO! is int the right type for index?
T &operator()(ssize_t i) { return data_[i]; }
const T &operator()(ssize_t i) const { return data_[i]; }
T &operator[](ssize_t i) { return data_[i]; }
const T &operator[](ssize_t i) const { return data_[i]; }
T *data() { return data_; }
std::byte *buffer() { return reinterpret_cast<std::byte *>(data_); }
ssize_t size() const { return static_cast<ssize_t>(size_); }
size_t total_bytes() const { return size_ * sizeof(T); }
std::array<ssize_t, Ndim> shape() const noexcept { return shape_; }
ssize_t shape(ssize_t i) const noexcept { return shape_[i]; }
std::array<ssize_t, Ndim> strides() const noexcept { return strides_; }
size_t bitdepth() const noexcept { return sizeof(T) * 8; }
std::array<ssize_t, Ndim> byte_strides() const noexcept {
auto byte_strides = strides_;
for (auto &val : byte_strides)
val *= sizeof(T);
return byte_strides;
}
/**
@@ -513,12 +225,10 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
*/
NDView<T, Ndim> view() const { return NDView<T, Ndim>{data_, shape_}; }
private:
/**
* @brief Reset the NDArray to an empty state. Dropping the ownership of
* the data. Used internally for move operations to avoid double free or
* dangling pointers.
*/
void Print();
void Print_all();
void Print_some();
void reset() {
data_ = nullptr;
size_ = 0;
@@ -527,10 +237,167 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
}
};
///////////////////////////////////////////////////////////////////////////////
// Free functions closely related to NDArray
//
///////////////////////////////////////////////////////////////////////////////
// Move assign
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &
NDArray<T, Ndim>::operator=(NDArray<T, Ndim> &&other) noexcept {
if (this != &other) {
delete[] data_;
data_ = other.data_;
shape_ = other.shape_;
size_ = other.size_;
strides_ = other.strides_;
other.reset();
}
return *this;
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const NDArray<T, Ndim> &other) {
// check shape
if (shape_ == other.shape_) {
for (size_t i = 0; i < size_; ++i) {
data_[i] += other.data_[i];
}
return *this;
}
throw(std::runtime_error("Shape of ImageDatas must match"));
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator-=(const NDArray<T, Ndim> &other) {
// check shape
if (shape_ == other.shape_) {
for (uint32_t i = 0; i < size_; ++i) {
data_[i] -= other.data_[i];
}
return *this;
}
throw(std::runtime_error("Shape of ImageDatas must match"));
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator*=(const NDArray<T, Ndim> &other) {
// check shape
if (shape_ == other.shape_) {
for (uint32_t i = 0; i < size_; ++i) {
data_[i] *= other.data_[i];
}
return *this;
}
throw(std::runtime_error("Shape of ImageDatas must match"));
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator&=(const T &mask) {
for (auto it = begin(); it != end(); ++it)
*it &= mask;
return *this;
}
template <typename T, ssize_t Ndim>
NDArray<bool, Ndim> NDArray<T, Ndim>::operator>(const NDArray &other) {
if (shape_ == other.shape_) {
NDArray<bool, Ndim> result{shape_};
for (int i = 0; i < size_; ++i) {
result(i) = (data_[i] > other.data_[i]);
}
return result;
}
throw(std::runtime_error("Shape of ImageDatas must match"));
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator=(const NDArray<T, Ndim> &other) {
if (this != &other) {
delete[] data_;
shape_ = other.shape_;
strides_ = other.strides_;
size_ = other.size_;
data_ = new T[size_];
std::copy(other.data_, other.data_ + size_, data_);
}
return *this;
}
template <typename T, ssize_t Ndim>
bool NDArray<T, Ndim>::operator==(const NDArray<T, Ndim> &other) const {
if (shape_ != other.shape_)
return false;
for (uint32_t i = 0; i != size_; ++i)
if (data_[i] != other.data_[i])
return false;
return true;
}
template <typename T, ssize_t Ndim>
bool NDArray<T, Ndim>::operator!=(const NDArray<T, Ndim> &other) const {
return !((*this) == other);
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator++() {
for (uint32_t i = 0; i < size_; ++i)
data_[i] += 1;
return *this;
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator=(const T &value) {
std::fill_n(data_, size_, value);
return *this;
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const T &value) {
for (uint32_t i = 0; i < size_; ++i)
data_[i] += value;
return *this;
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator+(const T &value) {
NDArray result = *this;
result += value;
return result;
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator-=(const T &value) {
for (uint32_t i = 0; i < size_; ++i)
data_[i] -= value;
return *this;
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator-(const T &value) {
NDArray result = *this;
result -= value;
return result;
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator/=(const T &value) {
for (uint32_t i = 0; i < size_; ++i)
data_[i] /= value;
return *this;
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator/(const T &value) {
NDArray result = *this;
result /= value;
return result;
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator*=(const T &value) {
for (uint32_t i = 0; i < size_; ++i)
data_[i] *= value;
return *this;
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator*(const T &value) {
NDArray result = *this;
result *= value;
return result;
}
template <typename T, ssize_t Ndim>
std::ostream &operator<<(std::ostream &os, const NDArray<T, Ndim> &arr) {
@@ -544,9 +411,27 @@ std::ostream &operator<<(std::ostream &os, const NDArray<T, Ndim> &arr) {
return os;
}
template <typename T, ssize_t Ndim> void NDArray<T, Ndim>::Print_all() {
for (auto row = 0; row < shape_[0]; ++row) {
for (auto col = 0; col < shape_[1]; ++col) {
std::cout << std::setw(3);
std::cout << (*this)(row, col) << " ";
}
std::cout << "\n";
}
}
template <typename T, ssize_t Ndim> void NDArray<T, Ndim>::Print_some() {
for (auto row = 0; row < 5; ++row) {
for (auto col = 0; col < 5; ++col) {
std::cout << std::setw(7);
std::cout << (*this)(row, col) << " ";
}
std::cout << "\n";
}
}
template <typename T, ssize_t Ndim>
[[deprecated("Saving of raw arrays without metadata is deprecated")]] void
save(NDArray<T, Ndim> &img, std::string &pathname) {
void save(NDArray<T, Ndim> &img, std::string &pathname) {
std::ofstream f;
f.open(pathname, std::ios::binary);
f.write(img.buffer(), img.size() * sizeof(T));
@@ -554,9 +439,8 @@ save(NDArray<T, Ndim> &img, std::string &pathname) {
}
template <typename T, ssize_t Ndim>
[[deprecated(
"Loading of raw arrays without metadata is deprecated")]] NDArray<T, Ndim>
load(const std::string &pathname, std::array<ssize_t, Ndim> shape) {
NDArray<T, Ndim> load(const std::string &pathname,
std::array<ssize_t, Ndim> shape) {
NDArray<T, Ndim> img{shape};
std::ifstream f;
f.open(pathname, std::ios::binary);
@@ -565,20 +449,6 @@ load(const std::string &pathname, std::array<ssize_t, Ndim> shape) {
return img;
}
/**
* @brief Free function to safely divide two NDArrays elementwise, handling
* division by zero. Uses static_cast to convert types as needed.
*
* @tparam RT Result type
* @tparam NT Numerator type
* @tparam DT Denominator type
* @tparam Ndim Number of dimensions
* @param numerator The numerator NDArray
* @param denominator The denominator NDArray
* @return NDArray<RT, Ndim> Resulting NDArray after safe division
* @throws std::runtime_error if the shapes of the numerator and denominator do
* not match
*/
template <typename RT, typename NT, typename DT, ssize_t Ndim>
NDArray<RT, Ndim> safe_divide(const NDArray<NT, Ndim> &numerator,
const NDArray<DT, Ndim> &denominator) {

View File

@@ -6,7 +6,6 @@
#include <fmt/format.h>
#include <fstream>
#include <optional>
#include <chrono>
#include <nlohmann/json.hpp>
using json = nlohmann::json;
@@ -85,9 +84,6 @@ class RawMasterFile {
size_t m_bitdepth{};
uint8_t m_quad = 0;
std::optional<std::chrono::nanoseconds> m_exptime;
std::chrono::nanoseconds m_period{0};
xy m_geometry{};
xy m_udp_interfaces_per_module{1, 1};
@@ -113,7 +109,6 @@ class RawMasterFile {
public:
RawMasterFile(const std::filesystem::path &fpath);
RawMasterFile(std::istream &is, const std::string &fname); // for testing
std::filesystem::path data_fname(size_t mod_id, size_t file_id) const;
@@ -145,12 +140,9 @@ class RawMasterFile {
ScanParameters scan_parameters() const;
std::optional<std::chrono::nanoseconds> exptime() const { return m_exptime; }
std::chrono::nanoseconds period() const { return m_period; }
private:
void parse_json(std::istream &is);
void parse_raw(std::istream &is);
void parse_json(const std::filesystem::path &fpath);
void parse_raw(const std::filesystem::path &fpath);
void retrieve_geometry();
};

View File

@@ -125,7 +125,7 @@ template <typename T> int VarClusterFinder<T>::check_neighbours(int i, int j) {
const auto row = i + di[k];
const auto col = j + dj[k];
if (row >= 0 && col >= 0 && row < shape_[0] && col < shape_[1]) {
auto tmp = labeled_(i + di[k], j + dj[k]);
auto tmp = labeled_.value(i + di[k], j + dj[k]);
if (tmp != 0)
neighbour_labels.push_back(tmp);
}

View File

@@ -1,12 +1,11 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/defs.hpp"
#include <aare/NDView.hpp>
#include <cstdint>
#include <vector>
namespace aare {
uint16_t adc_sar_05_decode64to16(uint64_t input);
uint16_t adc_sar_04_decode64to16(uint64_t input);
void adc_sar_05_decode64to16(NDView<uint64_t, 2> input,
@@ -14,25 +13,6 @@ void adc_sar_05_decode64to16(NDView<uint64_t, 2> input,
void adc_sar_04_decode64to16(NDView<uint64_t, 2> input,
NDView<uint16_t, 2> output);
/**
* @brief Called with a 32 bit unsigned integer, shift by offset
* and then return the lower 24 bits as an 32 bit integer
* @param input 32-ibt input value
* @param offset (should be in range 0-7 to allow for full 24 bits)
* @return uint32_t
*/
uint32_t mask32to24bits(uint32_t input, BitOffset offset={});
/**
* @brief Expand 24 bit values in a 8bit buffer to 32bit unsigned integers
* Used for detectors with 24bit counters in combination with CTB
*
* @param input View of the 24 bit data as uint8_t (no 24bit native data type exists)
* @param output Destination of the expanded data (32bit, unsigned)
* @param offset Offset within the first byte to where the data starts (0-7 bits)
*/
void expand24to32bit(NDView<uint8_t,1> input, NDView<uint32_t,1> output, BitOffset offset={});
/**
* @brief Apply custom weights to a 16-bit input value. Will sum up
* weights[i]**i for each bit i that is set in the input value.

View File

@@ -348,21 +348,28 @@ enum class corner : int {
enum class TimingMode { Auto, Trigger };
enum class FrameDiscardPolicy { NoDiscard, Discard, DiscardPartial };
template <class T> T StringTo(const std::string &arg) { return T(arg); }
template <class T> std::string ToString(T arg) { return T(arg); }
template <> DetectorType StringTo(const std::string & /*name*/);
template <> std::string ToString(DetectorType arg);
template <> TimingMode StringTo(const std::string & /*mode*/);
template <> FrameDiscardPolicy StringTo(const std::string & /*mode*/);
using DataTypeVariants = std::variant<uint16_t, uint32_t>;
constexpr uint16_t ADC_MASK =
0x3FFF; // used to mask out the gain bits in Jungfrau
class BitOffset{
uint8_t m_offset{};
public:
BitOffset() = default;
explicit BitOffset(uint32_t offset);
uint8_t value() const {return m_offset;}
bool operator==(const BitOffset& other) const;
bool operator<(const BitOffset& other) const;
};
/**
* @brief Convert a string to a DACIndex
* @param arg string representation of the dacIndex
* @return DACIndex
* @throw invalid argument error if the string does not match any DACIndex
*/
template <> DACIndex StringTo(const std::string &arg);
} // namespace aare

View File

@@ -106,7 +106,7 @@ class Logger {
}
std::ostringstream &Get() {
os << Color(m_level) << "- " << Timestamp() << " " << Logger::ToString(m_level)
os << Color(m_level) << "- " << Timestamp() << " " << ToString(m_level)
<< ": ";
return os;
}

View File

@@ -39,6 +39,8 @@ set( PYTHON_FILES
aare/transform.py
aare/ScanParameters.py
aare/utils.py
aare/experimental/__init__.py
aare/experimental/custom_io.py
)

View File

View File

@@ -0,0 +1,58 @@
import numpy as np
n_counters = 64*3
bitfield_size = 64
header_dt = [('frameNumber',np.uint64),
('expLength',np.uint32),
('packetNumber', np.uint32),
('bunchId', np.uint64),
('timestamp', np.uint64),
('modId', np.uint16),
('row', np.uint16),
('col', np.uint16),
('reserved', np.uint16),
('debug', np.uint32),
('roundRNumber', np.uint16),
('detType', np.uint8),
('version', np.uint8)]
def ExtractBits(raw_data, dr=24, bits = (17,6)):
bits = np.uint64(bits)
data = np.zeros(0, dtype = np.uint64)
for bit in bits:
tmp = (raw_data >> bit) & np.uint64(1)
data = np.hstack((data, tmp))
#Shift the bits to the righ place
for i in np.arange(dr, dtype = np.uint64):
data[i::dr] = data[i::dr] << i
data = data.reshape(data.size//dr, dr)
return data.sum(axis = 1)
def read_my302_file(fname, dr=24, bits = (17,6),
offset=48, tail = 72, n_frames=1):
header = np.zeros(n_frames, header_dt)
data = np.zeros((n_frames, n_counters), dtype = np.uint64)
with open(fname, 'rb') as f:
for i in range(n_frames):
header[i], raw_data = _read_my302_frame(f, offset, tail, dr)
data[i] = ExtractBits(raw_data, dr=dr, bits = bits)
return header, data
def _read_my302_frame(f, offset, tail, dr):
header = np.fromfile(f, count=1, dtype = header_dt)
f.seek(bitfield_size+offset, 1)
data = np.fromfile(f, count = int(n_counters*dr/2), dtype = np.uint64)
f.seek(tail, 1)
return header, data

View File

@@ -49,43 +49,6 @@ class Matterhorn02Transform:
else:
return np.take(data.view(np.uint16), self.pixel_map[0:counters])
class Mythen302Transform:
"""
Transform Mythen 302 test chip data from a buffer of bytes (uint8_t)
to a uint32 numpy array of [64,3] representing channels and counters.
Assumes data taken with rx_dbitlist 17 6, rx_dbitreorder 1 and Digital
Samples = 2310 [(64x3x24)/2 + some extra]
.. note::
The offset is in number of bits 0-7
"""
_n_channels = 64
_n_counters = 3
def __init__(self, offset=4):
self.offset = offset
def __call__(self, data : np.ndarray):
"""
Transform buffer of data to a [64,3] np.ndarray of uint32.
Parameters
----------
data : np.ndarray
Expected dtype: uint8
Returns
----------
image : np.ndarray
uint32 array of size 64, 3
"""
res = _aare.decode_my302(data, self.offset)
res = res.reshape(
Mythen302Transform._n_channels, Mythen302Transform._n_counters
)
return res
#on import generate the pixel maps to avoid doing it every time
moench05 = Moench05Transform()

View File

@@ -96,69 +96,6 @@ void define_ctb_raw_file_io_bindings(py::module &m) {
return output;
});
m.def("expand24to32bit",
[](py::array_t<uint8_t, py::array::c_style | py::array::forcecast>
&input, uint32_t offset){
aare::BitOffset bitoff(offset);
py::buffer_info buf = input.request();
constexpr uint32_t bytes_per_channel = 3; //24 bit
py::array_t<uint32_t> output(buf.size/bytes_per_channel);
NDView<uint8_t, 1> input_view(input.mutable_data(),
{input.size()});
NDView<uint32_t, 1> output_view(output.mutable_data(),
{output.size()});
aare::expand24to32bit(input_view, output_view, bitoff);
return output;
});
m.def("decode_my302",
[](py::array_t<uint8_t, py::array::c_style | py::array::forcecast>
&input, uint32_t offset){
// Physical layout of the chip
constexpr size_t channels = 64;
constexpr size_t counters = 3;
constexpr size_t bytes_per_channel = 3; //24 bit
constexpr int n_outputs = 2;
ssize_t expected_size = channels*counters*bytes_per_channel;
//If whe have an offset we need one extra byte per output
aare::BitOffset bitoff(offset);
if(bitoff.value())
expected_size += n_outputs;
if (input.size() != expected_size) {
throw std::runtime_error(
fmt::format("{} Expected an input size of {} bytes. Called "
"with input size of {}",
LOCATION, expected_size, input.size()));
}
py::buffer_info buf = input.request();
py::array_t<uint32_t> output(channels * counters);
for (int i = 0; i!=n_outputs; ++i){
auto step = input.size()/n_outputs;
auto out_step = output.size()/n_outputs;
NDView<uint8_t, 1> input_view(input.mutable_data()+step*i,
{input.size()/n_outputs});
NDView<uint32_t, 1> output_view(output.mutable_data()+out_step*i,
{output.size()/n_outputs});
aare::expand24to32bit(input_view, output_view, bitoff);
}
return output;
});
py::class_<CtbRawFile>(m, "CtbRawFile")
.def(py::init<const std::filesystem::path &>())
.def("read_frame",

View File

@@ -66,7 +66,7 @@ void define_file_io_bindings(py::module &m) {
.def_property_readonly("bytes_per_pixel", &File::bytes_per_pixel)
.def_property_readonly(
"detector_type",
[](File &self) { return self.detector_type(); })
[](File &self) { return ToString(self.detector_type()); })
.def("read_frame",
[](File &self) {
const uint8_t item_size = self.bytes_per_pixel();
@@ -161,6 +161,21 @@ void define_file_io_bindings(py::module &m) {
}
});
py::class_<FileConfig>(m, "FileConfig")
.def(py::init<>())
.def_readwrite("rows", &FileConfig::rows)
.def_readwrite("cols", &FileConfig::cols)
.def_readwrite("version", &FileConfig::version)
.def_readwrite("geometry", &FileConfig::geometry)
.def_readwrite("detector_type", &FileConfig::detector_type)
.def_readwrite("max_frames_per_file", &FileConfig::max_frames_per_file)
.def_readwrite("total_frames", &FileConfig::total_frames)
.def_readwrite("dtype", &FileConfig::dtype)
.def("__eq__", &FileConfig::operator==)
.def("__ne__", &FileConfig::operator!=)
.def("__repr__", [](const FileConfig &a) {
return "<FileConfig: " + a.to_string() + ">";
});
py::class_<ScanParameters>(m, "ScanParameters")
.def(py::init<const std::string &>())

View File

@@ -48,19 +48,6 @@ void register_interpolate(py::class_<aare::Interpolator> &interpolator,
docstring.c_str(), py::arg("cluster_vector"));
}
template <typename Type>
void register_transform_eta_values(
py::class_<aare::Interpolator> &interpolator) {
interpolator.def(
"transform_eta_values",
[](Interpolator &self, const Eta2<Type> &eta) {
auto uniform_coord = self.transform_eta_values(eta);
return py::make_tuple(uniform_coord.x, uniform_coord.y);
},
R"(eta.x eta.y transformed to uniform coordinates based on CDF ietax, ietay)",
py::arg("Eta"));
}
void define_interpolation_bindings(py::module &m) {
PYBIND11_NUMPY_DTYPE(aare::Photon, x, y, energy);
@@ -153,10 +140,6 @@ void define_interpolation_bindings(py::module &m) {
REGISTER_INTERPOLATOR_ETA2(float, 2, 2, uint16_t);
REGISTER_INTERPOLATOR_ETA2(double, 2, 2, uint16_t);
register_transform_eta_values<int>(interpolator);
register_transform_eta_values<float>(interpolator);
register_transform_eta_values<double>(interpolator);
// TODO! Evaluate without converting to double
m.def(
"hej",

View File

@@ -10,13 +10,13 @@
#include "bind_ClusterFinderMT.hpp"
#include "bind_ClusterVector.hpp"
#include "bind_Eta.hpp"
#include "bind_Interpolator.hpp"
#include "bind_calibration.hpp"
// TODO! migrate the other names
#include "ctb_raw_file.hpp"
#include "file.hpp"
#include "fit.hpp"
#include "interpolation.hpp"
#include "jungfrau_data_file.hpp"
#include "pedestal.hpp"
#include "pixel_map.hpp"

View File

@@ -85,21 +85,5 @@ void define_raw_master_file_bindings(py::module &m) {
.def_property_readonly("quad", &RawMasterFile::quad)
.def_property_readonly("scan_parameters",
&RawMasterFile::scan_parameters)
.def_property_readonly("roi", &RawMasterFile::roi)
.def_property_readonly(
"exptime",
[](RawMasterFile &self) -> std::optional<double> {
if (self.exptime()) {
double seconds =
std::chrono::duration<double>(*self.exptime()).count();
return seconds;
} else {
return std::nullopt;
}
})
.def_property_readonly("period", [](RawMasterFile &self) {
double seconds =
std::chrono::duration<double>(self.period()).count();
return seconds;
});
.def_property_readonly("roi", &RawMasterFile::roi);
}

View File

@@ -44,18 +44,13 @@ def test_Interpolator():
etacube = np.zeros(shape=[29, 29, 19], dtype=np.float64)
interpolator = _aare.Interpolator(etacube, xbins, ybins, ebins)
assert interpolator.get_ietax().shape == (29,29,19)
assert interpolator.get_ietay().shape == (29,29,19)
assert interpolator.get_ietax().shape == (30,30,20)
assert interpolator.get_ietay().shape == (30,30,20)
clustervector = _aare.ClusterVector_Cluster3x3i()
cluster = _aare.Cluster3x3i(1,1, np.ones(9, dtype=np.int32))
clustervector.push_back(cluster)
[u,v] = interpolator.transform_eta_values(_aare.Etai())
assert u == 0
assert v == 0
interpolated_photons = interpolator.interpolate(clustervector)
assert interpolated_photons.size == 1

View File

@@ -0,0 +1,146 @@
# SPDX-License-Identifier: MPL-2.0
import pytest
import numpy as np
import boost_histogram as bh
import pickle
from scipy.stats import multivariate_normal
from aare import Interpolator, calculate_eta2
from aare._aare import ClusterVector_Cluster2x2d, Cluster2x2d, Cluster3x3d, ClusterVector_Cluster3x3d
from conftest import test_data_path
pixel_width = 1e-4
values = np.arange(0.5*pixel_width, 0.1, pixel_width)
num_pixels = values.size
X, Y = np.meshgrid(values, values)
data_points = np.stack([X.ravel(), Y.ravel()], axis=1)
variance = 10*pixel_width
covariance_matrix = np.array([[variance, 0],[0, variance]])
def create_photon_hit_with_gaussian_distribution(mean, covariance_matrix, data_points):
gaussian = multivariate_normal(mean=mean, cov=covariance_matrix)
probability_values = gaussian.pdf(data_points)
return (probability_values.reshape(X.shape)).round() #python bindings only support frame types of uint16_t
def create_2x2cluster_from_frame(frame, pixels_per_superpixel):
return Cluster2x2d(1, 1, np.array([frame[0:pixels_per_superpixel, 0:pixels_per_superpixel].sum(),
frame[0:pixels_per_superpixel, pixels_per_superpixel:2*pixels_per_superpixel].sum(),
frame[pixels_per_superpixel:2*pixels_per_superpixel, 0:pixels_per_superpixel].sum(),
frame[pixels_per_superpixel:2*pixels_per_superpixel, pixels_per_superpixel:2*pixels_per_superpixel].sum()], dtype=np.float64))
def create_3x3cluster_from_frame(frame, pixels_per_superpixel):
return Cluster3x3d(1, 1, np.array([frame[0:pixels_per_superpixel, 0:pixels_per_superpixel].sum(),
frame[0:pixels_per_superpixel, pixels_per_superpixel:2*pixels_per_superpixel].sum(),
frame[0:pixels_per_superpixel, 2*pixels_per_superpixel:3*pixels_per_superpixel].sum(),
frame[pixels_per_superpixel:2*pixels_per_superpixel, 0:pixels_per_superpixel].sum(),
frame[pixels_per_superpixel:2*pixels_per_superpixel, pixels_per_superpixel:2*pixels_per_superpixel].sum(),
frame[pixels_per_superpixel:2*pixels_per_superpixel, 2*pixels_per_superpixel:3*pixels_per_superpixel].sum(),
frame[2*pixels_per_superpixel:3*pixels_per_superpixel, 0:pixels_per_superpixel].sum(),
frame[2*pixels_per_superpixel:3*pixels_per_superpixel, pixels_per_superpixel:2*pixels_per_superpixel].sum(),
frame[2*pixels_per_superpixel:3*pixels_per_superpixel, 2*pixels_per_superpixel:3*pixels_per_superpixel].sum()], dtype=np.float64))
def calculate_eta_distribution(num_frames, pixels_per_superpixel, random_number_generator, bin_edges_x = bh.axis.Regular(100, -0.2, 1.2), bin_edges_y = bh.axis.Regular(100, -0.2, 1.2), cluster_2x2 = True):
hist = bh.Histogram(
bin_edges_x,
bin_edges_y, bh.axis.Regular(1, 0, num_pixels*num_pixels*1/(variance*2*np.pi)))
for _ in range(0, num_frames):
mean_x = random_number_generator.uniform(pixels_per_superpixel*pixel_width, 2*pixels_per_superpixel*pixel_width)
mean_y = random_number_generator.uniform(pixels_per_superpixel*pixel_width, 2*pixels_per_superpixel*pixel_width)
frame = create_photon_hit_with_gaussian_distribution(np.array([mean_x, mean_y]), variance, data_points)
cluster = None
if cluster_2x2:
cluster = create_2x2cluster_from_frame(frame, pixels_per_superpixel)
else:
cluster = create_3x3cluster_from_frame(frame, pixels_per_superpixel)
eta2 = calculate_eta2(cluster)
hist.fill(eta2.x, eta2.y, eta2.sum)
return hist
@pytest.mark.withdata
def test_interpolation_of_2x2_cluster(test_data_path):
"""Test Interpolation of 2x2 cluster from Photon hit with Gaussian Distribution"""
#TODO maybe better to compute in test instead of loading - depends on eta
"""
filename = test_data_path/"eta_distributions"/"eta_distribution_2x2cluster_gaussian.pkl"
with open(filename, "rb") as f:
eta_distribution = pickle.load(f)
"""
num_frames = 1000
pixels_per_superpixel = int(num_pixels*0.5)
random_number_generator = np.random.default_rng(42)
eta_distribution = calculate_eta_distribution(num_frames, pixels_per_superpixel, random_number_generator, bin_edges_x = bh.axis.Regular(100, -0.1, 0.6), bin_edges_y = bh.axis.Regular(100, -0.1, 0.6))
interpolation = Interpolator(eta_distribution, eta_distribution.axes[0].edges, eta_distribution.axes[1].edges, eta_distribution.axes[2].edges)
#actual photon hit
mean = 1.2*pixels_per_superpixel*pixel_width
mean = np.array([mean, mean])
frame = create_photon_hit_with_gaussian_distribution(mean, covariance_matrix, data_points)
cluster = create_2x2cluster_from_frame(frame, pixels_per_superpixel)
clustervec = ClusterVector_Cluster2x2d()
clustervec.push_back(cluster)
interpolated_photon = interpolation.interpolate(clustervec)
assert interpolated_photon.size == 1
cluster_center = 1.5*pixels_per_superpixel*pixel_width
scaled_photon_hit = (interpolated_photon[0][0]*pixels_per_superpixel*pixel_width, interpolated_photon[0][1]*pixels_per_superpixel*pixel_width)
assert (np.linalg.norm(scaled_photon_hit - mean) < np.linalg.norm(np.array([cluster_center, cluster_center] - mean)))
@pytest.mark.withdata
def test_interpolation_of_3x3_cluster(test_data_path):
"""Test Interpolation of 3x3 Cluster from Photon hit with Gaussian Distribution"""
#TODO maybe better to compute in test instead of loading - depends on eta
"""
filename = test_data_path/"eta_distributions"/"eta_distribution_3x3cluster_gaussian.pkl"
with open(filename, "rb") as f:
eta_distribution = pickle.load(f)
"""
num_frames = 1000
pixels_per_superpixel = int(num_pixels/3)
random_number_generator = np.random.default_rng(42)
eta_distribution = calculate_eta_distribution(num_frames, pixels_per_superpixel, random_number_generator, bin_edges_x = bh.axis.Regular(100, -0.1, 1.1), bin_edges_y = bh.axis.Regular(100, -0.1, 1.1), cluster_2x2 = False)
interpolation = Interpolator(eta_distribution, eta_distribution.axes[0].edges, eta_distribution.axes[1].edges, eta_distribution.axes[2].edges)
#actual photon hit
mean_x = (1 + 0.8)*pixels_per_superpixel*pixel_width
mean_y = (1 + 0.2)*pixels_per_superpixel*pixel_width
mean = np.array([mean_x, mean_y])
frame = create_photon_hit_with_gaussian_distribution(mean, covariance_matrix, data_points)
cluster = create_3x3cluster_from_frame(frame, pixels_per_superpixel)
clustervec = ClusterVector_Cluster3x3d()
clustervec.push_back(cluster)
interpolated_photon = interpolation.interpolate(clustervec)
assert interpolated_photon.size == 1
cluster_center = 1.5*pixels_per_superpixel*pixel_width
scaled_photon_hit = (interpolated_photon[0][0]*pixels_per_superpixel*pixel_width, interpolated_photon[0][1]*pixels_per_superpixel*pixel_width)
assert (np.linalg.norm(scaled_photon_hit - mean) < np.linalg.norm(np.array([cluster_center, cluster_center] - mean)))

View File

@@ -5,6 +5,8 @@ import pytest
import pytest_check as check
import numpy as np
import boost_histogram as bh
import pickle
from scipy.stats import multivariate_normal
from aare import Interpolator, calculate_eta2, calculate_cross_eta3, calculate_full_eta2, calculate_eta3
from aare import ClusterFile

View File

@@ -61,6 +61,7 @@ void CtbRawFile::find_subfiles() {
while (std::filesystem::exists(m_master.data_fname(0, m_num_subfiles)))
m_num_subfiles++;
fmt::print("Found {} subfiles\n", m_num_subfiles);
}
void CtbRawFile::open_data_file(size_t subfile_index) {

View File

@@ -11,9 +11,9 @@ using namespace aare;
TEST_CASE("Test new Interpolation API", "[Interpolation]") {
NDArray<double, 1> energy_bins(std::array<double, 2>{0.0, 100.0});
NDArray<double, 1> etax_bins(std::array<double, 4>{0.0, 0.3, 0.6, 1.0});
NDArray<double, 1> etay_bins(std::array<double, 4>{0.0, 0.3, 0.6, 1.0});
NDArray<double, 1> energy_bins(std::array<ssize_t, 1>{2});
NDArray<double, 1> etax_bins(std::array<ssize_t, 1>{4}, 0.0);
NDArray<double, 1> etay_bins(std::array<ssize_t, 1>{4}, 0.0);
NDArray<double, 3> eta_distribution(std::array<ssize_t, 3>{3, 3, 1}, 0.0);
Interpolator interpolator(eta_distribution.view(), etax_bins.view(),

View File

@@ -20,7 +20,7 @@ Interpolator::Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
m_ietax = NDArray<double, 3>(etacube);
m_ietay = NDArray<double, 3>(etacube);
// TODO: etacube should have different strides energy should come first
// prefix sum - conditional CDF
for (ssize_t i = 0; i < m_ietax.shape(0); i++) {
for (ssize_t j = 0; j < m_ietax.shape(1); j++) {

View File

@@ -105,67 +105,6 @@ TEST_CASE("Indexing of a 3D image") {
REQUIRE(img(2, 3, 1) == 23);
}
TEST_CASE("Access to data using a pointer"){
// This pattern is discouraged but sometimes useful
NDArray<int,2> img{{4,5},0};
int* data_ptr = img.data();
for(int i=0; i < img.size(); ++i){
data_ptr[i] = i*2;
}
// Cross check using operator[]
for(int i=0; i < img.size(); ++i){
REQUIRE(img[i] == i*2);
}
}
TEST_CASE("Access to data using a pointer for a const NDArray"){
// This pattern is discouraged but sometimes useful
// Using a lambda to create a const NDArray with known data
const NDArray<int,2> arr = [](){
NDArray<int,2> img{{4,5},0};
int* data_ptr = img.data();
for(int i=0; i < img.size(); ++i){
data_ptr[i] = i*3;
}
return img;
}();
// Cross check using data() pointer, if compiles we can get a const pointer
const int* const_data_ptr = arr.data();
for(int i=0; i < arr.size(); ++i){
REQUIRE(const_data_ptr[i] == i*3);
}
}
TEST_CASE("Use *buffer"){
// Another useful but discouraged pattern. But can be useful when getting data
// from external sources
Shape<2> shape{{4, 5}};
NDArray<int, 2> src(shape);
NDArray<int, 2> dst(shape);
for (uint32_t i = 0; i < src.size(); ++i) {
src(i) = static_cast<int>(i * 7);
}
std::memcpy(dst.buffer(), src.buffer(), src.total_bytes());
REQUIRE(src.data() != dst.data());
for (uint32_t i = 0; i < dst.size(); ++i) {
REQUIRE(dst(i) == src(i));
}
}
TEST_CASE("Increment elements using prefix ++ operator") {
NDArray<int, 1> a{{5}, 0};
++a;
for (const auto it : a) {
REQUIRE(it == 1);
}
}
TEST_CASE("Divide double by int") {
NDArray<double, 1> a{{5}, 5};
NDArray<int, 1> b{{5}, 5};
@@ -491,42 +430,6 @@ TEST_CASE("Construct an NDArray from an std::array") {
}
}
TEST_CASE("Copy construct an NDArray"){
NDArray<int,2> a({{3,4}},0);
a(1,1) = 42;
a(2,3) = 84;
NDArray<int,2> b(a);
REQUIRE(b.shape() == Shape<2>{3,4});
REQUIRE(b.size() == 12);
REQUIRE(b(1,1) == 42);
REQUIRE(b(2,3) == 84);
// Modifying b should not affect a
b(1,1) = 7;
REQUIRE(a(1,1) == 42);
REQUIRE(a.data() != b.data());
}
TEST_CASE("Move construct an NDArray"){
NDArray<int,2> a({{3,4}},0);
a(1,1) = 42;
a(2,3) = 84;
NDArray<int,2> b(std::move(a));
REQUIRE(b.shape() == Shape<2>{3,4});
REQUIRE(b.size() == 12);
REQUIRE(b(1,1) == 42);
REQUIRE(b(2,3) == 84);
// The moved from object should be in a unspecified but valid state.
// This means original array pointer should be null, and size zero
REQUIRE(a.size() == 0);
REQUIRE(a.shape() == Shape<2>{0,0});
REQUIRE(a.data() == nullptr);
}
TEST_CASE("Move construct from an array with Ndim + 1") {

View File

@@ -5,7 +5,6 @@
#include <iostream>
#include <numeric>
#include <vector>
#include <cstddef>
using aare::NDView;
using aare::Shape;
@@ -261,11 +260,3 @@ TEST_CASE("Create a view over a vector") {
REQUIRE(v[0] == 0);
REQUIRE(v[11] == 11);
}
TEST_CASE("NDView over byte"){
std::vector<std::byte> buf(5);
auto v = aare::make_view(buf);
REQUIRE(v.shape()[0] == 5);
REQUIRE(v[0] == std::byte{0});
}

View File

@@ -4,8 +4,6 @@
#include "aare/logger.hpp"
#include <sstream>
#include "to_string.hpp"
namespace aare {
RawFileNameComponents::RawFileNameComponents(
@@ -71,7 +69,7 @@ ScanParameters::ScanParameters(const bool enabled, const DACIndex dac,
const int start, const int stop, const int step,
const int64_t settleTime)
: m_enabled(enabled), m_dac(dac), m_start(start), m_stop(stop),
m_step(step), m_settleTime(settleTime) {};
m_step(step), m_settleTime(settleTime){};
// "[enabled\ndac dac 4\nstart 500\nstop 2200\nstep 5\nsettleTime 100us\n]"
ScanParameters::ScanParameters(const std::string &par) {
@@ -81,7 +79,7 @@ ScanParameters::ScanParameters(const std::string &par) {
if (line == "enabled") {
m_enabled = true;
} else if (line.find("dac") != std::string::npos) {
m_dac = string_to<DACIndex>(line.substr(4));
m_dac = StringTo<DACIndex>(line.substr(4));
} else if (line.find("start") != std::string::npos) {
m_start = std::stoi(line.substr(6));
} else if (line.find("stop") != std::string::npos) {
@@ -106,24 +104,10 @@ RawMasterFile::RawMasterFile(const std::filesystem::path &fpath)
throw std::runtime_error(fmt::format("{} File does not exist: {}",
LOCATION, fpath.string()));
}
std::ifstream ifs(fpath);
if (m_fnc.ext() == ".json") {
parse_json(ifs);
parse_json(fpath);
} else if (m_fnc.ext() == ".raw") {
parse_raw(ifs);
} else {
throw std::runtime_error(LOCATION + "Unsupported file type");
}
}
RawMasterFile::RawMasterFile(std::istream &is, const std::string &fname)
: m_fnc(fname) {
if (m_fnc.ext() == ".json") {
parse_json(is);
} else if (m_fnc.ext() == ".raw") {
parse_raw(is);
parse_raw(fpath);
} else {
throw std::runtime_error(LOCATION + "Unsupported file type");
}
@@ -195,15 +179,16 @@ ScanParameters RawMasterFile::scan_parameters() const {
std::optional<ROI> RawMasterFile::roi() const { return m_roi; }
void RawMasterFile::parse_json(std::istream &is) {
void RawMasterFile::parse_json(const std::filesystem::path &fpath) {
std::ifstream ifs(fpath);
json j;
is >> j;
ifs >> j;
double v = j["Version"];
m_version = fmt::format("{:.1f}", v);
m_type = string_to<DetectorType>(j["Detector Type"].get<std::string>());
m_timing_mode = string_to<TimingMode>(j["Timing Mode"].get<std::string>());
m_type = StringTo<DetectorType>(j["Detector Type"].get<std::string>());
m_timing_mode = StringTo<TimingMode>(j["Timing Mode"].get<std::string>());
m_geometry = {
j["Geometry"]["y"],
@@ -219,46 +204,24 @@ void RawMasterFile::parse_json(std::istream &is) {
m_max_frames_per_file = j["Max Frames Per File"];
// Before v8.0 we had Exptime instead of Exposure Time
// Mythen3 uses 3 exposure times and is not handled at the moment
if (j.contains("Exptime") && j["Exptime"].is_string()) {
m_exptime = string_to<std::chrono::nanoseconds>(
j["Exptime"].get<std::string>());
}
if (j.contains("Exposure Time") && j["Exposure Time"].is_string()) {
m_exptime = string_to<std::chrono::nanoseconds>(
j["Exposure Time"].get<std::string>());
}
// Before v8.0 we had Period instead of Acquisition Period
if (j.contains("Period") && j["Period"].is_string()) {
m_period =
string_to<std::chrono::nanoseconds>(j["Period"].get<std::string>());
}
if (j.contains("Acquisition Period") &&
j["Acquisition Period"].is_string()) {
m_period = string_to<std::chrono::nanoseconds>(
j["Acquisition Period"].get<std::string>());
}
// TODO! Not valid for CTB but not changing api right now!
// Not all detectors write the bitdepth but in case
// its not there it is 16
if(j.contains("Bit Depth") && j["Bit Depth"].is_number()){
m_bitdepth = j["Bit Depth"];
} else {
try {
m_bitdepth = j.at("Dynamic Range");
} catch (const json::out_of_range &e) {
m_bitdepth = 16;
}
m_total_frames_expected = j["Total Frames"];
m_frame_padding = j["Frame Padding"];
m_frame_discard_policy = string_to<FrameDiscardPolicy>(
m_frame_discard_policy = StringTo<FrameDiscardPolicy>(
j["Frame Discard Policy"].get<std::string>());
if(j.contains("Number of rows") && j["Number of rows"].is_number()){
m_number_of_rows = j["Number of rows"];
try {
m_number_of_rows = j.at("Number of rows");
} catch (const json::out_of_range &e) {
// keep the optional empty
}
// ----------------------------------------------------------------
// Special treatment of analog flag because of Moench03
try {
@@ -371,18 +334,19 @@ void RawMasterFile::parse_json(std::istream &is) {
} catch (const json::out_of_range &e) {
// leave the optional empty
}
if (j.contains("Counter Mask")) {
if (j["Counter Mask"].is_number())
m_counter_mask = j["Counter Mask"];
else if (j["Counter Mask"].is_string())
m_counter_mask =
std::stoi(j["Counter Mask"].get<std::string>(), nullptr, 16);
try {
// TODO: what is the best format to handle
m_counter_mask = j.at("Counter Mask");
} catch (const json::out_of_range &e) {
// leave the optional empty
}
// Update detector type for Moench
// TODO! How does this work with old .raw master files?
// Update detector type for Moench
// TODO! How does this work with old .raw master files?
#ifdef AARE_VERBOSE
fmt::print("Detecting Moench03: m_pixels_y: {}, m_analog_samples: {}\n",
m_pixels_y, m_analog_samples.value_or(0));
#endif
if (m_type == DetectorType::Moench && !m_analog_samples &&
m_pixels_y == 400) {
m_type = DetectorType::Moench03;
@@ -391,8 +355,10 @@ void RawMasterFile::parse_json(std::istream &is) {
m_type = DetectorType::Moench03_old;
}
}
void RawMasterFile::parse_raw(std::istream &is) {
for (std::string line; std::getline(is, line);) {
void RawMasterFile::parse_raw(const std::filesystem::path &fpath) {
std::ifstream ifs(fpath);
for (std::string line; std::getline(ifs, line);) {
if (line == "#Frame Header")
break;
auto pos = line.find(':');
@@ -417,12 +383,12 @@ void RawMasterFile::parse_raw(std::istream &is) {
} else if (key == "TimeStamp") {
} else if (key == "Detector Type") {
m_type = string_to<DetectorType>(value);
m_type = StringTo<DetectorType>(value);
if (m_type == DetectorType::Moench) {
m_type = DetectorType::Moench03_old;
}
} else if (key == "Timing Mode") {
m_timing_mode = string_to<TimingMode>(value);
m_timing_mode = StringTo<TimingMode>(value);
} else if (key == "Image Size") {
m_image_size_in_bytes = std::stoi(value);
} else if (key == "Frame Padding") {
@@ -463,10 +429,6 @@ void RawMasterFile::parse_raw(std::istream &is) {
m_pixels_x = std::stoi(value.substr(0, pos));
} else if (key == "Total Frames") {
m_total_frames_expected = std::stoi(value);
} else if (key == "Exptime") {
m_exptime = string_to<std::chrono::nanoseconds>(value);
} else if (key == "Period") {
m_period = string_to<std::chrono::nanoseconds>(value);
} else if (key == "Dynamic Range") {
m_bitdepth = std::stoi(value);
} else if (key == "Quad") {

View File

@@ -4,7 +4,6 @@
#include "test_config.hpp"
#include <catch2/catch_test_macros.hpp>
#include <iostream>
#include <sstream>
using namespace aare;
@@ -165,7 +164,7 @@ TEST_CASE("Parse a master file in .raw format", "[.integration]") {
auto fpath =
test_data_path() /
"raw/moench04/"
"moench/"
"moench04_noise_200V_sto_both_100us_no_light_thresh_900_master_0.raw";
REQUIRE(std::filesystem::exists(fpath));
RawMasterFile f(fpath);
@@ -195,9 +194,7 @@ TEST_CASE("Parse a master file in .raw format", "[.integration]") {
// Total Frames : 100
REQUIRE(f.total_frames_expected() == 100);
// Exptime : 100us
REQUIRE(f.exptime() == std::chrono::microseconds(100));
// Period : 4ms
REQUIRE(f.period() == std::chrono::milliseconds(4));
// Ten Giga : 1
// ADC Mask : 0xffffffff
// Analog Flag : 1
@@ -258,13 +255,13 @@ TEST_CASE("Parse a master file in new .json format",
auto roi = f.roi().value();
REQUIRE(roi.xmin == 0);
REQUIRE(roi.xmax == 2560);
REQUIRE(roi.ymin == 0);
REQUIRE(roi.ymax == 1);
REQUIRE(roi.xmax == 2559);
REQUIRE(roi.ymin == -1);
REQUIRE(roi.ymax == -1);
}
TEST_CASE("Read eiger master file", "[.integration]") {
auto fpath = test_data_path() / "raw/eiger/eiger_500k_32bit_master_0.json";
auto fpath = test_data_path() / "eiger" / "eiger_500k_32bit_master_0.json";
REQUIRE(std::filesystem::exists(fpath));
RawMasterFile f(fpath);
@@ -306,9 +303,7 @@ TEST_CASE("Read eiger master file", "[.integration]") {
// "Dynamic Range": 32,
// "Ten Giga": 0,
// "Exptime": "5s",
REQUIRE(f.exptime() == std::chrono::seconds(5));
// "Period": "1s",
REQUIRE(f.period() == std::chrono::seconds(1));
// "Threshold Energy": -1,
// "Sub Exptime": "2.62144ms",
// "Sub Period": "2.62144ms",
@@ -334,392 +329,3 @@ TEST_CASE("Read eiger master file", "[.integration]") {
// }
// }
}
TEST_CASE("Parse EIGER 7.2 master from string stream") {
std::string master_content = R"({
"Version": 7.2,
"Timestamp": "Tue Mar 26 17:24:34 2024",
"Detector Type": "Eiger",
"Timing Mode": "auto",
"Geometry": {
"x": 2,
"y": 2
},
"Image Size in bytes": 524288,
"Pixels": {
"x": 512,
"y": 256
},
"Max Frames Per File": 10000,
"Frame Discard Policy": "nodiscard",
"Frame Padding": 1,
"Scan Parameters": "[disabled]",
"Total Frames": 3,
"Receiver Roi": {
"xmin": 4294967295,
"xmax": 4294967295,
"ymin": 4294967295,
"ymax": 4294967295
},
"Dynamic Range": 32,
"Ten Giga": 0,
"Exptime": "5s",
"Period": "1s",
"Threshold Energy": -1,
"Sub Exptime": "2.62144ms",
"Sub Period": "2.62144ms",
"Quad": 0,
"Number of rows": 256,
"Rate Corrections": "[0, 0]",
"Frames in File": 3,
"Frame Header Format": {
"Frame Number": "8 bytes",
"SubFrame Number/ExpLength": "4 bytes",
"Packet Number": "4 bytes",
"Bunch ID": "8 bytes",
"Timestamp": "8 bytes",
"Module Id": "2 bytes",
"Row": "2 bytes",
"Column": "2 bytes",
"Reserved": "2 bytes",
"Debug": "4 bytes",
"Round Robin Number": "2 bytes",
"Detector Type": "1 byte",
"Header Version": "1 byte",
"Packets Caught Mask": "64 bytes"
}
})";
std::istringstream iss(master_content);
RawMasterFile f(iss, "test_master_0.json");
REQUIRE(f.version() == "7.2");
REQUIRE(f.detector_type() == DetectorType::Eiger);
REQUIRE(f.timing_mode() == TimingMode::Auto);
REQUIRE(f.geometry().col == 2);
REQUIRE(f.geometry().row == 2);
REQUIRE(f.image_size_in_bytes() == 524288);
REQUIRE(f.pixels_x() == 512);
REQUIRE(f.pixels_y() == 256);
REQUIRE(f.max_frames_per_file() == 10000);
REQUIRE(f.frame_discard_policy() == FrameDiscardPolicy::NoDiscard);
REQUIRE(f.frame_padding() == 1);
REQUIRE(f.total_frames_expected() == 3);
REQUIRE(f.exptime() == std::chrono::seconds(5));
REQUIRE(f.period() == std::chrono::seconds(1));
}
TEST_CASE("Parse JUNGFRAU 7.2 master from string stream") {
std::string master_content = R"({
"Version": 7.2,
"Timestamp": "Tue Feb 20 08:29:19 2024",
"Detector Type": "Jungfrau",
"Timing Mode": "auto",
"Geometry": {
"x": 1,
"y": 2
},
"Image Size in bytes": 524288,
"Pixels": {
"x": 1024,
"y": 256
},
"Max Frames Per File": 3,
"Frame Discard Policy": "nodiscard",
"Frame Padding": 1,
"Scan Parameters": "[disabled]",
"Total Frames": 10,
"Receiver Roi": {
"xmin": 4294967295,
"xmax": 4294967295,
"ymin": 4294967295,
"ymax": 4294967295
},
"Exptime": "10us",
"Period": "1ms",
"Number of UDP Interfaces": 2,
"Number of rows": 512,
"Frames in File": 10,
"Frame Header Format": {
"Frame Number": "8 bytes",
"SubFrame Number/ExpLength": "4 bytes",
"Packet Number": "4 bytes",
"Bunch ID": "8 bytes",
"Timestamp": "8 bytes",
"Module Id": "2 bytes",
"Row": "2 bytes",
"Column": "2 bytes",
"Reserved": "2 bytes",
"Debug": "4 bytes",
"Round Robin Number": "2 bytes",
"Detector Type": "1 byte",
"Header Version": "1 byte",
"Packets Caught Mask": "64 bytes"
}
})";
std::istringstream iss(master_content);
RawMasterFile f(iss, "test_master_0.json");
REQUIRE(f.version() == "7.2");
REQUIRE(f.detector_type() == DetectorType::Jungfrau);
REQUIRE(f.timing_mode() == TimingMode::Auto);
REQUIRE(f.geometry().col == 1);
REQUIRE(f.geometry().row == 2);
REQUIRE(f.n_modules() == 2);
REQUIRE(f.image_size_in_bytes() == 524288);
REQUIRE(f.pixels_x() == 1024);
REQUIRE(f.pixels_y() == 256);
REQUIRE(f.max_frames_per_file() == 3);
REQUIRE(f.frame_discard_policy() == FrameDiscardPolicy::NoDiscard);
REQUIRE(f.frame_padding() == 1);
REQUIRE(f.total_frames_expected() == 10);
REQUIRE(f.exptime() == std::chrono::microseconds(10));
REQUIRE(f.period() == std::chrono::milliseconds(1));
REQUIRE(f.number_of_rows() == 512);
REQUIRE(f.frames_in_file() == 10);
REQUIRE(f.udp_interfaces_per_module() == xy{2, 1});
}
TEST_CASE("Parse a CTB file from stream"){
std::string master_content = R"({
"Version": 8.0,
"Timestamp": "Mon Dec 15 10:57:27 2025",
"Detector Type": "ChipTestBoard",
"Timing Mode": "auto",
"Geometry": {
"x": 1,
"y": 1
},
"Image Size": 18432,
"Pixels": {
"x": 2,
"y": 1
},
"Max Frames Per File": 20000,
"Frame Discard Policy": "nodiscard",
"Frame Padding": 1,
"Scan Parameters": {
"enable": 0,
"dacInd": 0,
"start offset": 0,
"stop offset": 0,
"step size": 0,
"dac settle time ns": 0
},
"Total Frames": 1,
"Exposure Time": "0.25s",
"Acquisition Period": "10ms",
"Ten Giga": 1,
"ADC Mask": 4294967295,
"Analog Flag": 0,
"Analog Samples": 1,
"Digital Flag": 0,
"Digital Samples": 1,
"Dbit Offset": 0,
"Dbit Reorder": 1,
"Dbit Bitset": 0,
"Transceiver Mask": 3,
"Transceiver Flag": 1,
"Transceiver Samples": 1152,
"Frames in File": 40,
"Additional JSON Header": {}
})";
std::istringstream iss(master_content);
RawMasterFile f(iss, "test_master_0.json");
REQUIRE(f.version() == "8.0");
REQUIRE(f.detector_type() == DetectorType::ChipTestBoard);
REQUIRE(f.timing_mode() == TimingMode::Auto);
REQUIRE(f.geometry().col == 1);
REQUIRE(f.geometry().row == 1);
REQUIRE(f.image_size_in_bytes() == 18432);
REQUIRE(f.pixels_x() == 2);
REQUIRE(f.pixels_y() == 1);
REQUIRE(f.max_frames_per_file() == 20000);
// CTB does not have bitdepth in master file, but for the moment we write 16
// TODO! refactor using std::optional
// REQUIRE(f.bitdepth() == std::nullopt);
REQUIRE(f.n_modules() == 1);
REQUIRE(f.quad() == 0);
REQUIRE(f.frame_discard_policy() == FrameDiscardPolicy::NoDiscard);
REQUIRE(f.frame_padding() == 1);
REQUIRE(f.total_frames_expected() == 1); //This is Total Frames in the master file
REQUIRE(f.exptime() == std::chrono::milliseconds(250));
REQUIRE(f.period() == std::chrono::milliseconds(10));
REQUIRE(f.analog_samples() == std::nullopt); //Analog Flag is 0
REQUIRE(f.digital_samples() == std::nullopt); //Digital Flag is 0
REQUIRE(f.transceiver_samples() == 1152);
REQUIRE(f.frames_in_file() == 40);
}
TEST_CASE("Parse v8.0 MYTHEN3 from stream"){
std::string master_content = R"({
"Version": 8.0,
"Timestamp": "Wed Oct 1 14:37:26 2025",
"Detector Type": "Mythen3",
"Timing Mode": "auto",
"Geometry": {
"x": 2,
"y": 1
},
"Image Size": 5120,
"Pixels": {
"x": 1280,
"y": 1
},
"Max Frames Per File": 10000,
"Frame Discard Policy": "nodiscard",
"Frame Padding": 1,
"Scan Parameters": {
"enable": 0,
"dacInd": 0,
"start offset": 0,
"stop offset": 0,
"step size": 0,
"dac settle time ns": 0
},
"Total Frames": 1,
"Receiver Rois": [
{
"xmin": 0,
"xmax": 2559,
"ymin": -1,
"ymax": -1
}
],
"Dynamic Range": 32,
"Ten Giga": 1,
"Acquisition Period": "0ns",
"Counter Mask": 4,
"Exposure Times": [
"5s",
"5s",
"5s"
],
"Gate Delays": [
"0ns",
"0ns",
"0ns"
],
"Gates": 1,
"Threshold Energies": [
-1,
-1,
-1
],
"Readout Speed": "half_speed",
"Frames in File": 1,
"Additional JSON Header": {}
})";
std::istringstream iss(master_content);
RawMasterFile f(iss, "test_master_0.json");
REQUIRE(f.version() == "8.0");
REQUIRE(f.detector_type() == DetectorType::Mythen3);
REQUIRE(f.timing_mode() == TimingMode::Auto);
REQUIRE(f.geometry().col == 2);
REQUIRE(f.geometry().row == 1);
REQUIRE(f.image_size_in_bytes() == 5120);
REQUIRE(f.pixels_x() == 1280);
REQUIRE(f.pixels_y() == 1);
REQUIRE(f.max_frames_per_file() == 10000);
REQUIRE(f.n_modules() == 2);
REQUIRE(f.quad() == 0);
REQUIRE(f.frame_discard_policy() == FrameDiscardPolicy::NoDiscard);
REQUIRE(f.frame_padding() == 1);
REQUIRE(f.total_frames_expected() == 1); //This is Total Frames in the master file
REQUIRE(f.counter_mask() == 4);
// Mythen3 has three exposure times, but for the moment we don't handle them
REQUIRE(f.exptime() == std::nullopt);
// Period is ok though
REQUIRE(f.period() == std::chrono::nanoseconds(0));
}
TEST_CASE("Parse a v7.1 Mythen3 from stream"){
std::string master_content = R"({
"Version": 7.1,
"Timestamp": "Wed Sep 21 13:48:10 2022",
"Detector Type": "Mythen3",
"Timing Mode": "auto",
"Geometry": {
"x": 1,
"y": 1
},
"Image Size in bytes": 15360,
"Pixels": {
"x": 3840,
"y": 1
},
"Max Frames Per File": 10000,
"Frame Discard Policy": "nodiscard",
"Frame Padding": 1,
"Scan Parameters": "[disabled]",
"Total Frames": 1,
"Receiver Roi": {
"xmin": 4294967295,
"xmax": 4294967295,
"ymin": 4294967295,
"ymax": 4294967295
},
"Dynamic Range": 32,
"Ten Giga": 1,
"Period": "2ms",
"Counter Mask": "0x7",
"Exptime1": "0.1s",
"Exptime2": "0.1s",
"Exptime3": "0.1s",
"GateDelay1": "0ns",
"GateDelay2": "0ns",
"GateDelay3": "0ns",
"Gates": 1,
"Threshold Energies": "[-1, -1, -1]",
"Frames in File": 1,
"Frame Header Format": {
"Frame Number": "8 bytes",
"SubFrame Number/ExpLength": "4 bytes",
"Packet Number": "4 bytes",
"Bunch ID": "8 bytes",
"Timestamp": "8 bytes",
"Module Id": "2 bytes",
"Row": "2 bytes",
"Column": "2 bytes",
"Reserved": "2 bytes",
"Debug": "4 bytes",
"Round Robin Number": "2 bytes",
"Detector Type": "1 byte",
"Header Version": "1 byte",
"Packets Caught Mask": "64 bytes"
}
})";
std::istringstream iss(master_content);
RawMasterFile f(iss, "test_master_0.json");
REQUIRE(f.version() == "7.1");
REQUIRE(f.detector_type() == DetectorType::Mythen3);
REQUIRE(f.timing_mode() == TimingMode::Auto);
REQUIRE(f.geometry().col == 1);
REQUIRE(f.geometry().row == 1);
REQUIRE(f.image_size_in_bytes() == 15360);
REQUIRE(f.pixels_x() == 3840);
REQUIRE(f.pixels_y() == 1);
REQUIRE(f.max_frames_per_file() == 10000);
REQUIRE(f.n_modules() == 1);
REQUIRE(f.quad() == 0);
REQUIRE(f.frame_discard_policy() == FrameDiscardPolicy::NoDiscard);
REQUIRE(f.frame_padding() == 1);
REQUIRE(f.total_frames_expected() == 1); //This is Total Frames in the master file
REQUIRE(f.counter_mask() == 0x7);
// Mythen3 has three exposure times, but for the moment we don't handle them
REQUIRE(f.exptime() == std::nullopt);
// Period is ok though
REQUIRE(f.period() == std::chrono::milliseconds(2));
}

View File

@@ -1,6 +1,5 @@
// SPDX-License-Identifier: MPL-2.0
#include "aare/decode.hpp"
#include <fmt/format.h>
#include <cmath>
namespace aare {
@@ -106,49 +105,4 @@ void apply_custom_weights(NDView<uint16_t, 1> input, NDView<double, 1> output,
}
}
uint32_t mask32to24bits(uint32_t input, BitOffset offset){
constexpr uint32_t mask24bits{0xFFFFFF};
return (input >> offset.value()) & mask24bits;
}
void expand24to32bit(NDView<uint8_t,1> input, NDView<uint32_t,1> output, BitOffset bit_offset){
ssize_t bytes_per_channel = 3; //24bit
ssize_t min_input_size = output.size()*bytes_per_channel;
//if we have an offset we need one more byte in the input data
if (bit_offset.value())
min_input_size += 1;
if (input.size() < min_input_size)
throw std::runtime_error(fmt::format(
"{} Mismatch between input and output size. Output "
"size of {} with bit offset {} requires an input of at least {} "
"bytes. Called with input size: {} output size: {}",
LOCATION, output.size(), bit_offset.value(), min_input_size, input.size(), output.size()));
auto* in = input.data();
if(bit_offset.value()){
//If there is a bit_offset we copy 4 bytes and then
//mask out the correct ones.
for (auto& v : output){
uint32_t val{};
std::memcpy(&val, in, sizeof(val));
v = mask32to24bits(val, bit_offset);
in += bytes_per_channel;
}
}else{
//If there is no offset we can directly copy the bits
//without masking
for (auto& v : output){
uint32_t val{};
std::memcpy(&val, in, 3);
v = val;
in += bytes_per_channel;
}
}
}
} // namespace aare

View File

@@ -7,8 +7,6 @@
using Catch::Matchers::WithinAbs;
#include <vector>
using aare::BitOffset;
TEST_CASE("test_adc_sar_05_decode64to16") {
uint64_t input = 0;
uint16_t output = aare::adc_sar_05_decode64to16(input);
@@ -74,93 +72,3 @@ TEST_CASE("test_apply_custom_weights") {
output = aare::apply_custom_weights(input, weights);
CHECK_THAT(output, WithinAbs(6.34, 0.001));
}
TEST_CASE("Mask 32 bit unsigned integer to 24 bit"){
//any number less than 2**24 (16777216) should be the same
CHECK(aare::mask32to24bits(0)==0);
CHECK(aare::mask32to24bits(19)==19);
CHECK(aare::mask32to24bits(29875)==29875);
CHECK(aare::mask32to24bits(1092177)==1092177);
CHECK(aare::mask32to24bits(0xFFFF)==0xFFFF);
CHECK(aare::mask32to24bits(0xFFFFFFFF)==0xFFFFFF);
// Offset specifies that the should ignore 0-7 bits
// at the start
CHECK(aare::mask32to24bits(0xFFFF, BitOffset(4))==0xFFF);
CHECK(aare::mask32to24bits(0xFF0000d9)==0xd9);
CHECK(aare::mask32to24bits(0xFF000d9F, BitOffset(4))==0xF000d9);
CHECK(aare::mask32to24bits(16777217)==1);
CHECK(aare::mask32to24bits(15,BitOffset(7))==0);
//Highest bit set to 1 should just be excluded
//lowest 4 bits set to 1
CHECK(aare::mask32to24bits(0x8000000f,BitOffset(7))==0);
}
TEST_CASE("Expand container with 24 bit data to 32"){
{
uint8_t buffer[] = {
0x00, 0x00, 0x00,
0x00, 0x00, 0x00,
0x00, 0x00, 0x00,
};
aare::NDView<uint8_t, 1> input(&buffer[0], {9});
aare::NDArray<uint32_t, 1> out({3});
aare::expand24to32bit(input, out.view());
CHECK(out(0) == 0);
CHECK(out(1) == 0);
CHECK(out(2) == 0);
}
{
uint8_t buffer[] = {
0x0F, 0x00, 0x00,
0xFF, 0x00, 0x00,
0xFF, 0xFF, 0xFF,
};
aare::NDView<uint8_t, 1> input(&buffer[0], {9});
aare::NDArray<uint32_t, 1> out({3});
aare::expand24to32bit(input, out.view());
CHECK(out(0) == 0xF);
CHECK(out(1) == 0xFF);
CHECK(out(2) == 0xFFFFFF);
}
{
uint8_t buffer[] = {
0x00, 0x00, 0xFF,
0xFF, 0xFF, 0x00,
0x00, 0xFF, 0x00,
};
aare::NDView<uint8_t, 1> input(&buffer[0], {9});
aare::NDArray<uint32_t, 1> out({3});
aare::expand24to32bit(input, out.view());
CHECK(out(0) == 0xFF0000);
CHECK(out(1) == 0xFFFF);
CHECK(out(2) == 0xFF00);
REQUIRE_THROWS(aare::expand24to32bit(input, out.view(), BitOffset(4)));
}
{
//For use with offset we need an extra byte
uint8_t buffer[] = {
0x00, 0x00, 0xFF,
0xFF, 0xFF, 0x00,
0x00, 0xFF, 0x00, 0x00
};
aare::NDView<uint8_t, 1> input(&buffer[0], {10});
aare::NDArray<uint32_t, 1> out({3}); //still output.size == 3
aare::expand24to32bit(input, out.view(), BitOffset(4));
CHECK(out(0) == 0xFFF000);
CHECK(out(1) == 0xFFF);
CHECK(out(2) == 0xFF0);
}
}

View File

@@ -11,63 +11,291 @@ void assert_failed(const std::string &msg) {
exit(1);
}
// /**
// * @brief Convert a DetectorType to a string
// * @param type DetectorType
// * @return string representation of the DetectorType
// */
// template <> std::string ToString(DetectorType arg) {
// switch (arg) {
// case DetectorType::Generic:
// return "Generic";
// case DetectorType::Eiger:
// return "Eiger";
// case DetectorType::Gotthard:
// return "Gotthard";
// case DetectorType::Jungfrau:
// return "Jungfrau";
// case DetectorType::ChipTestBoard:
// return "ChipTestBoard";
// case DetectorType::Moench:
// return "Moench";
// case DetectorType::Mythen3:
// return "Mythen3";
// case DetectorType::Gotthard2:
// return "Gotthard2";
// case DetectorType::Xilinx_ChipTestBoard:
// return "Xilinx_ChipTestBoard";
/**
* @brief Convert a DetectorType to a string
* @param type DetectorType
* @return string representation of the DetectorType
*/
template <> std::string ToString(DetectorType arg) {
switch (arg) {
case DetectorType::Generic:
return "Generic";
case DetectorType::Eiger:
return "Eiger";
case DetectorType::Gotthard:
return "Gotthard";
case DetectorType::Jungfrau:
return "Jungfrau";
case DetectorType::ChipTestBoard:
return "ChipTestBoard";
case DetectorType::Moench:
return "Moench";
case DetectorType::Mythen3:
return "Mythen3";
case DetectorType::Gotthard2:
return "Gotthard2";
case DetectorType::Xilinx_ChipTestBoard:
return "Xilinx_ChipTestBoard";
// // Custom ones
// case DetectorType::Moench03:
// return "Moench03";
// case DetectorType::Moench03_old:
// return "Moench03_old";
// case DetectorType::Unknown:
// return "Unknown";
// // no default case to trigger compiler warning if not all
// // enum values are handled
// }
// throw std::runtime_error("Could not decode detector to string");
// }
BitOffset::BitOffset(uint32_t offset){
if (offset>7)
throw std::runtime_error(fmt::format("{} BitOffset needs to be <8: Called with {}", LOCATION, offset));
m_offset = static_cast<uint8_t>(offset);
// Custom ones
case DetectorType::Moench03:
return "Moench03";
case DetectorType::Moench03_old:
return "Moench03_old";
case DetectorType::Unknown:
return "Unknown";
// no default case to trigger compiler warning if not all
// enum values are handled
}
throw std::runtime_error("Could not decode detector to string");
}
bool BitOffset::operator==(const BitOffset& other) const {
return m_offset == other.m_offset;
}
/**
* @brief Convert a string to a DetectorType
* @param name string representation of the DetectorType
* @return DetectorType
* @throw runtime_error if the string does not match any DetectorType
*/
template <> DetectorType StringTo(const std::string &arg) {
if (arg == "Generic")
return DetectorType::Generic;
if (arg == "Eiger")
return DetectorType::Eiger;
if (arg == "Gotthard")
return DetectorType::Gotthard;
if (arg == "Jungfrau")
return DetectorType::Jungfrau;
if (arg == "ChipTestBoard")
return DetectorType::ChipTestBoard;
if (arg == "Moench")
return DetectorType::Moench;
if (arg == "Mythen3")
return DetectorType::Mythen3;
if (arg == "Gotthard2")
return DetectorType::Gotthard2;
if (arg == "Xilinx_ChipTestBoard")
return DetectorType::Xilinx_ChipTestBoard;
bool BitOffset::operator<(const BitOffset& other) const {
return m_offset < other.m_offset;
}
// Custom ones
if (arg == "Moench03")
return DetectorType::Moench03;
if (arg == "Moench03_old")
return DetectorType::Moench03_old;
if (arg == "Unknown")
return DetectorType::Unknown;
throw std::runtime_error("Could not decode detector from: \"" + arg + "\"");
}
/**
* @brief Convert a string to a TimingMode
* @param mode string representation of the TimingMode
* @return TimingMode
* @throw runtime_error if the string does not match any TimingMode
*/
template <> TimingMode StringTo(const std::string &arg) {
if (arg == "auto")
return TimingMode::Auto;
if (arg == "trigger")
return TimingMode::Trigger;
throw std::runtime_error("Could not decode timing mode from: \"" + arg +
"\"");
}
template <> FrameDiscardPolicy StringTo(const std::string &arg) {
if (arg == "nodiscard")
return FrameDiscardPolicy::NoDiscard;
if (arg == "discard")
return FrameDiscardPolicy::Discard;
if (arg == "discardpartial")
return FrameDiscardPolicy::DiscardPartial;
throw std::runtime_error("Could not decode frame discard policy from: \"" +
arg + "\"");
}
// template <> TimingMode StringTo<TimingMode>(std::string mode);
template <> DACIndex StringTo(const std::string &arg) {
if (arg == "dac 0")
return DACIndex::DAC_0;
else if (arg == "dac 1")
return DACIndex::DAC_1;
else if (arg == "dac 2")
return DACIndex::DAC_2;
else if (arg == "dac 3")
return DACIndex::DAC_3;
else if (arg == "dac 4")
return DACIndex::DAC_4;
else if (arg == "dac 5")
return DACIndex::DAC_5;
else if (arg == "dac 6")
return DACIndex::DAC_6;
else if (arg == "dac 7")
return DACIndex::DAC_7;
else if (arg == "dac 8")
return DACIndex::DAC_8;
else if (arg == "dac 9")
return DACIndex::DAC_9;
else if (arg == "dac 10")
return DACIndex::DAC_10;
else if (arg == "dac 11")
return DACIndex::DAC_11;
else if (arg == "dac 12")
return DACIndex::DAC_12;
else if (arg == "dac 13")
return DACIndex::DAC_13;
else if (arg == "dac 14")
return DACIndex::DAC_14;
else if (arg == "dac 15")
return DACIndex::DAC_15;
else if (arg == "dac 16")
return DACIndex::DAC_16;
else if (arg == "dac 17")
return DACIndex::DAC_17;
else if (arg == "vsvp")
return DACIndex::VSVP;
else if (arg == "vtrim")
return DACIndex::VTRIM;
else if (arg == "vrpreamp")
return DACIndex::VRPREAMP;
else if (arg == "vrshaper")
return DACIndex::VRSHAPER;
else if (arg == "vsvn")
return DACIndex::VSVN;
else if (arg == "vtgstv")
return DACIndex::VTGSTV;
else if (arg == "vcmp_ll")
return DACIndex::VCMP_LL;
else if (arg == "vcmp_lr")
return DACIndex::VCMP_LR;
else if (arg == "vcal")
return DACIndex::VCAL;
else if (arg == "vcmp_rl")
return DACIndex::VCMP_RL;
else if (arg == "rxb_rb")
return DACIndex::RXB_RB;
else if (arg == "rxb_lb")
return DACIndex::RXB_LB;
else if (arg == "vcmp_rr")
return DACIndex::VCMP_RR;
else if (arg == "vcp")
return DACIndex::VCP;
else if (arg == "vcn")
return DACIndex::VCN;
else if (arg == "vishaper")
return DACIndex::VISHAPER;
else if (arg == "vthreshold")
return DACIndex::VTHRESHOLD;
else if (arg == "vref_ds")
return DACIndex::VREF_DS;
else if (arg == "vout_cm")
return DACIndex::VOUT_CM;
else if (arg == "vin_cm")
return DACIndex::VIN_CM;
else if (arg == "vref_comp")
return DACIndex::VREF_COMP;
else if (arg == "vb_comp")
return DACIndex::VB_COMP;
else if (arg == "vdd_prot")
return DACIndex::VDD_PROT;
else if (arg == "vin_com")
return DACIndex::VIN_COM;
else if (arg == "vref_prech")
return DACIndex::VREF_PRECH;
else if (arg == "vb_pixbuf")
return DACIndex::VB_PIXBUF;
else if (arg == "vb_ds")
return DACIndex::VB_DS;
else if (arg == "vref_h_adc")
return DACIndex::VREF_H_ADC;
else if (arg == "vb_comp_fe")
return DACIndex::VB_COMP_FE;
else if (arg == "vb_comp_adc")
return DACIndex::VB_COMP_ADC;
else if (arg == "vcom_cds")
return DACIndex::VCOM_CDS;
else if (arg == "vref_rstore")
return DACIndex::VREF_RSTORE;
else if (arg == "vb_opa_1st")
return DACIndex::VB_OPA_1ST;
else if (arg == "vref_comp_fe")
return DACIndex::VREF_COMP_FE;
else if (arg == "vcom_adc1")
return DACIndex::VCOM_ADC1;
else if (arg == "vref_l_adc")
return DACIndex::VREF_L_ADC;
else if (arg == "vref_cds")
return DACIndex::VREF_CDS;
else if (arg == "vb_cs")
return DACIndex::VB_CS;
else if (arg == "vb_opa_fd")
return DACIndex::VB_OPA_FD;
else if (arg == "vcom_adc2")
return DACIndex::VCOM_ADC2;
else if (arg == "vcassh")
return DACIndex::VCASSH;
else if (arg == "vth2")
return DACIndex::VTH2;
else if (arg == "vrshaper_n")
return DACIndex::VRSHAPER_N;
else if (arg == "vipre_out")
return DACIndex::VIPRE_OUT;
else if (arg == "vth3")
return DACIndex::VTH3;
else if (arg == "vth1")
return DACIndex::VTH1;
else if (arg == "vicin")
return DACIndex::VICIN;
else if (arg == "vcas")
return DACIndex::VCAS;
else if (arg == "vcal_n")
return DACIndex::VCAL_N;
else if (arg == "vipre")
return DACIndex::VIPRE;
else if (arg == "vcal_p")
return DACIndex::VCAL_P;
else if (arg == "vdcsh")
return DACIndex::VDCSH;
else if (arg == "vbp_colbuf")
return DACIndex::VBP_COLBUF;
else if (arg == "vb_sda")
return DACIndex::VB_SDA;
else if (arg == "vcasc_sfp")
return DACIndex::VCASC_SFP;
else if (arg == "vipre_cds")
return DACIndex::VIPRE_CDS;
else if (arg == "ibias_sfp")
return DACIndex::IBIAS_SFP;
else if (arg == "trimbits")
return DACIndex::TRIMBIT_SCAN;
else if (arg == "highvoltage")
return DACIndex::HIGH_VOLTAGE;
else if (arg == "iodelay")
return DACIndex::IO_DELAY;
else if (arg == "temp_adc")
return DACIndex::TEMPERATURE_ADC;
else if (arg == "temp_fpga")
return DACIndex::TEMPERATURE_FPGA;
else if (arg == "temp_fpgaext")
return DACIndex::TEMPERATURE_FPGAEXT;
else if (arg == "temp_10ge")
return DACIndex::TEMPERATURE_10GE;
else if (arg == "temp_dcdc")
return DACIndex::TEMPERATURE_DCDC;
else if (arg == "temp_sodl")
return DACIndex::TEMPERATURE_SODL;
else if (arg == "temp_sodr")
return DACIndex::TEMPERATURE_SODR;
else if (arg == "temp_fpgafl")
return DACIndex::TEMPERATURE_FPGA2;
else if (arg == "temp_fpgafr")
return DACIndex::TEMPERATURE_FPGA3;
else if (arg == "temp_slowadc")
return DACIndex::SLOW_ADC_TEMP;
else
throw std::invalid_argument("Could not decode DACIndex from: \"" + arg +
"\"");
}
} // namespace aare

View File

@@ -4,7 +4,52 @@
#include <catch2/catch_test_macros.hpp>
#include <string>
using aare::StringTo;
using aare::ToString;
TEST_CASE("Enum to string conversion") {
// TODO! By the way I don't think the enum string conversions should be in
// the defs.hpp file but let's use this to show a test
REQUIRE(ToString(aare::DetectorType::Generic) == "Generic");
REQUIRE(ToString(aare::DetectorType::Eiger) == "Eiger");
REQUIRE(ToString(aare::DetectorType::Gotthard) == "Gotthard");
REQUIRE(ToString(aare::DetectorType::Jungfrau) == "Jungfrau");
REQUIRE(ToString(aare::DetectorType::ChipTestBoard) == "ChipTestBoard");
REQUIRE(ToString(aare::DetectorType::Moench) == "Moench");
REQUIRE(ToString(aare::DetectorType::Mythen3) == "Mythen3");
REQUIRE(ToString(aare::DetectorType::Gotthard2) == "Gotthard2");
REQUIRE(ToString(aare::DetectorType::Xilinx_ChipTestBoard) ==
"Xilinx_ChipTestBoard");
REQUIRE(ToString(aare::DetectorType::Moench03) == "Moench03");
REQUIRE(ToString(aare::DetectorType::Moench03_old) == "Moench03_old");
REQUIRE(ToString(aare::DetectorType::Unknown) == "Unknown");
}
TEST_CASE("String to enum") {
REQUIRE(StringTo<aare::DetectorType>("Generic") ==
aare::DetectorType::Generic);
REQUIRE(StringTo<aare::DetectorType>("Eiger") == aare::DetectorType::Eiger);
REQUIRE(StringTo<aare::DetectorType>("Gotthard") ==
aare::DetectorType::Gotthard);
REQUIRE(StringTo<aare::DetectorType>("Jungfrau") ==
aare::DetectorType::Jungfrau);
REQUIRE(StringTo<aare::DetectorType>("ChipTestBoard") ==
aare::DetectorType::ChipTestBoard);
REQUIRE(StringTo<aare::DetectorType>("Moench") ==
aare::DetectorType::Moench);
REQUIRE(StringTo<aare::DetectorType>("Mythen3") ==
aare::DetectorType::Mythen3);
REQUIRE(StringTo<aare::DetectorType>("Gotthard2") ==
aare::DetectorType::Gotthard2);
REQUIRE(StringTo<aare::DetectorType>("Xilinx_ChipTestBoard") ==
aare::DetectorType::Xilinx_ChipTestBoard);
REQUIRE(StringTo<aare::DetectorType>("Moench03") ==
aare::DetectorType::Moench03);
REQUIRE(StringTo<aare::DetectorType>("Moench03_old") ==
aare::DetectorType::Moench03_old);
REQUIRE(StringTo<aare::DetectorType>("Unknown") ==
aare::DetectorType::Unknown);
}
TEST_CASE("Enum values") {
// Since some of the enums are written to file we need to make sure
@@ -38,23 +83,21 @@ TEST_CASE("DynamicCluster creation") {
REQUIRE(c2.data() != nullptr);
}
TEST_CASE("Basic ops on BitOffset"){
REQUIRE_THROWS(aare::BitOffset(10));
// TEST_CASE("cluster set and get data") {
aare::BitOffset offset(5);
REQUIRE(offset.value()==5);
// aare::DynamicCluster c2(33, 44, aare::Dtype(typeid(double)));
// REQUIRE(c2.cluster_sizeX == 33);
// REQUIRE(c2.cluster_sizeY == 44);
// REQUIRE(c2.dt == aare::Dtype::DOUBLE);
// double v = 3.14;
// c2.set<double>(0, v);
// double v2 = c2.get<double>(0);
// REQUIRE(aare::compare_floats<double>(v, v2));
aare::BitOffset offset2;
REQUIRE(offset2.value()==0);
aare::BitOffset offset3(offset);
REQUIRE(offset3.value()==5);
REQUIRE(offset==offset3);
//Now assign offset to offset2 which should get the value 5
offset2 = offset;
REQUIRE(offset2.value()==5);
REQUIRE(offset2==offset);
}
// c2.set<double>(33 * 44 - 1, 123.11);
// double v3 = c2.get<double>(33 * 44 - 1);
// REQUIRE(aare::compare_floats<double>(123.11, v3));
// REQUIRE_THROWS_AS(c2.set(0, 1), std::invalid_argument); // set int to
// double
// }

View File

@@ -1,276 +0,0 @@
#include "to_string.hpp"
namespace aare {
template <> DetectorType string_to(const std::string &arg) {
if (arg == "Generic")
return DetectorType::Generic;
if (arg == "Eiger")
return DetectorType::Eiger;
if (arg == "Gotthard")
return DetectorType::Gotthard;
if (arg == "Jungfrau")
return DetectorType::Jungfrau;
if (arg == "ChipTestBoard")
return DetectorType::ChipTestBoard;
if (arg == "Moench")
return DetectorType::Moench;
if (arg == "Mythen3")
return DetectorType::Mythen3;
if (arg == "Gotthard2")
return DetectorType::Gotthard2;
if (arg == "Xilinx_ChipTestBoard")
return DetectorType::Xilinx_ChipTestBoard;
// Custom ones
if (arg == "Moench03")
return DetectorType::Moench03;
if (arg == "Moench03_old")
return DetectorType::Moench03_old;
if (arg == "Unknown")
return DetectorType::Unknown;
throw std::runtime_error("Could not decode detector from: \"" + arg + "\"");
}
template <> TimingMode string_to(const std::string &arg) {
if (arg == "auto")
return TimingMode::Auto;
if (arg == "trigger")
return TimingMode::Trigger;
throw std::runtime_error("Could not decode timing mode from: \"" + arg +
"\"");
}
template <> FrameDiscardPolicy string_to(const std::string &arg) {
if (arg == "nodiscard")
return FrameDiscardPolicy::NoDiscard;
if (arg == "discard")
return FrameDiscardPolicy::Discard;
if (arg == "discardpartial")
return FrameDiscardPolicy::DiscardPartial;
throw std::runtime_error("Could not decode frame discard policy from: \"" +
arg + "\"");
}
template <> DACIndex string_to(const std::string &arg) {
if (arg == "dac 0")
return DACIndex::DAC_0;
else if (arg == "dac 1")
return DACIndex::DAC_1;
else if (arg == "dac 2")
return DACIndex::DAC_2;
else if (arg == "dac 3")
return DACIndex::DAC_3;
else if (arg == "dac 4")
return DACIndex::DAC_4;
else if (arg == "dac 5")
return DACIndex::DAC_5;
else if (arg == "dac 6")
return DACIndex::DAC_6;
else if (arg == "dac 7")
return DACIndex::DAC_7;
else if (arg == "dac 8")
return DACIndex::DAC_8;
else if (arg == "dac 9")
return DACIndex::DAC_9;
else if (arg == "dac 10")
return DACIndex::DAC_10;
else if (arg == "dac 11")
return DACIndex::DAC_11;
else if (arg == "dac 12")
return DACIndex::DAC_12;
else if (arg == "dac 13")
return DACIndex::DAC_13;
else if (arg == "dac 14")
return DACIndex::DAC_14;
else if (arg == "dac 15")
return DACIndex::DAC_15;
else if (arg == "dac 16")
return DACIndex::DAC_16;
else if (arg == "dac 17")
return DACIndex::DAC_17;
else if (arg == "vsvp")
return DACIndex::VSVP;
else if (arg == "vtrim")
return DACIndex::VTRIM;
else if (arg == "vrpreamp")
return DACIndex::VRPREAMP;
else if (arg == "vrshaper")
return DACIndex::VRSHAPER;
else if (arg == "vsvn")
return DACIndex::VSVN;
else if (arg == "vtgstv")
return DACIndex::VTGSTV;
else if (arg == "vcmp_ll")
return DACIndex::VCMP_LL;
else if (arg == "vcmp_lr")
return DACIndex::VCMP_LR;
else if (arg == "vcal")
return DACIndex::VCAL;
else if (arg == "vcmp_rl")
return DACIndex::VCMP_RL;
else if (arg == "rxb_rb")
return DACIndex::RXB_RB;
else if (arg == "rxb_lb")
return DACIndex::RXB_LB;
else if (arg == "vcmp_rr")
return DACIndex::VCMP_RR;
else if (arg == "vcp")
return DACIndex::VCP;
else if (arg == "vcn")
return DACIndex::VCN;
else if (arg == "vishaper")
return DACIndex::VISHAPER;
else if (arg == "vthreshold")
return DACIndex::VTHRESHOLD;
else if (arg == "vref_ds")
return DACIndex::VREF_DS;
else if (arg == "vout_cm")
return DACIndex::VOUT_CM;
else if (arg == "vin_cm")
return DACIndex::VIN_CM;
else if (arg == "vref_comp")
return DACIndex::VREF_COMP;
else if (arg == "vb_comp")
return DACIndex::VB_COMP;
else if (arg == "vdd_prot")
return DACIndex::VDD_PROT;
else if (arg == "vin_com")
return DACIndex::VIN_COM;
else if (arg == "vref_prech")
return DACIndex::VREF_PRECH;
else if (arg == "vb_pixbuf")
return DACIndex::VB_PIXBUF;
else if (arg == "vb_ds")
return DACIndex::VB_DS;
else if (arg == "vref_h_adc")
return DACIndex::VREF_H_ADC;
else if (arg == "vb_comp_fe")
return DACIndex::VB_COMP_FE;
else if (arg == "vb_comp_adc")
return DACIndex::VB_COMP_ADC;
else if (arg == "vcom_cds")
return DACIndex::VCOM_CDS;
else if (arg == "vref_rstore")
return DACIndex::VREF_RSTORE;
else if (arg == "vb_opa_1st")
return DACIndex::VB_OPA_1ST;
else if (arg == "vref_comp_fe")
return DACIndex::VREF_COMP_FE;
else if (arg == "vcom_adc1")
return DACIndex::VCOM_ADC1;
else if (arg == "vref_l_adc")
return DACIndex::VREF_L_ADC;
else if (arg == "vref_cds")
return DACIndex::VREF_CDS;
else if (arg == "vb_cs")
return DACIndex::VB_CS;
else if (arg == "vb_opa_fd")
return DACIndex::VB_OPA_FD;
else if (arg == "vcom_adc2")
return DACIndex::VCOM_ADC2;
else if (arg == "vcassh")
return DACIndex::VCASSH;
else if (arg == "vth2")
return DACIndex::VTH2;
else if (arg == "vrshaper_n")
return DACIndex::VRSHAPER_N;
else if (arg == "vipre_out")
return DACIndex::VIPRE_OUT;
else if (arg == "vth3")
return DACIndex::VTH3;
else if (arg == "vth1")
return DACIndex::VTH1;
else if (arg == "vicin")
return DACIndex::VICIN;
else if (arg == "vcas")
return DACIndex::VCAS;
else if (arg == "vcal_n")
return DACIndex::VCAL_N;
else if (arg == "vipre")
return DACIndex::VIPRE;
else if (arg == "vcal_p")
return DACIndex::VCAL_P;
else if (arg == "vdcsh")
return DACIndex::VDCSH;
else if (arg == "vbp_colbuf")
return DACIndex::VBP_COLBUF;
else if (arg == "vb_sda")
return DACIndex::VB_SDA;
else if (arg == "vcasc_sfp")
return DACIndex::VCASC_SFP;
else if (arg == "vipre_cds")
return DACIndex::VIPRE_CDS;
else if (arg == "ibias_sfp")
return DACIndex::IBIAS_SFP;
else if (arg == "trimbits")
return DACIndex::TRIMBIT_SCAN;
else if (arg == "highvoltage")
return DACIndex::HIGH_VOLTAGE;
else if (arg == "iodelay")
return DACIndex::IO_DELAY;
else if (arg == "temp_adc")
return DACIndex::TEMPERATURE_ADC;
else if (arg == "temp_fpga")
return DACIndex::TEMPERATURE_FPGA;
else if (arg == "temp_fpgaext")
return DACIndex::TEMPERATURE_FPGAEXT;
else if (arg == "temp_10ge")
return DACIndex::TEMPERATURE_10GE;
else if (arg == "temp_dcdc")
return DACIndex::TEMPERATURE_DCDC;
else if (arg == "temp_sodl")
return DACIndex::TEMPERATURE_SODL;
else if (arg == "temp_sodr")
return DACIndex::TEMPERATURE_SODR;
else if (arg == "temp_fpgafl")
return DACIndex::TEMPERATURE_FPGA2;
else if (arg == "temp_fpgafr")
return DACIndex::TEMPERATURE_FPGA3;
else if (arg == "temp_slowadc")
return DACIndex::SLOW_ADC_TEMP;
else
throw std::invalid_argument("Could not decode DACIndex from: \"" + arg +
"\"");
}
std::string remove_unit(std::string &str) {
auto it = str.begin();
while (it != str.end()) {
if (std::isalpha(*it)) {
// Check if this is scientific notation (e or E followed by optional sign and digits)
if (((*it == 'e' || *it == 'E') && (it + 1) != str.end())) {
auto next = it + 1;
// Skip optional sign
if (*next == '+' || *next == '-') {
++next;
}
// Check if followed by at least one digit
if (next != str.end() && std::isdigit(*next)) {
// This is scientific notation, continue scanning
it = next;
while (it != str.end() && std::isdigit(*it)) {
++it;
}
continue;
}
}
// Not scientific notation, this is the start of the unit
break;
}
++it;
}
auto pos = it - str.begin();
auto unit = str.substr(pos);
str.erase(it, end(str));
// Strip trailing whitespace
while (!str.empty() && std::isspace(str.back())) {
str.pop_back();
}
return unit;
}
} // namespace aare

View File

@@ -1,90 +0,0 @@
// SPDX-License-Identifier: MPL-2.0
/*
*The file to_string.hpp contains conversion to and from string for various aare
* types. It is intentionally not a part of the public API.
* Only include the conversion that are actually needed!!!
*/
#include "aare/defs.hpp" //enums
#include <chrono> //time conversions
namespace aare {
std::string remove_unit(std::string &str);
template <typename T>
T string_to(const std::string &t, const std::string &unit) {
double tval{0};
try {
tval = std::stod(t);
} catch (const std::invalid_argument &e) {
throw std::runtime_error("Could not convert string to time");
}
using std::chrono::duration;
using std::chrono::duration_cast;
if (unit == "ns") {
return duration_cast<T>(duration<double, std::nano>(tval));
} else if (unit == "us") {
return duration_cast<T>(duration<double, std::micro>(tval));
} else if (unit == "ms") {
return duration_cast<T>(duration<double, std::milli>(tval));
} else if (unit == "s" || unit.empty()) {
return duration_cast<T>(std::chrono::duration<double>(tval));
} else {
throw std::runtime_error(
"Invalid unit in conversion from string to std::chrono::duration");
}
}
// if T has a constructor that takes a string, lets use it.
// template <class T> T string_to(const std::string &arg) { return T{arg}; }
template <typename T> T string_to(const std::string &arg) {
std::string tmp{arg};
auto unit = remove_unit(tmp);
return string_to<T>(tmp, unit);
}
/**
* @brief Convert a string to DetectorType
* @param name string representation of the DetectorType
* @return DetectorType
* @throw runtime_error if the string does not match any DetectorType
*/
template <> DetectorType string_to(const std::string &arg);
/**
* @brief Convert a string to TimingMode
* @param mode string representation of the TimingMode
* @return TimingMode
* @throw runtime_error if the string does not match any TimingMode
*/
template <> TimingMode string_to(const std::string &arg);
/**
* @brief Convert a string to FrameDiscardPolicy
* @param mode string representation of the FrameDiscardPolicy
* @return FrameDiscardPolicy
* @throw runtime_error if the string does not match any FrameDiscardPolicy
*/
template <> FrameDiscardPolicy string_to(const std::string &arg);
/**
* @brief Convert a string to a DACIndex
* @param arg string representation of the dacIndex
* @return DACIndex
* @throw invalid argument error if the string does not match any DACIndex
*/
template <> DACIndex string_to(const std::string &arg);
} // namespace aare

View File

@@ -1,266 +0,0 @@
// SPDX-License-Identifier: MPL-2.0
#include "aare/defs.hpp"
#include <catch2/catch_test_macros.hpp>
#include <string>
#include "to_string.hpp"
using aare::string_to;
TEST_CASE("DetectorType string to enum") {
REQUIRE(string_to<aare::DetectorType>("Generic") ==
aare::DetectorType::Generic);
REQUIRE(string_to<aare::DetectorType>("Eiger") == aare::DetectorType::Eiger);
REQUIRE(string_to<aare::DetectorType>("Gotthard") ==
aare::DetectorType::Gotthard);
REQUIRE(string_to<aare::DetectorType>("Jungfrau") ==
aare::DetectorType::Jungfrau);
REQUIRE(string_to<aare::DetectorType>("ChipTestBoard") ==
aare::DetectorType::ChipTestBoard);
REQUIRE(string_to<aare::DetectorType>("Moench") ==
aare::DetectorType::Moench);
REQUIRE(string_to<aare::DetectorType>("Mythen3") ==
aare::DetectorType::Mythen3);
REQUIRE(string_to<aare::DetectorType>("Gotthard2") ==
aare::DetectorType::Gotthard2);
REQUIRE(string_to<aare::DetectorType>("Xilinx_ChipTestBoard") ==
aare::DetectorType::Xilinx_ChipTestBoard);
REQUIRE(string_to<aare::DetectorType>("Moench03") ==
aare::DetectorType::Moench03);
REQUIRE(string_to<aare::DetectorType>("Moench03_old") ==
aare::DetectorType::Moench03_old);
REQUIRE(string_to<aare::DetectorType>("Unknown") ==
aare::DetectorType::Unknown);
REQUIRE_THROWS(string_to<aare::DetectorType>("invalid_detector"));
}
TEST_CASE("TimingMode string to enum") {
REQUIRE(string_to<aare::TimingMode>("auto") == aare::TimingMode::Auto);
REQUIRE(string_to<aare::TimingMode>("trigger") == aare::TimingMode::Trigger);
REQUIRE_THROWS(string_to<aare::TimingMode>("invalid_mode"));
}
TEST_CASE("FrameDiscardPolicy string to enum") {
REQUIRE(string_to<aare::FrameDiscardPolicy>("nodiscard") ==
aare::FrameDiscardPolicy::NoDiscard);
REQUIRE(string_to<aare::FrameDiscardPolicy>("discard") ==
aare::FrameDiscardPolicy::Discard);
REQUIRE(string_to<aare::FrameDiscardPolicy>("discardpartial") ==
aare::FrameDiscardPolicy::DiscardPartial);
REQUIRE_THROWS(string_to<aare::FrameDiscardPolicy>("invalid_policy"));
}
TEST_CASE("DACIndex string to enum") {
// Numeric DACs
REQUIRE(string_to<aare::DACIndex>("dac 0") == aare::DACIndex::DAC_0);
REQUIRE(string_to<aare::DACIndex>("dac 1") == aare::DACIndex::DAC_1);
REQUIRE(string_to<aare::DACIndex>("dac 2") == aare::DACIndex::DAC_2);
REQUIRE(string_to<aare::DACIndex>("dac 3") == aare::DACIndex::DAC_3);
REQUIRE(string_to<aare::DACIndex>("dac 4") == aare::DACIndex::DAC_4);
REQUIRE(string_to<aare::DACIndex>("dac 5") == aare::DACIndex::DAC_5);
REQUIRE(string_to<aare::DACIndex>("dac 6") == aare::DACIndex::DAC_6);
REQUIRE(string_to<aare::DACIndex>("dac 7") == aare::DACIndex::DAC_7);
REQUIRE(string_to<aare::DACIndex>("dac 8") == aare::DACIndex::DAC_8);
REQUIRE(string_to<aare::DACIndex>("dac 9") == aare::DACIndex::DAC_9);
REQUIRE(string_to<aare::DACIndex>("dac 10") == aare::DACIndex::DAC_10);
REQUIRE(string_to<aare::DACIndex>("dac 11") == aare::DACIndex::DAC_11);
REQUIRE(string_to<aare::DACIndex>("dac 12") == aare::DACIndex::DAC_12);
REQUIRE(string_to<aare::DACIndex>("dac 13") == aare::DACIndex::DAC_13);
REQUIRE(string_to<aare::DACIndex>("dac 14") == aare::DACIndex::DAC_14);
REQUIRE(string_to<aare::DACIndex>("dac 15") == aare::DACIndex::DAC_15);
REQUIRE(string_to<aare::DACIndex>("dac 16") == aare::DACIndex::DAC_16);
REQUIRE(string_to<aare::DACIndex>("dac 17") == aare::DACIndex::DAC_17);
// Named DACs
REQUIRE(string_to<aare::DACIndex>("vsvp") == aare::DACIndex::VSVP);
REQUIRE(string_to<aare::DACIndex>("vtrim") == aare::DACIndex::VTRIM);
REQUIRE(string_to<aare::DACIndex>("vrpreamp") == aare::DACIndex::VRPREAMP);
REQUIRE(string_to<aare::DACIndex>("vrshaper") == aare::DACIndex::VRSHAPER);
REQUIRE(string_to<aare::DACIndex>("vsvn") == aare::DACIndex::VSVN);
REQUIRE(string_to<aare::DACIndex>("vtgstv") == aare::DACIndex::VTGSTV);
REQUIRE(string_to<aare::DACIndex>("vcmp_ll") == aare::DACIndex::VCMP_LL);
REQUIRE(string_to<aare::DACIndex>("vcmp_lr") == aare::DACIndex::VCMP_LR);
REQUIRE(string_to<aare::DACIndex>("vcal") == aare::DACIndex::VCAL);
REQUIRE(string_to<aare::DACIndex>("vcmp_rl") == aare::DACIndex::VCMP_RL);
REQUIRE(string_to<aare::DACIndex>("rxb_rb") == aare::DACIndex::RXB_RB);
REQUIRE(string_to<aare::DACIndex>("rxb_lb") == aare::DACIndex::RXB_LB);
REQUIRE(string_to<aare::DACIndex>("vcmp_rr") == aare::DACIndex::VCMP_RR);
REQUIRE(string_to<aare::DACIndex>("vcp") == aare::DACIndex::VCP);
REQUIRE(string_to<aare::DACIndex>("vcn") == aare::DACIndex::VCN);
REQUIRE(string_to<aare::DACIndex>("vishaper") == aare::DACIndex::VISHAPER);
REQUIRE(string_to<aare::DACIndex>("vthreshold") == aare::DACIndex::VTHRESHOLD);
REQUIRE(string_to<aare::DACIndex>("vref_ds") == aare::DACIndex::VREF_DS);
REQUIRE(string_to<aare::DACIndex>("vout_cm") == aare::DACIndex::VOUT_CM);
REQUIRE(string_to<aare::DACIndex>("vin_cm") == aare::DACIndex::VIN_CM);
REQUIRE(string_to<aare::DACIndex>("vref_comp") == aare::DACIndex::VREF_COMP);
REQUIRE(string_to<aare::DACIndex>("vb_comp") == aare::DACIndex::VB_COMP);
REQUIRE(string_to<aare::DACIndex>("vdd_prot") == aare::DACIndex::VDD_PROT);
REQUIRE(string_to<aare::DACIndex>("vin_com") == aare::DACIndex::VIN_COM);
REQUIRE(string_to<aare::DACIndex>("vref_prech") == aare::DACIndex::VREF_PRECH);
REQUIRE(string_to<aare::DACIndex>("vb_pixbuf") == aare::DACIndex::VB_PIXBUF);
REQUIRE(string_to<aare::DACIndex>("vb_ds") == aare::DACIndex::VB_DS);
REQUIRE(string_to<aare::DACIndex>("vref_h_adc") == aare::DACIndex::VREF_H_ADC);
REQUIRE(string_to<aare::DACIndex>("vb_comp_fe") == aare::DACIndex::VB_COMP_FE);
REQUIRE(string_to<aare::DACIndex>("vb_comp_adc") == aare::DACIndex::VB_COMP_ADC);
REQUIRE(string_to<aare::DACIndex>("vcom_cds") == aare::DACIndex::VCOM_CDS);
REQUIRE(string_to<aare::DACIndex>("vref_rstore") == aare::DACIndex::VREF_RSTORE);
REQUIRE(string_to<aare::DACIndex>("vb_opa_1st") == aare::DACIndex::VB_OPA_1ST);
REQUIRE(string_to<aare::DACIndex>("vref_comp_fe") == aare::DACIndex::VREF_COMP_FE);
REQUIRE(string_to<aare::DACIndex>("vcom_adc1") == aare::DACIndex::VCOM_ADC1);
REQUIRE(string_to<aare::DACIndex>("vref_l_adc") == aare::DACIndex::VREF_L_ADC);
REQUIRE(string_to<aare::DACIndex>("vref_cds") == aare::DACIndex::VREF_CDS);
REQUIRE(string_to<aare::DACIndex>("vb_cs") == aare::DACIndex::VB_CS);
REQUIRE(string_to<aare::DACIndex>("vb_opa_fd") == aare::DACIndex::VB_OPA_FD);
REQUIRE(string_to<aare::DACIndex>("vcom_adc2") == aare::DACIndex::VCOM_ADC2);
REQUIRE(string_to<aare::DACIndex>("vcassh") == aare::DACIndex::VCASSH);
REQUIRE(string_to<aare::DACIndex>("vth2") == aare::DACIndex::VTH2);
REQUIRE(string_to<aare::DACIndex>("vrshaper_n") == aare::DACIndex::VRSHAPER_N);
REQUIRE(string_to<aare::DACIndex>("vipre_out") == aare::DACIndex::VIPRE_OUT);
REQUIRE(string_to<aare::DACIndex>("vth3") == aare::DACIndex::VTH3);
REQUIRE(string_to<aare::DACIndex>("vth1") == aare::DACIndex::VTH1);
REQUIRE(string_to<aare::DACIndex>("vicin") == aare::DACIndex::VICIN);
REQUIRE(string_to<aare::DACIndex>("vcas") == aare::DACIndex::VCAS);
REQUIRE(string_to<aare::DACIndex>("vcal_n") == aare::DACIndex::VCAL_N);
REQUIRE(string_to<aare::DACIndex>("vipre") == aare::DACIndex::VIPRE);
REQUIRE(string_to<aare::DACIndex>("vcal_p") == aare::DACIndex::VCAL_P);
REQUIRE(string_to<aare::DACIndex>("vdcsh") == aare::DACIndex::VDCSH);
REQUIRE(string_to<aare::DACIndex>("vbp_colbuf") == aare::DACIndex::VBP_COLBUF);
REQUIRE(string_to<aare::DACIndex>("vb_sda") == aare::DACIndex::VB_SDA);
REQUIRE(string_to<aare::DACIndex>("vcasc_sfp") == aare::DACIndex::VCASC_SFP);
REQUIRE(string_to<aare::DACIndex>("vipre_cds") == aare::DACIndex::VIPRE_CDS);
REQUIRE(string_to<aare::DACIndex>("ibias_sfp") == aare::DACIndex::IBIAS_SFP);
REQUIRE(string_to<aare::DACIndex>("trimbits") == aare::DACIndex::TRIMBIT_SCAN);
REQUIRE(string_to<aare::DACIndex>("highvoltage") == aare::DACIndex::HIGH_VOLTAGE);
REQUIRE(string_to<aare::DACIndex>("iodelay") == aare::DACIndex::IO_DELAY);
REQUIRE(string_to<aare::DACIndex>("temp_adc") == aare::DACIndex::TEMPERATURE_ADC);
REQUIRE(string_to<aare::DACIndex>("temp_fpga") == aare::DACIndex::TEMPERATURE_FPGA);
REQUIRE(string_to<aare::DACIndex>("temp_fpgaext") == aare::DACIndex::TEMPERATURE_FPGAEXT);
REQUIRE(string_to<aare::DACIndex>("temp_10ge") == aare::DACIndex::TEMPERATURE_10GE);
REQUIRE(string_to<aare::DACIndex>("temp_dcdc") == aare::DACIndex::TEMPERATURE_DCDC);
REQUIRE(string_to<aare::DACIndex>("temp_sodl") == aare::DACIndex::TEMPERATURE_SODL);
REQUIRE(string_to<aare::DACIndex>("temp_sodr") == aare::DACIndex::TEMPERATURE_SODR);
REQUIRE(string_to<aare::DACIndex>("temp_fpgafl") == aare::DACIndex::TEMPERATURE_FPGA2);
REQUIRE(string_to<aare::DACIndex>("temp_fpgafr") == aare::DACIndex::TEMPERATURE_FPGA3);
REQUIRE(string_to<aare::DACIndex>("temp_slowadc") == aare::DACIndex::SLOW_ADC_TEMP);
REQUIRE_THROWS(string_to<aare::DACIndex>("invalid_dac"));
}
TEST_CASE("Remove unit from string") {
using aare::remove_unit;
// Test basic numeric value with unit
{
std::string input = "123.45 V";
std::string unit = remove_unit(input);
REQUIRE(unit == "V");
REQUIRE(input == "123.45");
}
// Test integer value with unit
{
std::string input = "42 Hz";
std::string unit = remove_unit(input);
REQUIRE(unit == "Hz");
REQUIRE(input == "42");
}
// Test negative value with unit
{
std::string input = "-50.5 mV";
std::string unit = remove_unit(input);
REQUIRE(unit == "mV");
REQUIRE(input == "-50.5");
}
// Test value with no unit (only numbers)
{
std::string input = "123.45";
std::string unit = remove_unit(input);
REQUIRE(unit == "");
REQUIRE(input == "123.45");
}
// Test value with only unit (letters at start)
{
std::string input = "kHz";
std::string unit = remove_unit(input);
REQUIRE(unit == "kHz");
REQUIRE(input == "");
}
// Test with multiple word units
{
std::string input = "100 degrees Celsius";
std::string unit = remove_unit(input);
REQUIRE(unit == "degrees Celsius");
REQUIRE(input == "100");
}
// Test with scientific notation
{
std::string input = "1.23e-5 A";
std::string unit = remove_unit(input);
REQUIRE(unit == "A");
REQUIRE(input == "1.23e-5");
}
// Another test with scientific notation
{
std::string input = "-4.56E6 m/s";
std::string unit = remove_unit(input);
REQUIRE(unit == "m/s");
REQUIRE(input == "-4.56E6");
}
// Test with scientific notation uppercase
{
std::string input = "5.67E+3 Hz";
std::string unit = remove_unit(input);
REQUIRE(unit == "Hz");
REQUIRE(input == "5.67E+3");
}
// Test with leading zeros
{
std::string input = "00123 ohm";
std::string unit = remove_unit(input);
REQUIRE(unit == "ohm");
REQUIRE(input == "00123");
}
// Test with leading zeros no space
{
std::string input = "00123ohm";
std::string unit = remove_unit(input);
REQUIRE(unit == "ohm");
REQUIRE(input == "00123");
}
// Test empty string
{
std::string input = "";
std::string unit = remove_unit(input);
REQUIRE(unit == "");
REQUIRE(input == "");
}
}
TEST_CASE("Conversions from time string to chrono durations") {
using namespace std::chrono;
using aare::string_to;
REQUIRE(string_to<nanoseconds>("100 ns") == nanoseconds(100));
REQUIRE(string_to<nanoseconds>("1s") == nanoseconds(1000000000));
REQUIRE(string_to<microseconds>("200 us") == microseconds(200));
REQUIRE(string_to<milliseconds>("300 ms") == milliseconds(300));
REQUIRE(string_to<seconds>("5 s") == seconds(5));
REQUIRE(string_to<nanoseconds>("1.5 us") == nanoseconds(1500));
REQUIRE(string_to<microseconds>("2.5 ms") == microseconds(2500));
REQUIRE(string_to<milliseconds>("3.5 s") == milliseconds(3500));
REQUIRE(string_to<seconds>("2") == seconds(2)); // No unit defaults to seconds
REQUIRE_THROWS(string_to<seconds>("10 min")); // Unsupported unit
}