test implementation

- 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()
- Comment for documentation
- Some extra tests

.github/workflows/build_and_deploy_conda.yml

@@ -10,7 +10,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        platform: [ubuntu-latest, ] # macos-12, windows-2019]
+        platform: [ubuntu-latest] # macos-12, windows-2019]
         python-version: ["3.12",]
 
     runs-on: ${{ matrix.platform }}

.github/workflows/build_with_docs.yml

@@ -2,11 +2,15 @@ name: Build the package using cmake then documentation
 
 on:
   workflow_dispatch:
+  push:
   pull_request:
   release:
     types:
       - published
     
+env:
+  # Customize the CMake build type here (Release, Debug, RelWithDebInfo, etc.)
+  BUILD_TYPE: Debug
 
 permissions:
   contents: read
@@ -18,7 +22,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        platform: [ubuntu-latest, ] 
+        platform: [ubuntu-latest, macos-latest] 
         python-version: ["3.12",]
 
     runs-on: ${{ matrix.platform }}
@@ -39,15 +43,20 @@ jobs:
           channels: conda-forge
           conda-remove-defaults: "true"
 
-      - name: Build library 
+      - name: Build library and docs
         run: |
             mkdir build 
             cd build 
-            cmake .. -DAARE_SYSTEM_LIBRARIES=ON -DAARE_PYTHON_BINDINGS=ON -DAARE_DOCS=ON 
-            make -j 2 
+            cmake .. -DCMAKE_BUILD_TYPE=${{env.BUILD_TYPE}} -DAARE_SYSTEM_LIBRARIES=ON -DAARE_PYTHON_BINDINGS=ON -DAARE_DOCS=ON -DAARE_TESTS=ON
+            make -j 4 
             make docs
 
+      - name: C++ unit tests
+        working-directory: ${{github.workspace}}/build
+        run: ctest -C ${{env.BUILD_TYPE}} -j4
+
       - name: Upload static files as artifact
+        if: matrix.platform == 'ubuntu-latest'
         id: deployment
         uses: actions/upload-pages-artifact@v3 
         with:

CMakeLists.txt

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 cmake_minimum_required(VERSION 3.15)
 
 project(aare 
@@ -386,6 +387,7 @@ 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
 ) 
@@ -442,6 +444,7 @@ if(AARE_TESTS)
             ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/DetectorGeometry.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/Interpolation.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/NDArray.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/NDView.test.cpp
@@ -458,6 +461,7 @@ 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} )

LICENSE

@@ -0,0 +1,373 @@
+Mozilla Public License Version 2.0
+==================================
+
+1. Definitions
+--------------
+
+1.1. "Contributor"
+    means each individual or legal entity that creates, contributes to
+    the creation of, or owns Covered Software.
+
+1.2. "Contributor Version"
+    means the combination of the Contributions of others (if any) used
+    by a Contributor and that particular Contributor's Contribution.
+
+1.3. "Contribution"
+    means Covered Software of a particular Contributor.
+
+1.4. "Covered Software"
+    means Source Code Form to which the initial Contributor has attached
+    the notice in Exhibit A, the Executable Form of such Source Code
+    Form, and Modifications of such Source Code Form, in each case
+    including portions thereof.
+
+1.5. "Incompatible With Secondary Licenses"
+    means
+
+    (a) that the initial Contributor has attached the notice described
+        in Exhibit B to the Covered Software; or
+
+    (b) that the Covered Software was made available under the terms of
+        version 1.1 or earlier of the License, but not also under the
+        terms of a Secondary License.
+
+1.6. "Executable Form"
+    means any form of the work other than Source Code Form.
+
+1.7. "Larger Work"
+    means a work that combines Covered Software with other material, in 
+    a separate file or files, that is not Covered Software.
+
+1.8. "License"
+    means this document.
+
+1.9. "Licensable"
+    means having the right to grant, to the maximum extent possible,
+    whether at the time of the initial grant or subsequently, any and
+    all of the rights conveyed by this License.
+
+1.10. "Modifications"
+    means any of the following:
+
+    (a) any file in Source Code Form that results from an addition to,
+        deletion from, or modification of the contents of Covered
+        Software; or
+
+    (b) any new file in Source Code Form that contains any Covered
+        Software.
+
+1.11. "Patent Claims" of a Contributor
+    means any patent claim(s), including without limitation, method,
+    process, and apparatus claims, in any patent Licensable by such
+    Contributor that would be infringed, but for the grant of the
+    License, by the making, using, selling, offering for sale, having
+    made, import, or transfer of either its Contributions or its
+    Contributor Version.
+
+1.12. "Secondary License"
+    means either the GNU General Public License, Version 2.0, the GNU
+    Lesser General Public License, Version 2.1, the GNU Affero General
+    Public License, Version 3.0, or any later versions of those
+    licenses.
+
+1.13. "Source Code Form"
+    means the form of the work preferred for making modifications.
+
+1.14. "You" (or "Your")
+    means an individual or a legal entity exercising rights under this
+    License. For legal entities, "You" includes any entity that
+    controls, is controlled by, or is under common control with You. For
+    purposes of this definition, "control" means (a) the power, direct
+    or indirect, to cause the direction or management of such entity,
+    whether by contract or otherwise, or (b) ownership of more than
+    fifty percent (50%) of the outstanding shares or beneficial
+    ownership of such entity.
+
+2. License Grants and Conditions
+--------------------------------
+
+2.1. Grants
+
+Each Contributor hereby grants You a world-wide, royalty-free,
+non-exclusive license:
+
+(a) under intellectual property rights (other than patent or trademark)
+    Licensable by such Contributor to use, reproduce, make available,
+    modify, display, perform, distribute, and otherwise exploit its
+    Contributions, either on an unmodified basis, with Modifications, or
+    as part of a Larger Work; and
+
+(b) under Patent Claims of such Contributor to make, use, sell, offer
+    for sale, have made, import, and otherwise transfer either its
+    Contributions or its Contributor Version.
+
+2.2. Effective Date
+
+The licenses granted in Section 2.1 with respect to any Contribution
+become effective for each Contribution on the date the Contributor first
+distributes such Contribution.
+
+2.3. Limitations on Grant Scope
+
+The licenses granted in this Section 2 are the only rights granted under
+this License. No additional rights or licenses will be implied from the
+distribution or licensing of Covered Software under this License.
+Notwithstanding Section 2.1(b) above, no patent license is granted by a
+Contributor:
+
+(a) for any code that a Contributor has removed from Covered Software;
+    or
+
+(b) for infringements caused by: (i) Your and any other third party's
+    modifications of Covered Software, or (ii) the combination of its
+    Contributions with other software (except as part of its Contributor
+    Version); or
+
+(c) under Patent Claims infringed by Covered Software in the absence of
+    its Contributions.
+
+This License does not grant any rights in the trademarks, service marks,
+or logos of any Contributor (except as may be necessary to comply with
+the notice requirements in Section 3.4).
+
+2.4. Subsequent Licenses
+
+No Contributor makes additional grants as a result of Your choice to
+distribute the Covered Software under a subsequent version of this
+License (see Section 10.2) or under the terms of a Secondary License (if
+permitted under the terms of Section 3.3).
+
+2.5. Representation
+
+Each Contributor represents that the Contributor believes its
+Contributions are its original creation(s) or it has sufficient rights
+to grant the rights to its Contributions conveyed by this License.
+
+2.6. Fair Use
+
+This License is not intended to limit any rights You have under
+applicable copyright doctrines of fair use, fair dealing, or other
+equivalents.
+
+2.7. Conditions
+
+Sections 3.1, 3.2, 3.3, and 3.4 are conditions of the licenses granted
+in Section 2.1.
+
+3. Responsibilities
+-------------------
+
+3.1. Distribution of Source Form
+
+All distribution of Covered Software in Source Code Form, including any
+Modifications that You create or to which You contribute, must be under
+the terms of this License. You must inform recipients that the Source
+Code Form of the Covered Software is governed by the terms of this
+License, and how they can obtain a copy of this License. You may not
+attempt to alter or restrict the recipients' rights in the Source Code
+Form.
+
+3.2. Distribution of Executable Form
+
+If You distribute Covered Software in Executable Form then:
+
+(a) such Covered Software must also be made available in Source Code
+    Form, as described in Section 3.1, and You must inform recipients of
+    the Executable Form how they can obtain a copy of such Source Code
+    Form by reasonable means in a timely manner, at a charge no more
+    than the cost of distribution to the recipient; and
+
+(b) You may distribute such Executable Form under the terms of this
+    License, or sublicense it under different terms, provided that the
+    license for the Executable Form does not attempt to limit or alter
+    the recipients' rights in the Source Code Form under this License.
+
+3.3. Distribution of a Larger Work
+
+You may create and distribute a Larger Work under terms of Your choice,
+provided that You also comply with the requirements of this License for
+the Covered Software. If the Larger Work is a combination of Covered
+Software with a work governed by one or more Secondary Licenses, and the
+Covered Software is not Incompatible With Secondary Licenses, this
+License permits You to additionally distribute such Covered Software
+under the terms of such Secondary License(s), so that the recipient of
+the Larger Work may, at their option, further distribute the Covered
+Software under the terms of either this License or such Secondary
+License(s).
+
+3.4. Notices
+
+You may not remove or alter the substance of any license notices
+(including copyright notices, patent notices, disclaimers of warranty,
+or limitations of liability) contained within the Source Code Form of
+the Covered Software, except that You may alter any license notices to
+the extent required to remedy known factual inaccuracies.
+
+3.5. Application of Additional Terms
+
+You may choose to offer, and to charge a fee for, warranty, support,
+indemnity or liability obligations to one or more recipients of Covered
+Software. However, You may do so only on Your own behalf, and not on
+behalf of any Contributor. You must make it absolutely clear that any
+such warranty, support, indemnity, or liability obligation is offered by
+You alone, and You hereby agree to indemnify every Contributor for any
+liability incurred by such Contributor as a result of warranty, support,
+indemnity or liability terms You offer. You may include additional
+disclaimers of warranty and limitations of liability specific to any
+jurisdiction.
+
+4. Inability to Comply Due to Statute or Regulation
+---------------------------------------------------
+
+If it is impossible for You to comply with any of the terms of this
+License with respect to some or all of the Covered Software due to
+statute, judicial order, or regulation then You must: (a) comply with
+the terms of this License to the maximum extent possible; and (b)
+describe the limitations and the code they affect. Such description must
+be placed in a text file included with all distributions of the Covered
+Software under this License. Except to the extent prohibited by statute
+or regulation, such description must be sufficiently detailed for a
+recipient of ordinary skill to be able to understand it.
+
+5. Termination
+--------------
+
+5.1. The rights granted under this License will terminate automatically
+if You fail to comply with any of its terms. However, if You become
+compliant, then the rights granted under this License from a particular
+Contributor are reinstated (a) provisionally, unless and until such
+Contributor explicitly and finally terminates Your grants, and (b) on an
+ongoing basis, if such Contributor fails to notify You of the
+non-compliance by some reasonable means prior to 60 days after You have
+come back into compliance. Moreover, Your grants from a particular
+Contributor are reinstated on an ongoing basis if such Contributor
+notifies You of the non-compliance by some reasonable means, this is the
+first time You have received notice of non-compliance with this License
+from such Contributor, and You become compliant prior to 30 days after
+Your receipt of the notice.
+
+5.2. If You initiate litigation against any entity by asserting a patent
+infringement claim (excluding declaratory judgment actions,
+counter-claims, and cross-claims) alleging that a Contributor Version
+directly or indirectly infringes any patent, then the rights granted to
+You by any and all Contributors for the Covered Software under Section
+2.1 of this License shall terminate.
+
+5.3. In the event of termination under Sections 5.1 or 5.2 above, all
+end user license agreements (excluding distributors and resellers) which
+have been validly granted by You or Your distributors under this License
+prior to termination shall survive termination.
+
+************************************************************************
+*                                                                      *
+*  6. Disclaimer of Warranty                                           *
+*  -------------------------                                           *
+*                                                                      *
+*  Covered Software is provided under this License on an "as is"       *
+*  basis, without warranty of any kind, either expressed, implied, or  *
+*  statutory, including, without limitation, warranties that the       *
+*  Covered Software is free of defects, merchantable, fit for a        *
+*  particular purpose or non-infringing. The entire risk as to the     *
+*  quality and performance of the Covered Software is with You.        *
+*  Should any Covered Software prove defective in any respect, You     *
+*  (not any Contributor) assume the cost of any necessary servicing,   *
+*  repair, or correction. This disclaimer of warranty constitutes an   *
+*  essential part of this License. No use of any Covered Software is   *
+*  authorized under this License except under this disclaimer.         *
+*                                                                      *
+************************************************************************
+
+************************************************************************
+*                                                                      *
+*  7. Limitation of Liability                                          *
+*  --------------------------                                          *
+*                                                                      *
+*  Under no circumstances and under no legal theory, whether tort      *
+*  (including negligence), contract, or otherwise, shall any           *
+*  Contributor, or anyone who distributes Covered Software as          *
+*  permitted above, be liable to You for any direct, indirect,         *
+*  special, incidental, or consequential damages of any character      *
+*  including, without limitation, damages for lost profits, loss of    *
+*  goodwill, work stoppage, computer failure or malfunction, or any    *
+*  and all other commercial damages or losses, even if such party      *
+*  shall have been informed of the possibility of such damages. This   *
+*  limitation of liability shall not apply to liability for death or   *
+*  personal injury resulting from such party's negligence to the       *
+*  extent applicable law prohibits such limitation. Some               *
+*  jurisdictions do not allow the exclusion or limitation of           *
+*  incidental or consequential damages, so this exclusion and          *
+*  limitation may not apply to You.                                    *
+*                                                                      *
+************************************************************************
+
+8. Litigation
+-------------
+
+Any litigation relating to this License may be brought only in the
+courts of a jurisdiction where the defendant maintains its principal
+place of business and such litigation shall be governed by laws of that
+jurisdiction, without reference to its conflict-of-law provisions.
+Nothing in this Section shall prevent a party's ability to bring
+cross-claims or counter-claims.
+
+9. Miscellaneous
+----------------
+
+This License represents the complete agreement concerning the subject
+matter hereof. If any provision of this License is held to be
+unenforceable, such provision shall be reformed only to the extent
+necessary to make it enforceable. Any law or regulation which provides
+that the language of a contract shall be construed against the drafter
+shall not be used to construe this License against a Contributor.
+
+10. Versions of the License
+---------------------------
+
+10.1. New Versions
+
+Mozilla Foundation is the license steward. Except as provided in Section
+10.3, no one other than the license steward has the right to modify or
+publish new versions of this License. Each version will be given a
+distinguishing version number.
+
+10.2. Effect of New Versions
+
+You may distribute the Covered Software under the terms of the version
+of the License under which You originally received the Covered Software,
+or under the terms of any subsequent version published by the license
+steward.
+
+10.3. Modified Versions
+
+If you create software not governed by this License, and you want to
+create a new license for such software, you may create and use a
+modified version of this License if you rename the license and remove
+any references to the name of the license steward (except to note that
+such modified license differs from this License).
+
+10.4. Distributing Source Code Form that is Incompatible With Secondary
+Licenses
+
+If You choose to distribute Source Code Form that is Incompatible With
+Secondary Licenses under the terms of this version of the License, the
+notice described in Exhibit B of this License must be attached.
+
+Exhibit A - Source Code Form License Notice
+-------------------------------------------
+
+  This Source Code Form is subject to the terms of the Mozilla Public
+  License, v. 2.0. If a copy of the MPL was not distributed with this
+  file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+If it is not possible or desirable to put the notice in a particular
+file, then You may include the notice in a location (such as a LICENSE
+file in a relevant directory) where a recipient would be likely to look
+for such a notice.
+
+You may add additional accurate notices of copyright ownership.
+
+Exhibit B - "Incompatible With Secondary Licenses" Notice
+---------------------------------------------------------
+
+  This Source Code Form is "Incompatible With Secondary Licenses", as
+  defined by the Mozilla Public License, v. 2.0.

