added custom io for my302

RELEASE.md

@@ -1,12 +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. 
 
 
 ### 2025.11.21

docs/CMakeLists.txt

@@ -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"

docs/src/Interpolation.rst

@@ -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:

docs/src/index.rst

@@ -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::

docs/src/pyCtbRawFile.rst

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

docs/src/pyInterpolation.rst

@@ -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:
+

docs/src/python/cluster/index.rst

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

docs/src/python/cluster/pyInterpolation.rst

@@ -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:
-

docs/src/python/file/index.rst

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

docs/src/python/file/pyCtbRawFile.rst

@@ -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:

docs/src/python/file/pyTransform.rst

@@ -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:

include/aare/CalculateEta.hpp

@@ -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{};

include/aare/Interpolator.hpp

@@ -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,28 +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);
-
-    // 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());
 
@@ -194,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;
@@ -226,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

include/aare/decode.hpp

@@ -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.

include/aare/defs.hpp

@@ -372,15 +372,4 @@ constexpr uint16_t ADC_MASK =
  */
 template <> DACIndex StringTo(const std::string &arg);
 
-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;
-
-};
-
 } // namespace aare

python/CMakeLists.txt

@@ -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
 )
 
 

python/aare/experimental/custom_io.py

@@ -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

python/aare/transform.py

@@ -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()

python/src/ctb_raw_file.hpp

@@ -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",

python/src/interpolation.hpp

@@ -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",

python/src/module.cpp

@@ -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"

python/tests/test_Cluster.py

@@ -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

python/tests/test_Interpolation.py

@@ -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)))
+

python/tests/test_interpolation_simulated.py

@@ -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

src/CtbRawFile.cpp

@@ -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) {

src/Interpolator.cpp

@@ -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++) {

src/NDView.test.cpp

@@ -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});
-
-}

src/decode.cpp

@@ -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

src/decode.test.cpp

@@ -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);
-    }
-
-}

src/defs.cpp

@@ -298,22 +298,4 @@ template <> DACIndex StringTo(const std::string &arg) {
                                     "\"");
 }
 
-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);
-
-}
-
-bool BitOffset::operator==(const BitOffset& other) const {
-        return m_offset == other.m_offset;
-    }
-
-bool BitOffset::operator<(const BitOffset& other) const {
-        return m_offset < other.m_offset;
-    }
-
-
-
 } // namespace aare

src/defs.test.cpp

@@ -83,28 +83,6 @@ TEST_CASE("DynamicCluster creation") {
     REQUIRE(c2.data() != nullptr);
 }
 
-TEST_CASE("Basic ops on BitOffset"){
-    REQUIRE_THROWS(aare::BitOffset(10));
-
-    aare::BitOffset offset(5);
-    REQUIRE(offset.value()==5);
-
-    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);
-}
-
-
-
 // TEST_CASE("cluster set and get data") {
 
 //     aare::DynamicCluster c2(33, 44, aare::Dtype(typeid(double)));