README.md

@@ -1,6 +1,14 @@
 # aare
 Data analysis library for PSI hybrid detectors
 
+## Documentation 
+
+Detailed documentation including installation can be found in [Documentation](https://slsdetectorgroup.github.io/aare/)
+
+## License
+
+This project is licensed under the MPL-2.0 license.
+See the LICENSE file or https://www.mozilla.org/en-US/MPL/ for details.
 
 ## Build and install
 

RELEASE.md

@@ -1,11 +1,52 @@
 # Release notes
 
-### 2025.10.1
+## head
 
-Bugfixes: 
+### 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
+
+### New Features: 
+
+- Added SPDX-License-Identifier: MPL-2.0 to source files
+- Calculate Eta3 supports all cluster types 
+- interpolation class supports using cross eta3x3 and eta3x3 on full cluster as well as eta2x2 on full cluster
+- interpolation class has option to calculate the rosenblatt transform 
+- reduction operations to reduce Clusters of general size to 2x2 or 3x3 clusters 
+- `max_sum_2x2` including index of subcluster with highest energy is now available from Python API 
+- interpolation supports bilinear interpolation of eta values for more fine grained transformed uniform coordinates
+- Interpolation is documented 
+
+- Added tell to ClusterFile. Returns position in bytes for debugging
+
+### Resolved Features: 
+
+- calculate_eta coincides with theoretical definition
+
+### Bugfixes: 
+
+- eta calculation assumes correct photon center 
+- eta transformation to uniform coordinates starts at 0
+- Bug in interpolation 
 - File supports reading new master json file format (multiple ROI's not supported yet)
 
+
+### API Changes: 
+
+- ClusterFinder for 2x2 Cluster disabled 
+- eta stores corner as enum class cTopLeft, cTopRight, BottomLeft, cBottomRight indicating 2x2 subcluster with largest energy relative to cluster center 
+- max_sum_2x2 returns corner as index 
+
 ### 2025.8.22
 
 Features:
@@ -45,6 +86,31 @@ Bugfixes:
 
 - Fixed crash when opening raw files with large number of data files
 
+## Download, Documentation & Support 
+
+### Download
+
+The Source Code: 
+https://github.com/slsdetectorgroup/aare
+
+
+### Documentation 
+
+
+Documentation including installation details: 
+https://github.com/slsdetectorgroup/aare 
+
+
+### Support
+
+
+erik.frojdh@psi.ch \
+alice.mazzoleni@psi.ch \
+dhanya.thattil@psi.ch
+
+
+
+
 
 
 

VERSION

@@ -1 +1 @@
-2025.8.22
+2025.11.21

benchmarks/CMakeLists.txt

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 
 include(FetchContent)
 

benchmarks/calculateeta_benchmark.cpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #include "aare/CalculateEta.hpp"
 #include "aare/ClusterFile.hpp"
 #include <benchmark/benchmark.h>
@@ -8,6 +9,7 @@ class ClusterFixture : public benchmark::Fixture {
   public:
     Cluster<int, 2, 2> cluster_2x2{};
     Cluster<int, 3, 3> cluster_3x3{};
+    Cluster<int, 4, 4> cluster_4x4{};
 
   private:
     using benchmark::Fixture::SetUp;
@@ -26,6 +28,13 @@ class ClusterFixture : public benchmark::Fixture {
 
         cluster_3x3.x = 0;
         cluster_3x3.y = 0;
+
+        int temp_data3[16] = {1, 2,  3,  4,  5,  6,  7,  8,
+                              9, 10, 11, 12, 13, 14, 15, 16};
+        std::copy(std::begin(temp_data3), std::end(temp_data3),
+                  std::begin(cluster_4x4.data));
+        cluster_4x4.x = 0;
+        cluster_4x4.y = 0;
     }
 
     // void TearDown(::benchmark::State& state) {
@@ -67,4 +76,29 @@ BENCHMARK_F(ClusterFixture, CalculateGeneralEtaFor3x3Cluster)
         benchmark::DoNotOptimize(eta);
     }
 }
+
+BENCHMARK_F(ClusterFixture, Calculate2x2Etawithreduction)
+(benchmark::State &st) {
+    for (auto _ : st) {
+        // This code gets timed
+        auto reduced_cluster = reduce_to_2x2(cluster_4x4);
+        Eta2 eta = calculate_eta2(reduced_cluster);
+        auto reduced_cluster_from_3x3 = reduce_to_2x2(cluster_3x3);
+        Eta2 eta2 = calculate_eta2(reduced_cluster_from_3x3);
+        benchmark::DoNotOptimize(eta);
+        benchmark::DoNotOptimize(eta2);
+    }
+}
+
+BENCHMARK_F(ClusterFixture, Calculate2x2Etawithoutreduction)
+(benchmark::State &st) {
+    for (auto _ : st) {
+        // This code gets timed
+        Eta2 eta = calculate_eta2(cluster_4x4);
+        Eta2 eta2 = calculate_eta2(cluster_3x3);
+        benchmark::DoNotOptimize(eta);
+        benchmark::DoNotOptimize(eta2);
+    }
+}
+
 // BENCHMARK_MAIN();

benchmarks/ndarray_benchmark.cpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #include "aare/NDArray.hpp"
 #include <benchmark/benchmark.h>
 

benchmarks/reduce_benchmark.cpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #include "aare/Cluster.hpp"
 #include <benchmark/benchmark.h>
 
@@ -33,8 +34,8 @@ class ClustersForReduceFixture : public benchmark::Fixture {
 };
 
 template <typename T>
-Cluster<T, 3, 3, int16_t> reduce_to_3x3(const Cluster<T, 5, 5, int16_t> &c) {
-    Cluster<T, 3, 3, int16_t> result;
+Cluster<T, 3, 3, uint16_t> reduce_to_3x3(const Cluster<T, 5, 5, uint16_t> &c) {
+    Cluster<T, 3, 3, uint16_t> result;
 
     // Write out the sums in the hope that the compiler can optimize this
     std::array<T, 9> sum_3x3_subclusters;
@@ -140,7 +141,7 @@ Cluster<T, 3, 3, int16_t> reduce_to_3x3(const Cluster<T, 5, 5, int16_t> &c) {
 BENCHMARK_F(ClustersForReduceFixture, Reduce2x2)(benchmark::State &st) {
     for (auto _ : st) {
         // This code gets timed
-        benchmark::DoNotOptimize(reduce_to_2x2<int, 3, 3, int16_t>(
+        benchmark::DoNotOptimize(reduce_to_2x2<int, 3, 3, uint16_t>(
             cluster_3x3)); // make sure compiler evaluates the expression
     }
 }
@@ -156,7 +157,7 @@ BENCHMARK_F(ClustersForReduceFixture, Reduce3x3)(benchmark::State &st) {
     for (auto _ : st) {
         // This code gets timed
         benchmark::DoNotOptimize(
-            reduce_to_3x3<int, 5, 5, int16_t>(cluster_5x5));
+            reduce_to_3x3<int, 5, 5, uint16_t>(cluster_5x5));
     }
 }
 

conda-recipe/conda_build_config.yaml

@@ -1,8 +1,7 @@
 python:
-    # - 3.11
-    # - 3.12
-    # - 3.13
-    - 3.14
+    - 3.11
+    - 3.12
+    - 3.13
 
 c_compiler:
   - gcc                        # [linux]

conda-recipe/meta.yaml

@@ -25,7 +25,7 @@ requirements:
   host:
     - python
     - pip
-    - numpy=2.3
+    - numpy=2.1
     - scikit-build-core
     - pybind11 >=2.13.0
     - matplotlib # needed in host to solve the environment for run
@@ -42,11 +42,12 @@ test:
     - aare
   requires:
     - pytest
-    # - boost-histogram 
+    - boost-histogram 
   source_files:
-    # - python/tests
+    - python/tests
   commands:
-    # - python -m pytest python/tests
+    - python -m pytest python/tests
 
 about:
+  license: SPDX-License-Identifier MPL-2.0
   summary: Data analysis library for hybrid pixel detectors from PSI 

docs/CMakeLists.txt

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 find_package(Doxygen REQUIRED)
 find_package(Sphinx REQUIRED)
 
@@ -11,15 +12,19 @@ set(SPHINX_SOURCE ${CMAKE_CURRENT_SOURCE_DIR}/src)
 set(SPHINX_BUILD ${CMAKE_CURRENT_BINARY_DIR})
 
 
-file(GLOB SPHINX_SOURCE_FILES CONFIGURE_DEPENDS "src/*.rst")
+file(GLOB_RECURSE SPHINX_SOURCE_FILES
+     CONFIGURE_DEPENDS
+     RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}/src"
+     "${CMAKE_CURRENT_SOURCE_DIR}/src/*.rst"
+)
 
+foreach(relpath IN LISTS SPHINX_SOURCE_FILES)
+    set(src  "${CMAKE_CURRENT_SOURCE_DIR}/src/${relpath}")
+    set(dst  "${SPHINX_BUILD}/src/${relpath}")
 
-
-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})
+    message(STATUS "Copying ${src} to ${dst}")
+    configure_file("${src}" "${dst}" COPYONLY)
+endforeach()
 
 configure_file(
     "${CMAKE_CURRENT_SOURCE_DIR}/conf.py.in"
@@ -27,6 +32,8 @@ configure_file(
     @ONLY
 )
 
+file(COPY "${CMAKE_CURRENT_SOURCE_DIR}/figures"
+     DESTINATION "${SPHINX_BUILD}") 
 
 configure_file(
     "${CMAKE_CURRENT_SOURCE_DIR}/static/extra.css"

docs/src/Cluster.rst

@@ -0,0 +1,15 @@
+Cluster
+========
+
+.. doxygenstruct:: aare::Cluster
+   :members:
+   :undoc-members:
+   :private-members:
+
+
+**Free Functions:** 
+
+.. doxygenfunction:: aare::reduce_to_3x3(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>&)
+   
+.. doxygenfunction:: aare::reduce_to_2x2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>&)
+   

docs/src/Interpolation.rst

@@ -0,0 +1,160 @@
+.. _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.
+
+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. 
+
+
+:math:`\eta`-Functions: 
+---------------------------
+
+.. doxygenstruct:: aare::Eta2
+    :members: 
+    :undoc-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: 
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. 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.
+
+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
+    :align: center
+    :alt: Eta2x2 Full 
+
+.. 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}}
+    \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
+    :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=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}}
+    \end{equation*}
+
+The :math:`\eta` values can range between -0.5,0.5.
+
+.. doxygenfunction:: aare::calculate_eta3(const ClusterVector<ClusterType>&)
+
+.. doxygenfunction:: aare::calculate_eta3(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>&)
+
+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. 
+
+.. doxygenfunction:: aare::calculate_cross_eta3(const ClusterVector<ClusterType>&)
+
+.. doxygenfunction:: aare::calculate_cross_eta3(const Cluster<T, ClusterSizeX, ClusterSizeY, 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:
+   :private-members:
+
+

docs/src/index.rst

@@ -22,19 +22,14 @@ AARE
 
 .. toctree::
     :caption: Python API
-    :maxdepth: 1
-    
-    pyFile
-    pycalibration
-    pyCtbRawFile
-    pyClusterFile
-    pyClusterVector
-    pyJungfrauDataFile
-    pyRawFile
-    pyRawMasterFile
-    pyVarClusterFinder
+    :maxdepth: 3
+    :hidden:
 
+    pycalibration
+    python/cluster/index
+    python/file/index
     pyFit
+    
 
 
 .. toctree::
@@ -47,10 +42,12 @@ AARE
     Frame
     File
     Dtype
+    Cluster
     ClusterFinder
     ClusterFinderMT
     ClusterFile
     ClusterVector
+    Interpolation
     JungfrauDataFile
     Pedestal
     RawFile

docs/src/pyCtbRawFile.rst

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

docs/src/python/cluster/index.rst

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

docs/src/python/cluster/pyCluster.rst

@@ -0,0 +1,23 @@
+Cluster
+========
+
+.. py:currentmodule:: aare
+
+.. autoclass:: Cluster
+    :members:
+    :undoc-members:
+    :inherited-members:
+
+
+Below is the API of a cluster of size :math:`3\times 3` and type ``int`` but all variants share the same API.
+
+    .. autoclass::  aare._aare.Cluster3x3i
+        :special-members: __init__
+        :members:
+        :undoc-members:
+        :show-inheritance:
+        :inherited-members:
+
+.. note:: 
+    More functions can be found in the :ref:`ClusterVector <py_clustervector>` documentation. Generally apply functions directly on the ``ClusterVector`` instead of looping over individual clusters.
+

docs/src/python/cluster/pyClusterVector.rst

@@ -1,3 +1,5 @@
+.. _py_clustervector:
+    
 ClusterVector
 ================
 
@@ -28,6 +30,13 @@ C++ functions that support the ClusterVector or to view it as a numpy array.
 
 .. py:currentmodule:: aare
 
+.. autoclass:: ClusterVector
+    :members:
+    :undoc-members:
+    :inherited-members:
+
+Below is the API of the ClusterVector_Cluster3x3i but all variants share the same API.
+
 .. autoclass::  aare._aare.ClusterVector_Cluster3x3i
     :special-members: __init__
     :members:

docs/src/python/cluster/pyInterpolation.rst

@@ -0,0 +1,124 @@
+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

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

docs/src/python/file/pyCtbRawFile.rst

@@ -0,0 +1,25 @@
+
+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

@@ -0,0 +1,27 @@
+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:

etc/add_license.py

@@ -0,0 +1,103 @@
+#!/usr/bin/env python3
+import argparse
+import fnmatch
+import os
+from pathlib import Path
+
+CPP_PATTERNS = ["*.h", "*.hpp", "*.cpp"]
+PY_PATTERNS = ["*.py"]
+CMAKE_PATTERNS = ["CMakeLists.txt"]
+
+FILE_PATTERNS = CPP_PATTERNS + PY_PATTERNS + CMAKE_PATTERNS
+LICENSE_TEXT = "SPDX-License-Identifier: MPL-2.0"
+
+
+def get_comment_prefix(filename: str) -> str | None:
+    if any(fnmatch.fnmatch(filename, p) for p in CPP_PATTERNS):
+        return "// "
+    if any(fnmatch.fnmatch(filename, p) for p in (PY_PATTERNS + CMAKE_PATTERNS)):
+        return "# "
+    return None
+
+
+def matches_pattern(filename: str) -> bool:
+    return any(fnmatch.fnmatch(filename, p) for p in FILE_PATTERNS)
+
+
+def process_file(filepath: Path) -> bool:
+    filename = filepath.name
+    prefix = get_comment_prefix(filename)
+    if not prefix:
+        return False
+
+    license_line = f"{prefix}{LICENSE_TEXT}\n"
+
+    try:
+        with filepath.open("r", encoding="utf-8") as f:
+            lines = f.readlines()
+    except Exception as e:
+        print(f"⚠️ Error reading {filepath}: {e}")
+        return False
+
+    # Skip if SPDX already present anywhere in the file
+    if any("SPDX-License-Identifier" in line for line in lines):
+        return False
+
+    insert_index = 0
+
+    # For Python, keep shebang on the very first line
+    if filename.endswith(".py") and lines:
+        if lines[0].startswith("#!"):
+            insert_index = 1
+
+    lines.insert(insert_index, license_line)
+
+    try:
+        with filepath.open("w", encoding="utf-8") as f:
+            f.writelines(lines)
+    except Exception as e:
+        print(f"⚠️ Error writing {filepath}: {e}")
+        return False
+
+    return True
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser(
+        description="Add SPDX-License-Identifier: MPL-2.0 to source files."
+    )
+    parser.add_argument(
+        "path",
+        help="Root directory to recursively process "
+             "(*.h, *.cpp, *.py, and CMakeLists.txt).",
+    )
+
+    args = parser.parse_args()
+
+    root_path = Path(args.path).expanduser().resolve()
+
+    if not root_path.exists():
+        print(f"Error: Path does not exist: {root_path}")
+        raise SystemExit(1)
+
+    if not root_path.is_dir():
+        print(f"Error: Path is not a directory: {root_path}")
+        raise SystemExit(1)
+
+    print(f"Processing directory: {root_path}")
+    modified = 0
+
+    for dirpath, _, files in os.walk(root_path):
+        dirpath = Path(dirpath)
+        for name in files:
+            if matches_pattern(name):
+                fullpath = dirpath / name
+                if process_file(fullpath):
+                    print(f"✔ Added SPDX: {fullpath}")
+                    modified += 1
+
+    print(f"\nDone. Updated {modified} file(s).")
+
+
+if __name__ == "__main__":
+    main()

etc/dev-env.yml

@@ -3,6 +3,7 @@ channels:
   - conda-forge
 dependencies:
   - anaconda-client
+  - catch2
   - conda-build
   - doxygen 
   - sphinx

etc/update_version.py

@@ -1,4 +1,4 @@
-# SPDX-License-Identifier: LGPL-3.0-or-other
+# SPDX-License-Identifier: MPL-2.0
 # Copyright (C) 2021 Contributors to the Aare Package
 """
 Script to update VERSION file with semantic versioning if provided as an argument, or with 0.0.0 if no argument is provided.
@@ -8,11 +8,12 @@ import sys
 import os
 import re
 from datetime import datetime
+from pathlib import Path
 
 from packaging.version import Version, InvalidVersion
 
 
-SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
+SCRIPT_DIR = Path(__file__).absolute().parent.parent
 
 def is_integer(value): 
     try:
@@ -46,7 +47,8 @@ def get_version():
     
 
 def write_version_to_file(version):
-    version_file_path = os.path.join(SCRIPT_DIR, "VERSION")
+    version_file_path =  SCRIPT_DIR/"VERSION"
+    print(version_file_path)
     with open(version_file_path, "w") as version_file:
         version_file.write(version)
     print(f"Version {version} written to VERSION file.")

include/aare/ArrayExpr.hpp

@@ -1,12 +1,53 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include "aare/defs.hpp"
 #include <array>
 #include <cassert>
 #include <cstddef>
 #include <cstdint>
+#include <type_traits>
 
 namespace aare {
 
+template <ssize_t Dim = 0, typename Strides>
+ssize_t element_offset(const Strides & /*unused*/) {
+    return 0;
+}
+
+template <ssize_t Dim = 0, typename Strides, typename... Ix>
+ssize_t element_offset(const Strides &strides, ssize_t i, Ix... index) {
+    return i * strides[Dim] + element_offset<Dim + 1>(strides, index...);
+}
+
+template <typename Derived, typename T, ssize_t Ndim>
+class NDIndexOps {
+  public:
+    template <typename... Ix>
+    std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) {
+        return derived().data()[element_offset(derived().strides(), index...)];
+    }
+
+    template <typename... Ix>
+    std::enable_if_t<sizeof...(Ix) == Ndim, const T &> operator()(Ix... index) const {
+        return derived().data()[element_offset(derived().strides(), index...)];
+    }
+
+    T &operator()(ssize_t i) {
+        return derived().data()[i];
+    }
+    
+    const T &operator()(ssize_t i) const {
+        return derived().data()[i];
+    }
+
+    T &operator[](ssize_t i) { return derived().data()[i]; }
+    const T &operator[](ssize_t i) const { return derived().data()[i]; }
+
+  private:
+    Derived &derived() { return static_cast<Derived &>(*this); }
+    const Derived &derived() const { return static_cast<const Derived &>(*this); }
+};
+
 template <typename E, ssize_t Ndim> class ArrayExpr {
   public:
     static constexpr bool is_leaf = false;
@@ -95,4 +136,4 @@ auto operator/(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
     return ArrayDiv<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
 }
 
-} // namespace aare
+} // namespace aare

include/aare/CalculateEta.hpp

@@ -1,18 +1,13 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 
 #include "aare/Cluster.hpp"
 #include "aare/ClusterVector.hpp"
 #include "aare/NDArray.hpp"
+#include "aare/defs.hpp"
 
 namespace aare {
 
-enum class corner : int {
-    cTopLeft = 0,
-    cTopRight = 1,
-    cBottomLeft = 2,
-    cBottomRight = 3
-};
-
 enum class pixel : int {
     pBottomLeft = 0,
     pBottom = 1,
@@ -25,37 +20,125 @@ enum class pixel : int {
     pTopRight = 8
 };
 
+// TODO: better to have sum after x,y
+/**
+ * eta struct
+ */
 template <typename T> struct Eta2 {
-    double x;
-    double y;
-    int c{0};
-    T sum;
+    /// @brief eta in x direction
+    double x{};
+    /// @brief eta in y direction
+    double y{};
+    /// @brief index of subcluster with highest energy value (given as corner
+    /// relative to cluster center)
+    corner c{0};
+    /// @brief photon energy (cluster sum)
+    T sum{};
 };
 
 /**
- * @brief Calculate the eta2 values for all clusters in a Clustervector
+ * @brief Calculate the eta2 values for all clusters in a ClusterVector
  */
 template <typename ClusterType,
           typename = std::enable_if_t<is_cluster_v<ClusterType>>>
-NDArray<double, 2> calculate_eta2(const ClusterVector<ClusterType> &clusters) {
-    NDArray<double, 2> eta2({static_cast<int64_t>(clusters.size()), 2});
+std::vector<Eta2<typename ClusterType::value_type>>
+calculate_eta2(const ClusterVector<ClusterType> &clusters) {
+
+    std::vector<Eta2<typename ClusterType::value_type>> eta2{};
+    eta2.reserve(clusters.size());
 
     for (size_t i = 0; i < clusters.size(); i++) {
         auto e = calculate_eta2(clusters[i]);
-        eta2(i, 0) = e.x;
-        eta2(i, 1) = e.y;
+        eta2.push_back(e);
     }
 
     return eta2;
 }
 
+/**
+ * @brief Calculate the full eta2 values for all clusters in a ClusterVector
+ */
+template <typename ClusterType,
+          typename = std::enable_if_t<is_cluster_v<ClusterType>>>
+std::vector<Eta2<typename ClusterType::value_type>>
+calculate_full_eta2(const ClusterVector<ClusterType> &clusters) {
+    std::vector<Eta2<typename ClusterType::value_type>> eta2{};
+    eta2.reserve(clusters.size());
+
+    for (size_t i = 0; i < clusters.size(); i++) {
+        auto e = calculate_full_eta2(clusters[i]);
+        eta2.push_back(e);
+    }
+
+    return eta2;
+}
+
+/**
+ * @brief Calculate eta3 for all 3x3 clusters in a ClusterVector
+ */
+template <typename ClusterType,
+          typename = std::enable_if_t<is_cluster_v<ClusterType>>>
+std::vector<Eta2<typename ClusterType::value_type>>
+calculate_eta3(const ClusterVector<ClusterType> &clusters) {
+    std::vector<Eta2<typename ClusterType::value_type>> eta2{};
+    eta2.reserve(clusters.size());
+
+    for (size_t i = 0; i < clusters.size(); i++) {
+        auto e = calculate_eta3(clusters[i]);
+        eta2.push_back(e);
+    }
+
+    return eta2;
+}
+
+/**
+ * @brief Calculate cross eta3 for all 3x3 clusters in a ClusterVector
+ */
+template <typename ClusterType,
+          typename = std::enable_if_t<is_cluster_v<ClusterType>>>
+std::vector<Eta2<typename ClusterType::value_type>>
+calculate_cross_eta3(const ClusterVector<ClusterType> &clusters) {
+    std::vector<Eta2<typename ClusterType::value_type>> eta2{};
+    eta2.reserve(clusters.size());
+
+    for (size_t i = 0; i < clusters.size(); i++) {
+        auto e = calculate_cross_eta3(clusters[i]);
+        eta2.push_back(e);
+    }
+
+    return eta2;
+}
+
+/**
+ * @brief helper function to calculate eta2 x and y values
+ * @param eta reference to the Eta2 object to update
+ * @param left_x value of the left pixel
+ * @param right_x value of the right pixel
+ * @param bottom_y value of the bottom pixel
+ * @param top_y value of the top pixel
+ */
+template <typename T>
+inline void calculate_eta2(Eta2<T> &eta, const T left_x, const T right_x,
+                           const T bottom_y, const T top_y) {
+    if ((right_x + left_x) != 0)
+        eta.x = static_cast<double>(right_x) /
+                static_cast<double>(right_x + left_x); // between (0,1) the
+                                                       // closer to zero left
+                                                       // value probably larger
+    if ((top_y + bottom_y) != 0)
+        eta.y = static_cast<double>(top_y) /
+                static_cast<double>(top_y + bottom_y); // between (0,1) the
+                                                       // closer to zero bottom
+                                                       // value probably larger
+}
+
 /**
  * @brief Calculate the eta2 values for a generic sized cluster and return them
- * in a Eta2 struct containing etay, etax and the index of the respective 2x2
- * subcluster.
+ * in a Eta2 struct containing etay, etax and the index (as corner) of the
+ * respective 2x2 subcluster relative to the cluster center.
  */
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType>
+          typename CoordType = uint16_t>
 Eta2<T>
 calculate_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
 
@@ -66,73 +149,42 @@ calculate_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
         (ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
 
     auto max_sum = cl.max_sum_2x2();
-    eta.sum = max_sum.first;
-    int c = max_sum.second;
+    eta.sum = max_sum.sum;
+    corner c = max_sum.index;
 
     // subcluster top right from center
-    switch (static_cast<corner>(c)) {
+    switch (c) {
     case (corner::cTopLeft):
-        if ((cl.data[cluster_center_index - 1] +
-             cl.data[cluster_center_index]) != 0)
-            eta.x = static_cast<double>(cl.data[cluster_center_index - 1]) /
-                    static_cast<double>(cl.data[cluster_center_index - 1] +
-                                        cl.data[cluster_center_index]);
-        if ((cl.data[cluster_center_index - ClusterSizeX] +
-             cl.data[cluster_center_index]) != 0)
-            eta.y = static_cast<double>(
-                        cl.data[cluster_center_index - ClusterSizeX]) /
-                    static_cast<double>(
-                        cl.data[cluster_center_index - ClusterSizeX] +
-                        cl.data[cluster_center_index]);
-
-        // dx = 0
-        // dy = 0
+        calculate_eta2(eta, cl.data[cluster_center_index - 1],
+                       cl.data[cluster_center_index],
+                       cl.data[cluster_center_index - ClusterSizeX],
+                       cl.data[cluster_center_index]);
+        // dx = -1
+        // dy = -1
         break;
     case (corner::cTopRight):
-        if (cl.data[cluster_center_index] + cl.data[cluster_center_index + 1] !=
-            0)
-            eta.x = static_cast<double>(cl.data[cluster_center_index]) /
-                    static_cast<double>(cl.data[cluster_center_index] +
-                                        cl.data[cluster_center_index + 1]);
-        if ((cl.data[cluster_center_index - ClusterSizeX] +
-             cl.data[cluster_center_index]) != 0)
-            eta.y = static_cast<double>(
-                        cl.data[cluster_center_index - ClusterSizeX]) /
-                    static_cast<double>(
-                        cl.data[cluster_center_index - ClusterSizeX] +
-                        cl.data[cluster_center_index]);
-        // dx = 1
-        // dy = 0
+        calculate_eta2(eta, cl.data[cluster_center_index],
+                       cl.data[cluster_center_index + 1],
+                       cl.data[cluster_center_index - ClusterSizeX],
+                       cl.data[cluster_center_index]);
+        // dx = 0
+        // dy = -1
         break;
     case (corner::cBottomLeft):
-        if ((cl.data[cluster_center_index - 1] +
-             cl.data[cluster_center_index]) != 0)
-            eta.x = static_cast<double>(cl.data[cluster_center_index - 1]) /
-                    static_cast<double>(cl.data[cluster_center_index - 1] +
-                                        cl.data[cluster_center_index]);
-        if ((cl.data[cluster_center_index] +
-             cl.data[cluster_center_index + ClusterSizeX]) != 0)
-            eta.y = static_cast<double>(cl.data[cluster_center_index]) /
-                    static_cast<double>(
-                        cl.data[cluster_center_index] +
-                        cl.data[cluster_center_index + ClusterSizeX]);
-        // dx = 0
-        // dy = 1
+        calculate_eta2(eta, cl.data[cluster_center_index - 1],
+                       cl.data[cluster_center_index],
+                       cl.data[cluster_center_index],
+                       cl.data[cluster_center_index + ClusterSizeX]);
+        // dx = -1
+        // dy = 0
         break;
     case (corner::cBottomRight):
-        if (cl.data[cluster_center_index] + cl.data[cluster_center_index + 1] !=
-            0)
-            eta.x = static_cast<double>(cl.data[cluster_center_index]) /
-                    static_cast<double>(cl.data[cluster_center_index] +
-                                        cl.data[cluster_center_index + 1]);
-        if ((cl.data[cluster_center_index] +
-             cl.data[cluster_center_index + ClusterSizeX]) != 0)
-            eta.y = static_cast<double>(cl.data[cluster_center_index]) /
-                    static_cast<double>(
-                        cl.data[cluster_center_index] +
-                        cl.data[cluster_center_index + ClusterSizeX]);
-        // dx = 1
-        // dy = 1
+        calculate_eta2(eta, cl.data[cluster_center_index],
+                       cl.data[cluster_center_index + 1],
+                       cl.data[cluster_center_index],
+                       cl.data[cluster_center_index + ClusterSizeX]);
+        // dx = 0
+        // dy = 0
         break;
     }
 
@@ -141,69 +193,255 @@ calculate_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
     return eta;
 }
 
-// TODO! Look up eta2 calculation - photon center should be bottom right corner
-template <typename T>
-Eta2<T> calculate_eta2(const Cluster<T, 2, 2, int16_t> &cl) {
+/**
+ * @brief Calculate the eta2 values for a generic sized cluster and return them
+ * in a Eta2 struct containing etay, etax and the index (as corner) of the
+ * respective 2x2 subcluster relative to the cluster center.
+ */
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType>
+Eta2<T> calculate_full_eta2(
+    const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
+
+    static_assert(ClusterSizeX > 1 && ClusterSizeY > 1);
     Eta2<T> eta{};
 
-    if ((cl.data[0] + cl.data[1]) != 0)
-        eta.x = static_cast<double>(cl.data[2]) /
-                (cl.data[2] + cl.data[3]); // between (0,1) the closer to zero
-                                           // left value probably larger
-    if ((cl.data[0] + cl.data[2]) != 0)
-        eta.y = static_cast<double>(cl.data[1]) /
-                (cl.data[1] + cl.data[3]); // between (0,1) the closer to zero
-                                           // bottom value probably larger
+    constexpr size_t cluster_center_index =
+        (ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
+
+    auto max_sum = cl.max_sum_2x2();
+    eta.sum = max_sum.sum;
+    corner c = max_sum.index;
+
+    // subcluster top right from center
+    switch (c) {
+    case (corner::cTopLeft):
+        if (eta.sum != 0) {
+            eta.x = static_cast<double>(
+                        cl.data[cluster_center_index] +
+                        cl.data[cluster_center_index - ClusterSizeX]) /
+                    static_cast<double>(eta.sum);
+
+            eta.y = static_cast<double>(cl.data[cluster_center_index - 1] +
+                                        cl.data[cluster_center_index]) /
+                    static_cast<double>(eta.sum);
+        }
+        // dx = -1
+        // dy = -1
+        break;
+    case (corner::cTopRight):
+        if (eta.sum != 0) {
+            eta.x = static_cast<double>(
+                        cl.data[cluster_center_index + 1] +
+                        cl.data[cluster_center_index - ClusterSizeX + 1]) /
+                    static_cast<double>(eta.sum);
+            eta.y = static_cast<double>(cl.data[cluster_center_index] +
+                                        cl.data[cluster_center_index + 1]) /
+                    static_cast<double>(eta.sum);
+        }
+        // dx = 0
+        // dy = -1
+        break;
+    case (corner::cBottomLeft):
+        if (eta.sum != 0) {
+            eta.x = static_cast<double>(
+                        cl.data[cluster_center_index] +
+                        cl.data[cluster_center_index + ClusterSizeX]) /
+                    static_cast<double>(eta.sum);
+            eta.y = static_cast<double>(
+                        cl.data[cluster_center_index + ClusterSizeX] +
+                        cl.data[cluster_center_index + ClusterSizeX - 1]) /
+                    static_cast<double>(eta.sum);
+        }
+        // dx = -1
+        // dy = 0
+        break;
+    case (corner::cBottomRight):
+        if (eta.sum != 0) {
+            eta.x = static_cast<double>(
+                        cl.data[cluster_center_index + 1] +
+                        cl.data[cluster_center_index + ClusterSizeX + 1]) /
+                    static_cast<double>(eta.sum);
+            eta.y = static_cast<double>(
+                        cl.data[cluster_center_index + ClusterSizeX] +
+                        cl.data[cluster_center_index + ClusterSizeX + 1]) /
+                    static_cast<double>(eta.sum);
+        }
+        // dx = 0
+        // dy = 0
+        break;
+    }
+
+    eta.c = c;
+
+    return eta;
+}
+
+template <typename T>
+Eta2<T> calculate_eta2(const Cluster<T, 2, 2, uint16_t> &cl) {
+    Eta2<T> eta{};
+
+    // TODO: maybe have as member function of cluster
+    const uint8_t photon_hit_index =
+        std::max_element(cl.data.begin(), cl.data.end()) - cl.data.begin();
+
+    eta.c = static_cast<corner>(3 - photon_hit_index);
+
+    switch (eta.c) {
+    case corner::cTopLeft:
+        calculate_eta2(eta, cl.data[2], cl.data[3], cl.data[1], cl.data[3]);
+        break;
+    case corner::cTopRight:
+        calculate_eta2(eta, cl.data[2], cl.data[3], cl.data[0], cl.data[2]);
+        break;
+    case corner::cBottomLeft:
+        calculate_eta2(eta, cl.data[0], cl.data[1], cl.data[1], cl.data[3]);
+        break;
+    case corner::cBottomRight:
+        calculate_eta2(eta, cl.data[0], cl.data[1], cl.data[0], cl.data[2]);
+        break;
+    }
+
     eta.sum = cl.sum();
 
     return eta;
 }
 
+template <typename T>
+Eta2<T> calculate_full_eta2(const Cluster<T, 2, 2, uint16_t> &cl) {
+
+    Eta2<T> eta{};
+
+    eta.sum = cl.sum();
+
+    const uint8_t photon_hit_index =
+        std::max_element(cl.data.begin(), cl.data.end()) - cl.data.begin();
+
+    eta.c = static_cast<corner>(3 - photon_hit_index);
+
+    if (eta.sum != 0) {
+        eta.x = static_cast<double>(cl.data[1] + cl.data[3]) /
+                static_cast<double>(eta.sum);
+        eta.y = static_cast<double>(cl.data[2] + cl.data[3]) /
+                static_cast<double>(eta.sum);
+    }
+
+    return eta;
+}
+
 // TODO generalize
 template <typename T>
-Eta2<T> calculate_eta2(const Cluster<T, 1, 2, int16_t> &cl) {
+Eta2<T> calculate_eta2(const Cluster<T, 1, 2, uint16_t> &cl) {
     Eta2<T> eta{};
 
     eta.x = 0;
-    eta.y = static_cast<double>(cl.data[0]) / cl.data[1];
+    eta.y = static_cast<double>(cl.data[1]) / cl.data[0];
     eta.sum = cl.sum();
 }
 
 template <typename T>
-Eta2<T> calculate_eta2(const Cluster<T, 2, 1, int16_t> &cl) {
+Eta2<T> calculate_eta2(const Cluster<T, 2, 1, uint16_t> &cl) {
     Eta2<T> eta{};
 
-    eta.x = static_cast<double>(cl.data[0]) / cl.data[1];
+    eta.x = static_cast<double>(cl.data[1]) / cl.data[0];
     eta.y = 0;
     eta.sum = cl.sum();
 }
 
-// calculates Eta3 for 3x3 cluster based on code from analyze_cluster
-// TODO only supported for 3x3 Clusters
-template <typename T> Eta2<T> calculate_eta3(const Cluster<T, 3, 3> &cl) {
+/**
+ * @brief calculates cross Eta3 for 3x3 cluster
+ * cross Eta3 calculates the eta by taking into account only the cross pixels
+ * {top, bottom, left, right, center}
+ */
+template <typename T, typename CoordType = uint16_t>
+Eta2<T> calculate_cross_eta3(const Cluster<T, 3, 3, CoordType> &cl) {
 
     Eta2<T> eta{};
 
-    T sum = 0;
+    T photon_energy = cl.sum();
 
-    std::for_each(std::begin(cl.data), std::end(cl.data),
-                  [&sum](T x) { sum += x; });
-
-    eta.sum = sum;
+    eta.sum = photon_energy;
 
     if ((cl.data[3] + cl.data[4] + cl.data[5]) != 0)
 
-        eta.x = static_cast<double>(-cl.data[3] + cl.data[3 + 2]) /
+        eta.x =
+            static_cast<double>(-cl.data[3] + cl.data[3 + 2]) /
 
-                (cl.data[3] + cl.data[4] + cl.data[5]); // (-1,1)
+            static_cast<double>(cl.data[3] + cl.data[4] + cl.data[5]); // (-1,1)
 
     if ((cl.data[1] + cl.data[4] + cl.data[7]) != 0)
 
         eta.y = static_cast<double>(-cl.data[1] + cl.data[2 * 3 + 1]) /
 
-                (cl.data[1] + cl.data[4] + cl.data[7]);
+                static_cast<double>(cl.data[1] + cl.data[4] + cl.data[7]);
 
     return eta;
 }
 
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = uint16_t>
+Eta2<T> calculate_cross_eta3(
+    const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
+
+    static_assert(ClusterSizeX > 2 && ClusterSizeY > 2,
+                  "calculate_eta3 only defined for clusters larger than 2x2");
+
+    if constexpr (ClusterSizeX != 3 || ClusterSizeY != 3) {
+        auto reduced_cluster = reduce_cluster_to_3x3(cl);
+        return calculate_cross_eta3(reduced_cluster);
+    } else {
+        return calculate_cross_eta3(cl);
+    }
+}
+
+/**
+ * @brief calculates Eta3 for 3x3 cluster
+ * It calculates the eta by taking into account all pixels in the 3x3 cluster
+ */
+template <typename T, typename CoordType = uint16_t>
+Eta2<T> calculate_eta3(const Cluster<T, 3, 3, CoordType> &cl) {
+
+    Eta2<T> eta{};
+
+    T photon_energy = cl.sum();
+
+    eta.sum = photon_energy;
+
+    // TODO: how do we handle potential arithmetic overflows? - T could be
+    // uint16
+    if (photon_energy != 0) {
+        std::array<T, 2> column_sums{
+            static_cast<T>(cl.data[0] + cl.data[3] + cl.data[6]),
+            static_cast<T>(cl.data[2] + cl.data[5] + cl.data[8])};
+
+        eta.x = static_cast<double>(-column_sums[0] + column_sums[1]) /
+                static_cast<double>(photon_energy);
+
+        std::array<T, 2> row_sums{
+            static_cast<T>(cl.data[0] + cl.data[1] + cl.data[2]),
+            static_cast<T>(cl.data[6] + cl.data[7] + cl.data[8])};
+
+        eta.y = static_cast<double>(-row_sums[0] + row_sums[1]) /
+                static_cast<double>(photon_energy);
+    }
+
+    return eta;
+}
+
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = uint16_t>
+Eta2<T>
+calculate_eta3(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
+
+    static_assert(ClusterSizeX > 2 && ClusterSizeY > 2,
+                  "calculate_eta3 only defined for clusters larger than 2x2");
+
+    if constexpr (ClusterSizeX != 3 || ClusterSizeY != 3) {
+        auto reduced_cluster = reduce_cluster_to_3x3(cl);
+        return calculate_eta3(reduced_cluster);
+    } else {
+        return calculate_eta3(cl);
+    }
+}
+
 } // namespace aare

include/aare/CircularFifo.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 
 #include <chrono>

include/aare/Cluster.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 
 /************************************************
  * @file Cluster.hpp
@@ -8,6 +9,7 @@
 
 #pragma once
 
+#include "defs.hpp"
 #include <algorithm>
 #include <array>
 #include <cstdint>
@@ -17,6 +19,10 @@
 namespace aare {
 
 // requires clause c++20 maybe update
+
+/**
+ * @brief Cluster struct
+ */
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
           typename CoordType = uint16_t>
 struct Cluster {
@@ -27,8 +33,11 @@ struct Cluster {
     static_assert(ClusterSizeX > 0 && ClusterSizeY > 0,
                   "Cluster sizes must be bigger than zero");
 
+    /// @brief Cluster center x coordinate (in pixel coordinates)
     CoordType x;
+    /// @brief Cluster center y coordinate (in pixel coordinates)
     CoordType y;
+    /// @brief Cluster data stored in row-major order starting from top-left
     std::array<T, ClusterSizeX * ClusterSizeY> data;
 
     static constexpr uint8_t cluster_size_x = ClusterSizeX;
@@ -36,16 +45,18 @@ struct Cluster {
     using value_type = T;
     using coord_type = CoordType;
 
+    /**
+     * @brief Sum of all elements in the cluster
+     */
     T sum() const { return std::accumulate(data.begin(), data.end(), T{}); }
 
     // TODO: handle 1 dimensional clusters
-    // TODO: change int to corner
     /**
      * @brief sum of 2x2 subcluster with highest energy
      * @return photon energy of subcluster, 2x2 subcluster index relative to
      * cluster center
      */
-    std::pair<T, int> max_sum_2x2() const {
+    Sum_index_pair<T, corner> max_sum_2x2() const {
 
         if constexpr (cluster_size_x == 3 && cluster_size_y == 3) {
             std::array<T, 4> sum_2x2_subclusters;
@@ -56,9 +67,11 @@ struct Cluster {
             int index = std::max_element(sum_2x2_subclusters.begin(),
                                          sum_2x2_subclusters.end()) -
                         sum_2x2_subclusters.begin();
-            return std::make_pair(sum_2x2_subclusters[index], index);
+            return Sum_index_pair<T, corner>{sum_2x2_subclusters[index],
+                                             corner{index}};
         } else if constexpr (cluster_size_x == 2 && cluster_size_y == 2) {
-            return std::make_pair(data[0] + data[1] + data[2] + data[3], 0);
+            return Sum_index_pair<T, corner>{
+                data[0] + data[1] + data[2] + data[3], corner{0}};
         } else {
             constexpr size_t cluster_center_index =
                 (ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
@@ -97,7 +110,8 @@ struct Cluster {
             int index = std::max_element(sum_2x2_subcluster.begin(),
                                          sum_2x2_subcluster.end()) -
                         sum_2x2_subcluster.begin();
-            return std::make_pair(sum_2x2_subcluster[index], index);
+            return Sum_index_pair<T, corner>{sum_2x2_subcluster[index],
+                                             corner{index}};
         }
     }
 };
@@ -107,66 +121,71 @@ struct Cluster {
  * highest sum.
  * @param c Cluster to reduce
  * @return reduced cluster
+ * @note The cluster is filled using row major ordering starting at the top-left
+ * (thus for a max subcluster in the top left cornern the photon hit is at
+ * the fourth position)
  */
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType = int16_t>
+          typename CoordType = uint16_t>
 Cluster<T, 2, 2, CoordType>
 reduce_to_2x2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &c) {
 
     static_assert(ClusterSizeX >= 2 && ClusterSizeY >= 2,
                   "Cluster sizes must be at least 2x2 for reduction to 2x2");
 
-    // TODO maybe add sanity check and check that center is in max subcluster
-    Cluster<T, 2, 2, CoordType> result;
+    Cluster<T, 2, 2, CoordType> result{};
 
     auto [sum, index] = c.max_sum_2x2();
 
-    int16_t cluster_center_index =
+    constexpr int16_t cluster_center_index =
         (ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
 
-    int16_t index_bottom_left_max_2x2_subcluster =
-        (int(index / (ClusterSizeX - 1))) * ClusterSizeX +
-        index % (ClusterSizeX - 1);
+    int16_t index_top_left_max_2x2_subcluster = cluster_center_index;
+    switch (index) {
+    case corner::cTopLeft:
+        index_top_left_max_2x2_subcluster -= (ClusterSizeX + 1);
+        break;
+    case corner::cTopRight:
+        index_top_left_max_2x2_subcluster -= ClusterSizeX;
+        break;
+    case corner::cBottomLeft:
+        index_top_left_max_2x2_subcluster -= 1;
+        break;
+    case corner::cBottomRight:
+        // no change needed
+        break;
+    }
 
-    result.x =
-        c.x + (index_bottom_left_max_2x2_subcluster - cluster_center_index) %
-                  ClusterSizeX;
+    result.x = c.x;
+    result.y = c.y;
 
-    result.y =
-        c.y - (index_bottom_left_max_2x2_subcluster - cluster_center_index) /
-                  ClusterSizeX;
     result.data = {
-        c.data[index_bottom_left_max_2x2_subcluster],
-        c.data[index_bottom_left_max_2x2_subcluster + 1],
-        c.data[index_bottom_left_max_2x2_subcluster + ClusterSizeX],
-        c.data[index_bottom_left_max_2x2_subcluster + ClusterSizeX + 1]};
+        c.data[index_top_left_max_2x2_subcluster],
+        c.data[index_top_left_max_2x2_subcluster + 1],
+        c.data[index_top_left_max_2x2_subcluster + ClusterSizeX],
+        c.data[index_top_left_max_2x2_subcluster + ClusterSizeX + 1]};
+
     return result;
 }
 
 template <typename T>
-Cluster<T, 2, 2, int16_t> reduce_to_2x2(const Cluster<T, 3, 3, int16_t> &c) {
-    Cluster<T, 2, 2, int16_t> result;
+Cluster<T, 2, 2, uint16_t> reduce_to_2x2(const Cluster<T, 3, 3, uint16_t> &c) {
+    Cluster<T, 2, 2, uint16_t> result{};
 
     auto [s, i] = c.max_sum_2x2();
+    result.x = c.x;
+    result.y = c.y;
     switch (i) {
-    case 0:
-        result.x = c.x - 1;
-        result.y = c.y + 1;
+    case corner::cTopLeft:
         result.data = {c.data[0], c.data[1], c.data[3], c.data[4]};
         break;
-    case 1:
-        result.x = c.x;
-        result.y = c.y + 1;
+    case corner::cTopRight:
         result.data = {c.data[1], c.data[2], c.data[4], c.data[5]};
         break;
-    case 2:
-        result.x = c.x - 1;
-        result.y = c.y;
+    case corner::cBottomLeft:
         result.data = {c.data[3], c.data[4], c.data[6], c.data[7]};
         break;
-    case 3:
-        result.x = c.x;
-        result.y = c.y;
+    case corner::cBottomRight:
         result.data = {c.data[4], c.data[5], c.data[7], c.data[8]};
         break;
     }
@@ -174,43 +193,8 @@ Cluster<T, 2, 2, int16_t> reduce_to_2x2(const Cluster<T, 3, 3, int16_t> &c) {
     return result;
 }
 
-template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType = int16_t>
-inline std::pair<T, uint16_t>
-max_3x3_sum(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cluster) {
-
-    if constexpr (ClusterSizeX == 3 && ClusterSizeY == 3) {
-        return std::make_pair(cluster.sum(), 0);
-    } else {
-
-        size_t index = 0;
-        T max_3x3_subcluster_sum = 0;
-        for (size_t i = 0; i < ClusterSizeY - 2; ++i) {
-            for (size_t j = 0; j < ClusterSizeX - 2; ++j) {
-
-                T sum = cluster.data[i * ClusterSizeX + j] +
-                        cluster.data[i * ClusterSizeX + j + 1] +
-                        cluster.data[i * ClusterSizeX + j + 2] +
-                        cluster.data[(i + 1) * ClusterSizeX + j] +
-                        cluster.data[(i + 1) * ClusterSizeX + j + 1] +
-                        cluster.data[(i + 1) * ClusterSizeX + j + 2] +
-                        cluster.data[(i + 2) * ClusterSizeX + j] +
-                        cluster.data[(i + 2) * ClusterSizeX + j + 1] +
-                        cluster.data[(i + 2) * ClusterSizeX + j + 2];
-                if (sum > max_3x3_subcluster_sum) {
-                    max_3x3_subcluster_sum = sum;
-                    index = i * (ClusterSizeX - 2) + j;
-                }
-            }
-        }
-
-        return std::make_pair(max_3x3_subcluster_sum, index);
-    }
-}
-
 /**
- * @brief Reduce a cluster to a 3x3 cluster by selecting the 3x3 block with the
- * highest sum.
+ * @brief Reduce a cluster to a 3x3 cluster
  * @param c Cluster to reduce
  * @return reduced cluster
  */
@@ -222,40 +206,24 @@ reduce_to_3x3(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &c) {
     static_assert(ClusterSizeX >= 3 && ClusterSizeY >= 3,
                   "Cluster sizes must be at least 3x3 for reduction to 3x3");
 
-    Cluster<T, 3, 3, CoordType> result;
-
-    // TODO maybe add sanity check and check that center is in max subcluster
-
-    auto [sum, index] = max_3x3_sum(c);
+    Cluster<T, 3, 3, CoordType> result{};
 
     int16_t cluster_center_index =
         (ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
 
-    int16_t index_center_max_3x3_subcluster =
-        (int(index / (ClusterSizeX - 2))) * ClusterSizeX + ClusterSizeX +
-        index % (ClusterSizeX - 2) + 1;
+    result.x = c.x;
+    result.y = c.y;
 
-    int16_t index_3x3_subcluster_cluster_center =
-        int((cluster_center_index - 1 - ClusterSizeX) / ClusterSizeX) *
-            (ClusterSizeX - 2) +
-        (cluster_center_index - 1 - ClusterSizeX) % ClusterSizeX;
+    result.data = {c.data[cluster_center_index - ClusterSizeX - 1],
+                   c.data[cluster_center_index - ClusterSizeX],
+                   c.data[cluster_center_index - ClusterSizeX + 1],
+                   c.data[cluster_center_index - 1],
+                   c.data[cluster_center_index],
+                   c.data[cluster_center_index + 1],
+                   c.data[cluster_center_index + ClusterSizeX - 1],
+                   c.data[cluster_center_index + ClusterSizeX],
+                   c.data[cluster_center_index + ClusterSizeX + 1]};
 
-    result.x =
-        c.x + (index % (ClusterSizeX - 2) -
-               (index_3x3_subcluster_cluster_center % (ClusterSizeX - 2)));
-    result.y =
-        c.y - (index / (ClusterSizeX - 2) -
-               (index_3x3_subcluster_cluster_center / (ClusterSizeX - 2)));
-
-    result.data = {c.data[index_center_max_3x3_subcluster - ClusterSizeX - 1],
-                   c.data[index_center_max_3x3_subcluster - ClusterSizeX],
-                   c.data[index_center_max_3x3_subcluster - ClusterSizeX + 1],
-                   c.data[index_center_max_3x3_subcluster - 1],
-                   c.data[index_center_max_3x3_subcluster],
-                   c.data[index_center_max_3x3_subcluster + 1],
-                   c.data[index_center_max_3x3_subcluster + ClusterSizeX - 1],
-                   c.data[index_center_max_3x3_subcluster + ClusterSizeX],
-                   c.data[index_center_max_3x3_subcluster + ClusterSizeX + 1]};
     return result;
 }
 

include/aare/ClusterCollector.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include <atomic>
 #include <thread>
@@ -5,6 +6,7 @@
 #include "aare/ClusterFinderMT.hpp"
 #include "aare/ClusterVector.hpp"
 #include "aare/ProducerConsumerQueue.hpp"
+#include "aare/defs.hpp"
 
 namespace aare {
 

include/aare/ClusterFile.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 
 #include "aare/Cluster.hpp"
@@ -189,6 +190,16 @@ class ClusterFile {
         }
     }
 
+    /**
+     * @brief Return the current position in the file (bytes)
+     */
+    int64_t tell() {
+        if (!fp) {
+            throw std::runtime_error(LOCATION + "File not opened");
+        }
+        return ftell(fp);
+    }
+
     /** @brief Open the file in specific mode
      *
      */
@@ -354,15 +365,20 @@ template <typename ClusterType, typename Enable>
 ClusterVector<ClusterType>
 ClusterFile<ClusterType, Enable>::read_frame_without_cut() {
     if (m_mode != "r") {
-        throw std::runtime_error("File not opened for reading");
+        throw std::runtime_error(LOCATION + "File not opened for reading");
     }
     if (m_num_left) {
         throw std::runtime_error(
-            "There are still photons left in the last frame");
+            LOCATION + "There are still photons left in the last frame");
     }
     int32_t frame_number;
     if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) {
-        throw std::runtime_error(LOCATION + "Could not read frame number");
+        if (feof(fp))
+            throw std::runtime_error(LOCATION + "Unexpected end of file");
+        else if (ferror(fp))
+            throw std::runtime_error(LOCATION + "Error reading from file");
+
+        throw std::runtime_error(LOCATION + "Unexpected error (not feof or ferror) when reading frame number");
     }
 
     int32_t n_clusters; // Saved as 32bit integer in the cluster file
@@ -438,8 +454,8 @@ bool ClusterFile<ClusterType, Enable>::is_selected(ClusterType &cl) {
 
     if (m_noise_map) {
         auto sum_1x1 = cl.data[cluster_center_index]; // central pixel
-        auto sum_2x2 = cl.max_sum_2x2().first; // highest sum of 2x2 subclusters
-        auto total_sum = cl.sum();             // sum of all pixels
+        auto sum_2x2 = cl.max_sum_2x2().sum; // highest sum of 2x2 subclusters
+        auto total_sum = cl.sum();           // sum of all pixels
 
         auto noise =
             (*m_noise_map)(cl.y, cl.x); // TODO! check if this is correct

include/aare/ClusterFileSink.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include <atomic>
 #include <filesystem>
@@ -10,7 +11,8 @@
 namespace aare {
 
 template <typename ClusterType,
-          typename = std::enable_if_t<is_cluster_v<ClusterType>>>
+          typename = std::enable_if_t<is_cluster_v<ClusterType>>,
+          typename = std::enable_if_t<no_2x2_cluster<ClusterType>::value>>
 class ClusterFileSink {
     ProducerConsumerQueue<ClusterVector<ClusterType>> *m_source;
     std::atomic<bool> m_stop_requested{false};

include/aare/ClusterFinder.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include "aare/ClusterFile.hpp"
 #include "aare/ClusterVector.hpp"
@@ -10,8 +11,16 @@
 
 namespace aare {
 
+template <typename ClusterType,
+          typename = std::enable_if_t<is_cluster_v<ClusterType>>>
+struct no_2x2_cluster {
+    constexpr static bool value =
+        ClusterType::cluster_size_x > 2 && ClusterType::cluster_size_y > 2;
+};
+
 template <typename ClusterType = Cluster<int32_t, 3, 3>,
-          typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double>
+          typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double,
+          typename = std::enable_if_t<no_2x2_cluster<ClusterType>::value>>
 class ClusterFinder {
     Shape<2> m_image_size;
     const PEDESTAL_TYPE m_nSigma;

include/aare/ClusterFinderMT.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include <atomic>
 #include <cstdint>
@@ -32,7 +33,8 @@ struct FrameWrapper {
  * @tparam CT type of the cluster data
  */
 template <typename ClusterType = Cluster<int32_t, 3, 3>,
-          typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double>
+          typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double,
+          typename = std::enable_if_t<no_2x2_cluster<ClusterType>::value>>
 class ClusterFinderMT {
 
   protected:

include/aare/ClusterVector.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include "aare/Cluster.hpp" //TODO maybe store in seperate file !!!
 #include <algorithm>
@@ -28,7 +29,7 @@ class ClusterVector; // Forward declaration
  * needed.
  * @tparam T data type of the pixels in the cluster
  * @tparam CoordType data type of the x and y coordinates of the cluster
- * (normally int16_t)
+ * (normally uint16_t)
  */
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
           typename CoordType>
@@ -86,15 +87,14 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
     /**
      * @brief Sum the pixels in the 2x2 subcluster with the biggest pixel sum in
      * each cluster
-     * @return std::vector<T> vector of sums for each cluster
+     * @return vector of sums index pairs for each cluster
      */
-    std::vector<T> sum_2x2() {
-        std::vector<T> sums_2x2(m_data.size());
+    std::vector<Sum_index_pair<T, corner>> sum_2x2() {
+        std::vector<Sum_index_pair<T, corner>> sums_2x2(m_data.size());
 
-        std::transform(m_data.begin(), m_data.end(), sums_2x2.begin(),
-                       [](const ClusterType &cluster) {
-                           return cluster.max_sum_2x2().first;
-                       });
+        std::transform(
+            m_data.begin(), m_data.end(), sums_2x2.begin(),
+            [](const ClusterType &cluster) { return cluster.max_sum_2x2(); });
 
         return sums_2x2;
     }
@@ -177,9 +177,12 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
  * highest sum.
  * @param cv Clustervector containing clusters to reduce
  * @return Clustervector with reduced clusters
+ * @note The cluster is filled using row major ordering starting at the top-left
+ * (thus for a max subcluster in the top left cornern the photon hit is at
+ * the fourth position)
  */
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType = uint16_t>
+          typename CoordType>
 ClusterVector<Cluster<T, 2, 2, CoordType>> reduce_to_2x2(
     const ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>
         &cv) {
@@ -191,13 +194,12 @@ ClusterVector<Cluster<T, 2, 2, CoordType>> reduce_to_2x2(
 }
 
 /**
- * @brief Reduce a cluster to a 3x3 cluster by selecting the 3x3 block with the
- * highest sum.
+ * @brief Reduce a cluster to a 3x3 cluster
  * @param cv Clustervector containing clusters to reduce
  * @return Clustervector with reduced clusters
  */
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType = uint16_t>
+          typename CoordType>
 ClusterVector<Cluster<T, 3, 3, CoordType>> reduce_to_3x3(
     const ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>
         &cv) {

include/aare/CtbRawFile.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 
 #include "aare/FileInterface.hpp"

include/aare/DetectorGeometry.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include "aare/RawMasterFile.hpp" //ROI refactor away
 #include "aare/defs.hpp"

include/aare/Dtype.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include <cstdint>
 #include <map>

include/aare/File.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include "aare/FileInterface.hpp"
 #include <memory>

include/aare/FileInterface.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include "aare/Dtype.hpp"
 #include "aare/Frame.hpp"
@@ -33,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) + " }";
+    // }
 };
 
 /**

include/aare/FilePtr.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include <cstdio>
 #include <filesystem>

include/aare/Fit.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 
 #include <cmath>

include/aare/Frame.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include "aare/Dtype.hpp"
 #include "aare/NDArray.hpp"

include/aare/GainMap.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 /************************************************
  * @file GainMap.hpp
  * @short function to apply gain map of image size to a vector of clusters -

include/aare/Interpolator.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 
 #include "aare/CalculateEta.hpp"
@@ -16,8 +17,27 @@ 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])
+     */
     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])
+     */
     NDArray<double, 3> m_ietay;
 
     NDArray<double, 1> m_etabinsx;
@@ -25,105 +45,207 @@ class Interpolator {
     NDArray<double, 1> m_energy_bins;
 
   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 xbins bin edges for etaX
+     * @param ybins bin edges for etaY
+     * @param ebins bin edges for photon energy
+     */
     Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
                  NDView<double, 1> ybins, NDView<double, 1> ebins);
+
+    /**
+     * @brief Constructor for the Interpolator class
+     * @param xbins bin edges for etaX
+     * @param ybins bin edges for etaY
+     * @param ebins bin edges for photon energy
+     */
+    Interpolator(NDView<double, 1> xbins, NDView<double, 1> ybins,
+                 NDView<double, 1> ebins);
+
+    /**
+     * @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)
+     */
+    void rosenblatttransform(NDView<double, 3> etacube);
+
     NDArray<double, 3> get_ietax() { return m_ietax; }
     NDArray<double, 3> get_ietay() { return m_ietay; }
 
-    template <typename ClusterType,
+    /**
+     * @brief interpolates the cluster centers for all clusters to a better
+     * precision
+     * @tparam ClusterType Type of Clusters to interpolate
+     * @tparam Etafunction Function object that calculates desired eta default:
+     * 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))
+     */
+    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);
+    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;
+
+  private:
+    /**
+     * @brief bilinear interpolation of the transformed eta values
+     * @param ix index of etaX bin
+     * @param iy index of etaY bin
+     * @param ie index of energy bin
+     * @return pair of interpolated transformed eta values (ietax, ietay)
+     */
+    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;
 };
 
-// TODO: generalize to support any clustertype!!! otherwise add std::enable_if_t
-// to only take Cluster2x2 and Cluster3x3
-template <typename ClusterType, typename Enable>
+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 {
+    auto next_index_y = static_cast<ssize_t>(iy + 1) >= m_ietax.shape(1)
+                            ? m_ietax.shape(1) - 1
+                            : iy + 1;
+    auto next_index_x = static_cast<ssize_t>(ix + 1) >= m_ietax.shape(0)
+                            ? m_ietax.shape(0) - 1
+                            : ix + 1;
+
+    // bilinear interpolation
+    double ietax_interp_left = linear_interpolation(
+        {m_etabinsy(iy), m_etabinsy(iy + 1)},
+        {m_ietax(ix, iy, ie), m_ietax(ix, next_index_y, ie)}, eta.y);
+    double ietax_interp_right =
+        linear_interpolation({m_etabinsy(iy), m_etabinsy(iy + 1)},
+                             {m_ietax(next_index_x, iy, ie),
+                              m_ietax(next_index_x, next_index_y, ie)},
+                             eta.y);
+
+    // transformed photon position x between [0,1]
+    double ietax_interpolated =
+        linear_interpolation({m_etabinsx(ix), m_etabinsx(ix + 1)},
+                             {ietax_interp_left, ietax_interp_right}, eta.x);
+
+    double ietay_interp_left = linear_interpolation(
+        {m_etabinsx(ix), m_etabinsx(ix + 1)},
+        {m_ietay(ix, iy, ie), m_ietay(next_index_x, iy, ie)}, eta.x);
+    double ietay_interp_right =
+        linear_interpolation({m_etabinsx(ix), m_etabinsx(ix + 1)},
+                             {m_ietay(ix, next_index_y, ie),
+                              m_ietay(next_index_x, next_index_y, ie)},
+                             eta.x);
+
+    // transformed photon position y between [0,1]
+    double ietay_interpolated =
+        linear_interpolation({m_etabinsy(iy), m_etabinsy(iy + 1)},
+                             {ietay_interp_left, ietay_interp_right}, eta.y);
+
+    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) {
+Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) const {
     std::vector<Photon> photons;
     photons.reserve(clusters.size());
 
-    if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) {
-        for (const ClusterType &cluster : clusters) {
+    for (const ClusterType &cluster : clusters) {
 
-            auto eta = calculate_eta2(cluster);
+        auto eta = EtaFunction(cluster);
 
-            Photon photon;
-            photon.x = cluster.x;
-            photon.y = cluster.y;
-            photon.energy = static_cast<decltype(photon.energy)>(eta.sum);
+        Photon photon;
+        photon.x = cluster.x;
+        photon.y = cluster.y;
+        photon.energy = static_cast<decltype(photon.energy)>(eta.sum);
 
-            // auto ie = nearest_index(m_energy_bins, photon.energy)-1;
-            // auto ix = nearest_index(m_etabinsx, eta.x)-1;
-            // auto iy = nearest_index(m_etabinsy, eta.y)-1;
-            // 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);
+        auto uniform_coordinates = transform_eta_values(eta);
 
-            // fmt::print("ex: {}, ix: {}, iy: {}\n", ie, ix, iy);
+        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,
+                                     ClusterType::cluster_size_x,
+                                     ClusterType::cluster_size_y,
+                                     typename ClusterType::coord_type>) {
+            double dX{}, dY{};
 
-            double dX, dY;
-            // cBottomLeft = 0,
-            // cBottomRight = 1,
-            // cTopLeft = 2,
-            // cTopRight = 3
             // TODO: could also chaneg the sign of the eta calculation
-            switch (static_cast<corner>(eta.c)) {
+            switch (eta.c) {
             case corner::cTopLeft:
-                dX = 0.0;
-                dY = 0.0;
+                dX = -1.0;
+                dY = -1.0;
                 break;
             case corner::cTopRight:;
-                dX = 1.0;
-                dY = 0.0;
+                dX = 0.0;
+                dY = -1.0;
                 break;
             case corner::cBottomLeft:
-                dX = 0.0;
-                dY = 1.0;
+                dX = -1.0;
+                dY = 0.0;
                 break;
             case corner::cBottomRight:
-                dX = 1.0;
-                dY = 1.0;
+                dX = 0.0;
+                dY = 0.0;
                 break;
             }
-            photon.x -= m_ietax(ix, iy, ie) - dX;
-            photon.y -= m_ietay(ix, iy, ie) - dY;
-            photons.push_back(photon);
-        }
-    } else if (clusters.cluster_size_x() == 2 ||
-               clusters.cluster_size_y() == 2) {
-        for (const ClusterType &cluster : clusters) {
-            auto eta = calculate_eta2(cluster);
-
-            Photon photon;
-            photon.x = cluster.x;
-            photon.y = cluster.y;
-            photon.energy = static_cast<decltype(photon.energy)>(eta.sum);
-
-            // Now do some actual interpolation.
-            // Find which energy bin the cluster is in
-            //  auto ie = nearest_index(m_energy_bins, photon.energy)-1;
-            //  auto ix = nearest_index(m_etabinsx, eta.x)-1;
-            //  auto iy = nearest_index(m_etabinsy, eta.y)-1;
-            // 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);
-
-            // TODO: why 2?
-            photon.x -=
-                m_ietax(ix, iy, ie); // eta goes between 0 and 1 but we could
-                                     // move the hit anywhere in the 2x2
-            photon.y -= m_ietay(ix, iy, ie);
-            photons.push_back(photon);
+            photon.x = photon.x + 0.5 + uniform_coordinates.x +
+                       dX; // use pixel center + 0.5
+            photon.y =
+                photon.y + 0.5 + uniform_coordinates.y +
+                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;
         }
 
-    } else {
-        throw std::runtime_error(
-            "Only 3x3 and 2x2 clusters are supported for interpolation");
+        photons.push_back(photon);
     }
 
     return photons;

include/aare/JungfrauDataFile.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include <cstdint>
 #include <filesystem>

include/aare/NDArray.hpp

@@ -1,12 +1,10 @@
+// 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"
 
@@ -22,15 +20,21 @@ TODO! Add expression templates for operators
 namespace aare {
 
 template <typename T, ssize_t Ndim = 2>
-class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
+class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim>,
+                public NDIndexOps<NDArray<T, Ndim>, T, 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. Will construct an empty NDArray.
+     * @brief Default constructor. Constructs an empty NDArray.
      *
      */
     NDArray() : shape_(), strides_(c_strides<Ndim>(shape_)), data_(nullptr) {};
@@ -43,8 +47,7 @@ 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.
@@ -56,6 +59,10 @@ 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.
@@ -66,44 +73,67 @@ 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());
     }
 
-    // Move constructor
+    /**
+     * @brief Move construct a new NDArray object. Cheap since it just
+     * reassigns the pointer and copy size/strides.
+     *
+     * @param other
+     */
     NDArray(NDArray &&other) noexcept
         : shape_(other.shape_), strides_(c_strides<Ndim>(shape_)),
           size_(other.size_), data_(other.data_) {
-        other.reset(); // TODO! is this necessary?
+        other.reset(); // Needed to avoid double free
     }
 
-
-    //Move constructor from an an array with Ndim + 1
+    /**
+     * @brief Move construct a new NDArray object from an array with Ndim + 1.
+     * Can be used to drop a dimension cheaply.
+     * @param other
+     */
     template <ssize_t M, typename = std::enable_if_t<(M == Ndim + 1)>>
-    NDArray(NDArray<T, M> &&other) 
+    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 +
-                                         "Size mismatch in move constructor of NDArray<T, Ndim-1>");
-            }
+        // 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();
     }
 
-    // Copy constructor
+    /**
+     * @brief Copy construct a new NDArray object from another NDArray.
+     *
+     * @param other
+     */
     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_);
     }
 
-    // Conversion operator from array expression to array
+    /**
+     * @brief Conversion from a ArrayExpr to an actual NDArray. Used when
+     * the expression is evaluated and data needed.
+     *
+     * @tparam E
+     * @param expr
+     */
     template <typename E>
     NDArray(ArrayExpr<E, Ndim> &&expr) : NDArray(expr.shape()) {
         for (size_t i = 0; i < size_; ++i) {
@@ -111,23 +141,83 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
         }
     }
 
+    /**
+     * @brief Destroy the NDArray object. Frees the allocated memory.
+     *
+     */
     ~NDArray() { delete[] data_; }
 
-    auto begin() { return data_; }
-    auto end() { return data_ + size_; }
+    ///////////////////////////////////////////////////////////////////////////////
+    // Iterators and indexing
+    //
+    ///////////////////////////////////////////////////////////////////////////////
 
-    auto begin() const { return data_; }
-    auto end() const { return data_ + size_; }
+    using NDIndexOps<NDArray<T, Ndim>, T, Ndim>::operator();
+    using NDIndexOps<NDArray<T, Ndim>, T, Ndim>::operator[];
+
+    auto *begin() { return data_; }
+    const auto *begin() const { return data_; }
+
+    auto *end() { return data_ + size_; }
+    const auto *end() const { return data_ + size_; }
+
+    /* @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;
+    }
 
     using value_type = T;
 
-    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);
+    ///////////////////////////////////////////////////////////////////////////////
+    // Assignments
+    //
+    ///////////////////////////////////////////////////////////////////////////////
 
-    // Write directly to the data array, or create a new one
+    /**
+     * @brief Copy to the NDArray from an std::array. If the size of the array
+     * is different we reallocate the data.
+     *
+     */
     template <size_t Size>
     NDArray<T, 1> &operator=(const std::array<T, Size> &other) {
         if (Size != size_) {
@@ -141,12 +231,94 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
         return *this;
     }
 
-    // NDArray& operator/=(const NDArray& other);
+    /**
+     * @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;
+    }
 
+    /**
+     * @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 (uint32_t i = 0; i < size_; ++i) {
+            for (size_t i = 0; i < size_; ++i) {
                 data_[i] /= other(i);
             }
             return *this;
@@ -154,67 +326,139 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
         throw(std::runtime_error("Shape of NDArray must match"));
     }
 
-    NDArray<bool, Ndim> operator>(const NDArray &other);
+    /**
+     * @brief Assign a scalar value to all elements in the NDArray.
+     */
+    NDArray &operator=(const T &value) {
+        std::fill_n(data_, size_, value);
+        return *this;
+    }
 
-    bool operator==(const NDArray &other) const;
-    bool operator!=(const NDArray &other) const;
+    /**
+     * @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;
+    }
 
-    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 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 & /*mask*/);
+    /**
+     * @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;
+    }
 
+    /**
+     * @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 (int i = 0; i < size_; ++i) {
+        for (size_t i = 0; i < size_; ++i) {
             data_[i] = std::sqrt(data_[i]);
         }
     }
 
-    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;
+    /*
+     * @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;
     }
 
     /**
@@ -224,10 +468,12 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
      */
     NDView<T, Ndim> view() const { return NDView<T, Ndim>{data_, shape_}; }
 
-    void Print();
-    void Print_all();
-    void Print_some();
-
+  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 reset() {
         data_ = nullptr;
         size_ = 0;
@@ -236,167 +482,10 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
     }
 };
 
-// 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;
-}
+///////////////////////////////////////////////////////////////////////////////
+// Free functions closely related to NDArray
+//
+///////////////////////////////////////////////////////////////////////////////
 
 template <typename T, ssize_t Ndim>
 std::ostream &operator<<(std::ostream &os, const NDArray<T, Ndim> &arr) {
@@ -410,27 +499,9 @@ 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>
-void save(NDArray<T, Ndim> &img, std::string &pathname) {
+[[deprecated("Saving of raw arrays without metadata is deprecated")]] 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));
@@ -438,8 +509,9 @@ void save(NDArray<T, Ndim> &img, std::string &pathname) {
 }
 
 template <typename T, ssize_t Ndim>
-NDArray<T, Ndim> load(const std::string &pathname,
-                      std::array<ssize_t, Ndim> shape) {
+[[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> img{shape};
     std::ifstream f;
     f.open(pathname, std::ios::binary);
@@ -448,6 +520,20 @@ NDArray<T, Ndim> load(const std::string &pathname,
     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) {
@@ -467,4 +553,4 @@ NDArray<RT, Ndim> safe_divide(const NDArray<NT, Ndim> &numerator,
     return result;
 }
 
-} // namespace aare
+} // namespace aare

include/aare/NDView.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include "aare/ArrayExpr.hpp"
 #include "aare/defs.hpp"
@@ -53,16 +54,6 @@ size_t num_elements(const Shape<Ndim> &shape) {
                            std::multiplies<size_t>());
 }
 
-template <ssize_t Dim = 0, typename Strides>
-ssize_t element_offset(const Strides & /*unused*/) {
-    return 0;
-}
-
-template <ssize_t Dim = 0, typename Strides, typename... Ix>
-ssize_t element_offset(const Strides &strides, ssize_t i, Ix... index) {
-    return i * strides[Dim] + element_offset<Dim + 1>(strides, index...);
-}
-
 template <ssize_t Ndim>
 std::array<ssize_t, Ndim> c_strides(const std::array<ssize_t, Ndim> &shape) {
     std::array<ssize_t, Ndim> strides{};
@@ -82,7 +73,8 @@ std::array<ssize_t, Ndim> make_array(const std::vector<ssize_t> &vec) {
 }
 
 template <typename T, ssize_t Ndim = 2>
-class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
+class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim>,
+               public NDIndexOps<NDView<T, Ndim>, T, Ndim> {
   public:
     NDView() = default;
     ~NDView() = default;
@@ -93,26 +85,9 @@ class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
         : buffer_(buffer), strides_(c_strides<Ndim>(shape)), shape_(shape),
           size_(std::accumulate(std::begin(shape), std::end(shape), 1,
                                 std::multiplies<>())) {}
-                   
-    template <typename... Ix>
-    std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) {
-        return buffer_[element_offset(strides_, index...)];
-    }
 
-    template <typename... Ix>
-    std::enable_if_t<sizeof...(Ix) == 1 && (Ndim > 1), NDView<T, Ndim - 1>> operator()(Ix... index) {
-        // return a view of the next dimension
-        std::array<ssize_t, Ndim - 1> new_shape{};
-        std::copy_n(shape_.begin() + 1, Ndim - 1, new_shape.begin());
-        return NDView<T, Ndim - 1>(&buffer_[element_offset(strides_, index...)],
-                                   new_shape);
-        
-    }
-
-    template <typename... Ix>
-    std::enable_if_t<sizeof...(Ix) == Ndim, const T &> operator()(Ix... index) const {
-        return buffer_[element_offset(strides_, index...)];
-    }
+    using NDIndexOps<NDView<T, Ndim>, T, Ndim>::operator();
+    using NDIndexOps<NDView<T, Ndim>, T, Ndim>::operator[];
 
 
     ssize_t size() const { return static_cast<ssize_t>(size_); }
@@ -128,16 +103,6 @@ class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
 
 
 
-    /**
-     * @brief Access element at index i.
-     */
-    T &operator[](ssize_t i)  { return buffer_[i]; }
-
-    /**
-     * @brief Access element at index i.
-     */
-    const T &operator[](ssize_t i) const { return buffer_[i]; }
-
     bool operator==(const NDView &other) const {
         if (size_ != other.size_)
             return false;
@@ -269,4 +234,4 @@ template <typename T> NDView<T, 1> make_view(std::vector<T> &vec) {
     return NDView<T, 1>(vec.data(), {static_cast<ssize_t>(vec.size())});
 }
 
-} // namespace aare
+} // namespace aare

include/aare/NumpyFile.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include "aare/Dtype.hpp"
 #include "aare/FileInterface.hpp"

include/aare/NumpyHelpers.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 
 #pragma once
 #include <algorithm>

include/aare/Pedestal.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include "aare/Frame.hpp"
 #include "aare/NDArray.hpp"

include/aare/PixelMap.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 
 #include "aare/NDArray.hpp"

include/aare/ProducerConsumerQueue.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 /*
  * Copyright (c) Meta Platforms, Inc. and affiliates.
  *

include/aare/RawFile.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include "aare/DetectorGeometry.hpp"
 #include "aare/FileInterface.hpp"

include/aare/RawMasterFile.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include "aare/defs.hpp"
 #include <algorithm>
@@ -5,6 +6,7 @@
 #include <fmt/format.h>
 #include <fstream>
 #include <optional>
+#include <chrono>
 
 #include <nlohmann/json.hpp>
 using json = nlohmann::json;
@@ -83,6 +85,9 @@ 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};
 
@@ -102,11 +107,13 @@ class RawMasterFile {
     std::optional<size_t> m_digital_samples;
     std::optional<size_t> m_transceiver_samples;
     std::optional<size_t> m_number_of_rows;
+    std::optional<uint8_t> m_counter_mask;
 
     std::optional<ROI> m_roi;
 
   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;
 
@@ -132,14 +139,18 @@ class RawMasterFile {
     std::optional<size_t> digital_samples() const;
     std::optional<size_t> transceiver_samples() const;
     std::optional<size_t> number_of_rows() const;
+    std::optional<uint8_t> counter_mask() const;
 
     std::optional<ROI> roi() const;
 
     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(const std::filesystem::path &fpath);
-    void parse_raw(const std::filesystem::path &fpath);
+    void parse_json(std::istream &is);
+    void parse_raw(std::istream &is);
     void retrieve_geometry();
 };
 

include/aare/RawSubFile.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include "aare/Frame.hpp"
 #include "aare/defs.hpp"

include/aare/VarClusterFinder.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 
 #include <algorithm>
@@ -124,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_.value(i + di[k], j + dj[k]);
+            auto tmp = labeled_(i + di[k], j + dj[k]);
             if (tmp != 0)
                 neighbour_labels.push_back(tmp);
         }

include/aare/algorithm.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 
 #pragma once
 #include <aare/NDArray.hpp>
@@ -109,4 +110,19 @@ template <typename Container> bool all_equal(const Container &c) {
     return false;
 }
 
+/**
+ * linear interpolation
+ * @param bin_edge left and right bin edges
+ * @param bin_values function values at bin edges
+ * @param coord coordinate to interpolate at
+ * @return interpolated value at coord
+ */
+inline double linear_interpolation(const std::pair<double, double> &bin_edge,
+                                   const std::pair<double, double> &bin_values,
+                                   const double coord) {
+    const double bin_width = bin_edge.second - bin_edge.first;
+    return bin_values.first * (1 - (coord - bin_edge.first) / bin_width) +
+           bin_values.second * (coord - bin_edge.first) / bin_width;
+}
+
 } // namespace aare

include/aare/calibration.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 
 #include "aare/NDArray.hpp"

include/aare/decode.hpp

@@ -1,10 +1,12 @@
+// 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,
@@ -12,6 +14,25 @@ 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

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 
 #include "aare/Dtype.hpp"
@@ -331,31 +332,37 @@ enum DACIndex {
     SLOW_ADC_TEMP
 };
 
+// helper pair class to easily expose in python
+template <typename T1, typename T2> struct Sum_index_pair {
+    T1 sum;
+    T2 index;
+};
+
+enum class corner : int {
+    cTopLeft = 0,
+    cTopRight = 1,
+    cBottomLeft = 2,
+    cBottomRight = 3
+};
+
 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
 
-/**
- * @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);
+
+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

include/aare/logger.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 /*Utility to log to console*/
 
@@ -105,7 +106,7 @@ class Logger {
     }
 
     std::ostringstream &Get() {
-        os << Color(m_level) << "- " << Timestamp() << " " << ToString(m_level)
+        os << Color(m_level) << "- " << Timestamp() << " " << Logger::ToString(m_level)
            << ": ";
         return os;
     }

include/aare/utils/ifstream_helpers.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 
 #include <fstream>

include/aare/utils/par.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include <thread>
 #include <utility>

include/aare/utils/task.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include <utility>
 #include <vector>

pyproject.toml

@@ -16,6 +16,7 @@ dependencies = [
     "numpy",
     "matplotlib",
 ]
+license = { file = "LICENSE" }
 
 
 [tool.cibuildwheel]

python/CMakeLists.txt

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 
 find_package (Python 3.10 COMPONENTS Interpreter Development.Module REQUIRED)
 set(PYBIND11_FINDPYTHON ON) # Needed for RH8
@@ -31,7 +32,7 @@ set( PYTHON_FILES
     aare/CtbRawFile.py
     aare/ClusterFinder.py
     aare/ClusterVector.py
-
+    aare/Cluster.py
     aare/calibration.py
     aare/func.py
     aare/RawFile.py

python/aare/Cluster.py

@@ -0,0 +1,24 @@
+from . import _aare 
+import numpy as np
+from .ClusterFinder import _type_to_char
+
+
+def Cluster(x : int, y : int, data, cluster_size=(3,3), dtype = np.int32):
+    """
+    Factory function to create a Cluster object. Provides a cleaner syntax for
+    the templated Cluster in C++.
+
+    .. code-block:: python
+
+        from aare import Cluster
+        
+        Cluster(cluster_size=(3,3), dtype=np.float64)
+    """
+
+    try:
+        class_name = f"Cluster{cluster_size[0]}x{cluster_size[1]}{_type_to_char(dtype)}"
+        cls = getattr(_aare, class_name)
+    except AttributeError:
+        raise ValueError(f"Unsupported combination of type and cluster size: {dtype}/{cluster_size} when requesting {class_name}")
+   
+    return cls(x, y, data) 

python/aare/ClusterFinder.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 from . import _aare
 import numpy as np
 
@@ -10,6 +11,8 @@ def _type_to_char(dtype):
         return 'f'
     elif dtype == np.float64:
         return 'd'
+    elif dtype == np.int16:
+        return 'i16'
     else:
         raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32, np.float32, and np.float64 are supported.")
 
@@ -26,7 +29,7 @@ def _get_class(name, cluster_size, dtype):
 
 
 
-def ClusterFinder(image_size, cluster_size, n_sigma=5, dtype = np.int32, capacity = 1024):
+def ClusterFinder(image_size, cluster_size=(3,3), n_sigma=5, dtype = np.int32, capacity = 1024):
     """
     Factory function to create a ClusterFinder object. Provides a cleaner syntax for 
     the templated ClusterFinder in C++.
@@ -65,7 +68,7 @@ def ClusterFileSink(clusterfindermt, cluster_file, dtype=np.int32):
     return cls(clusterfindermt, cluster_file)
 
 
-def ClusterFile(fname, cluster_size=(3,3), dtype=np.int32, chunk_size = 1000):
+def ClusterFile(fname, cluster_size=(3,3), dtype=np.int32, chunk_size = 1000, mode = "r"):
     """
     Factory function to create a ClusterFile object. Provides a cleaner syntax for
     the templated ClusterFile in C++.
@@ -83,4 +86,4 @@ def ClusterFile(fname, cluster_size=(3,3), dtype=np.int32, chunk_size = 1000):
     """
 
     cls = _get_class("ClusterFile", cluster_size, dtype)
-    return cls(fname, chunk_size=chunk_size)
+    return cls(fname, chunk_size=chunk_size, mode=mode)

python/aare/ClusterVector.py

@@ -1,11 +1,22 @@
+# SPDX-License-Identifier: MPL-2.0
 
 
-from ._aare import ClusterVector_Cluster3x3i
+from . import _aare 
 import numpy as np
+from .ClusterFinder import _get_class
 
-def ClusterVector(cluster_size, dtype = np.int32):
+def ClusterVector(cluster_size=(3,3), dtype = np.int32):
+    """
+    Factory function to create a ClusterVector object. Provides a cleaner syntax for
+    the templated ClusterVector in C++.
+
+    .. code-block:: python
+
+        from aare import ClusterVector
+        
+        ClusterVector(cluster_size=(3,3), dtype=np.float64)
+    """
+
+    cls = _get_class("ClusterVector", cluster_size, dtype)
+    return cls()
 
-    if dtype == np.int32 and cluster_size == (3,3):
-        return ClusterVector_Cluster3x3i()
-    else:
-        raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.")

python/aare/CtbRawFile.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 
 from . import _aare
 import numpy as np

python/aare/RawFile.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 from . import _aare
 import numpy as np
 from .ScanParameters import ScanParameters

python/aare/ScanParameters.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 from . import _aare
 
 class ScanParameters(_aare.ScanParameters):

python/aare/__init__.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 # Make the compiled classes that live in _aare available from aare.
 from . import _aare
 
@@ -7,16 +8,18 @@ from ._aare import Pedestal_d, Pedestal_f, ClusterFinder_Cluster3x3i, VarCluster
 from ._aare import DetectorType
 from ._aare import hitmap
 from ._aare import ROI
+from ._aare import corner 
 
 # from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i
 
 from .ClusterFinder import ClusterFinder, ClusterCollector, ClusterFinderMT, ClusterFileSink, ClusterFile
 from .ClusterVector import ClusterVector
+from .Cluster import Cluster
 
 
 from ._aare import fit_gaus, fit_pol1, fit_scurve, fit_scurve2
 from ._aare import Interpolator
-from ._aare import calculate_eta2
+from ._aare import calculate_eta2, calculate_eta3, calculate_cross_eta3, calculate_full_eta2
 from ._aare import reduce_to_2x2, reduce_to_3x3
 
 from ._aare import apply_custom_weights

python/aare/calibration.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 #Calibration related functions
 import numpy as np
 def load_calibration(fname, hg0=False):

python/aare/func.py

@@ -1 +1,2 @@
+# SPDX-License-Identifier: MPL-2.0
 from ._aare import gaus, pol1, scurve, scurve2

python/aare/transform.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 import numpy as np
 from . import _aare
 
@@ -48,6 +49,43 @@ 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/aare/utils.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 import numpy as np
 import matplotlib.pyplot as plt
 from mpl_toolkits.axes_grid1 import make_axes_locatable

python/examples/fits.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 import matplotlib.pyplot as plt
 import numpy as np
 from aare import fit_gaus, fit_pol1

python/examples/play.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 
 from aare import apply_calibration
 import numpy as np

python/src/bind_Cluster.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #include "aare/Cluster.hpp"
 
 #include <cstdint>
@@ -58,7 +59,17 @@ void define_Cluster(py::module &m, const std::string &typestr) {
                       &Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>::x)
 
         .def_readonly("y",
-                      &Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>::y);
+                      &Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>::y)
+
+        .def(
+            "max_sum_2x2",
+            [](Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType> &self) {
+                auto max_sum = self.max_sum_2x2();
+                return py::make_tuple(max_sum.sum,
+                                      static_cast<int>(max_sum.index));
+            },
+            R"(calculates sum of 2x2 subcluster with highest energy and index relative to cluster center 0: top_left, 1: top_right, 2: bottom_left, 3: bottom_right
+            )");
 }
 
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
@@ -70,9 +81,7 @@ void reduce_to_3x3(py::module &m) {
         [](const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
             return reduce_to_3x3(cl);
         },
-        py::return_value_policy::move,
-        "Reduce cluster to 3x3 subcluster by taking the 3x3 subcluster with "
-        "the highest photon energy.");
+        py::return_value_policy::move, R"(Reduce cluster to 3x3 subcluster)");
 }
 
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
@@ -85,8 +94,15 @@ void reduce_to_2x2(py::module &m) {
             return reduce_to_2x2(cl);
         },
         py::return_value_policy::move,
-        "Reduce cluster to 2x2 subcluster by taking the 2x2 subcluster with "
-        "the highest photon energy.");
+        R"(
+        Reduce cluster to 2x2 subcluster by taking the 2x2 subcluster with 
+        the highest photon energy.
+
+        RETURN: 
+        
+        reduced cluster (cluster is filled in row major ordering starting at the top left. Thus for a max subcluster in the top left corner the photon hit is at the fourth position.)
+        
+        )");
 }
 
 #pragma GCC diagnostic pop

python/src/bind_ClusterCollector.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #include "aare/ClusterCollector.hpp"
 #include "aare/ClusterFileSink.hpp"
 #include "aare/ClusterFinder.hpp"

python/src/bind_ClusterFile.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #include "aare/CalculateEta.hpp"
 #include "aare/ClusterFile.hpp"
 #include "aare/defs.hpp"
@@ -45,6 +46,7 @@ void define_ClusterFile(py::module &m, const std::string &typestr) {
                  return v;
              })
         .def("set_roi", &ClusterFile<ClusterType>::set_roi, py::arg("roi"))
+        .def("tell", &ClusterFile<ClusterType>::tell)
         .def(
             "set_noise_map",
             [](ClusterFile<ClusterType> &self, py::array_t<int32_t> noise_map) {
@@ -80,23 +82,4 @@ void define_ClusterFile(py::module &m, const std::string &typestr) {
         });
 }
 
-template <typename Type, uint8_t CoordSizeX, uint8_t CoordSizeY,
-          typename CoordType = uint16_t>
-void register_calculate_eta(py::module &m) {
-    using ClusterType = Cluster<Type, CoordSizeX, CoordSizeY, CoordType>;
-
-    m.def("calculate_eta2",
-          [](const aare::ClusterVector<ClusterType> &clusters) {
-              auto eta2 = new NDArray<double, 2>(calculate_eta2(clusters));
-              return return_image_data(eta2);
-          });
-
-    m.def("calculate_eta2", [](const aare::Cluster<Type, CoordSizeX, CoordSizeY,
-                                                   CoordType> &cluster) {
-        auto eta2 = calculate_eta2(cluster);
-        // TODO return proper eta class
-        return py::make_tuple(eta2.x, eta2.y, eta2.sum);
-    });
-}
-
 #pragma GCC diagnostic pop