Release 2025.11.21 (#249)

- Updated VERSION and script 
- Updated release notes

CMakeLists.txt

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 cmake_minimum_required(VERSION 3.15)
 
 project(aare 
@@ -388,7 +389,7 @@ set(SourceFiles
     ${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/utils/task.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/utils/ifstream_helpers.cpp
-)
+) 
 
 add_library(aare_core STATIC ${SourceFiles})
 target_include_directories(aare_core PUBLIC 
@@ -412,6 +413,8 @@ target_link_libraries(
 
 )
 
+set_property(TARGET aare_core PROPERTY POSITION_INDEPENDENT_CODE ON)
+
 if(AARE_TESTS)
     target_compile_definitions(aare_core PRIVATE AARE_TESTS)
 endif()
@@ -431,10 +434,6 @@ set_target_properties(aare_core PROPERTIES
     PUBLIC_HEADER "${PUBLICHEADERS}"
 )
 
-if (AARE_PYTHON_BINDINGS)
-    set_property(TARGET aare_core PROPERTY POSITION_INDEPENDENT_CODE ON)
-endif()
-
 if(AARE_TESTS)
     set(TestSources
             ${CMAKE_CURRENT_SOURCE_DIR}/src/algorithm.test.cpp
@@ -444,6 +443,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
@@ -465,6 +465,7 @@ if(AARE_TESTS)
     target_sources(tests PRIVATE ${TestSources} )
 endif()
 
+
 if(AARE_MASTER_PROJECT)
     install(TARGETS aare_core aare_compiler_flags 
         EXPORT "${TARGETS_EXPORT_NAME}"
@@ -474,7 +475,6 @@ if(AARE_MASTER_PROJECT)
     )
 endif()
 
-set(CMAKE_POSITION_INDEPENDENT_CODE ON)
 set(CMAKE_INSTALL_RPATH $ORIGIN)
 set(CMAKE_BUILD_WITH_INSTALL_RPATH FALSE)
 

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,16 +1,56 @@
 # Release notes
 
 
-### head
+
+### 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:
 
 - Apply calibration works in G0 if passes a 2D calibration and pedestal
 - count pixels that switch
 - calculate pedestal (also g0 version)
+- NDArray::view() needs an lvalue to reduce issues with the view outliving the array
 
 
-### 2025.07.18
+Bugfixes:
+
+- Now using glibc 2.17 in conda builds (was using the host)
+- Fixed shifted pixels in clusters close to the edge of a frame
+
+### 2025.7.18
 
 Features:
 
@@ -24,7 +64,7 @@ Bugfixes:
 - Removed unused file: ClusterFile.cpp 
 
 
-### 2025.05.22
+### 2025.5.22
 
 Features:
 
@@ -34,6 +74,34 @@ 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.7.18
+2025.11.21

benchmarks/CMakeLists.txt

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 
 include(FetchContent)
 
@@ -15,7 +16,7 @@ FetchContent_MakeAvailable(benchmark)
 
 add_executable(benchmarks)
 
-target_sources(benchmarks PRIVATE ndarray_benchmark.cpp calculateeta_benchmark.cpp)
+target_sources(benchmarks PRIVATE ndarray_benchmark.cpp calculateeta_benchmark.cpp reduce_benchmark.cpp)
 
 # Link Google Benchmark and other necessary libraries
 target_link_libraries(benchmarks PRIVATE benchmark::benchmark aare_core aare_compiler_flags)

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

@@ -0,0 +1,169 @@
+// SPDX-License-Identifier: MPL-2.0
+#include "aare/Cluster.hpp"
+#include <benchmark/benchmark.h>
+
+using namespace aare;
+
+class ClustersForReduceFixture : public benchmark::Fixture {
+  public:
+    Cluster<int, 5, 5> cluster_5x5{};
+    Cluster<int, 3, 3> cluster_3x3{};
+
+  private:
+    using benchmark::Fixture::SetUp;
+
+    void SetUp([[maybe_unused]] const benchmark::State &state) override {
+        int temp_data[25] = {1, 1, 1, 1, 1, 1, 1, 2, 1, 1,
+                             1, 2, 3, 1, 2, 1, 1, 1, 1, 2};
+        std::copy(std::begin(temp_data), std::end(temp_data),
+                  std::begin(cluster_5x5.data));
+
+        cluster_5x5.x = 5;
+        cluster_5x5.y = 5;
+
+        int temp_data2[9] = {1, 1, 1, 2, 3, 1, 2, 2, 1};
+        std::copy(std::begin(temp_data2), std::end(temp_data2),
+                  std::begin(cluster_3x3.data));
+
+        cluster_3x3.x = 5;
+        cluster_3x3.y = 5;
+    }
+
+    // void TearDown(::benchmark::State& state) {
+    // }
+};
+
+template <typename T>
+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;
+
+    // Write out the sums in the hope that the compiler can optimize this
+    sum_3x3_subclusters[0] = c.data[0] + c.data[1] + c.data[2] + c.data[5] +
+                             c.data[6] + c.data[7] + c.data[10] + c.data[11] +
+                             c.data[12];
+    sum_3x3_subclusters[1] = c.data[1] + c.data[2] + c.data[3] + c.data[6] +
+                             c.data[7] + c.data[8] + c.data[11] + c.data[12] +
+                             c.data[13];
+    sum_3x3_subclusters[2] = c.data[2] + c.data[3] + c.data[4] + c.data[7] +
+                             c.data[8] + c.data[9] + c.data[12] + c.data[13] +
+                             c.data[14];
+    sum_3x3_subclusters[3] = c.data[5] + c.data[6] + c.data[7] + c.data[10] +
+                             c.data[11] + c.data[12] + c.data[15] + c.data[16] +
+                             c.data[17];
+    sum_3x3_subclusters[4] = c.data[6] + c.data[7] + c.data[8] + c.data[11] +
+                             c.data[12] + c.data[13] + c.data[16] + c.data[17] +
+                             c.data[18];
+    sum_3x3_subclusters[5] = c.data[7] + c.data[8] + c.data[9] + c.data[12] +
+                             c.data[13] + c.data[14] + c.data[17] + c.data[18] +
+                             c.data[19];
+    sum_3x3_subclusters[6] = c.data[10] + c.data[11] + c.data[12] + c.data[15] +
+                             c.data[16] + c.data[17] + c.data[20] + c.data[21] +
+                             c.data[22];
+    sum_3x3_subclusters[7] = c.data[11] + c.data[12] + c.data[13] + c.data[16] +
+                             c.data[17] + c.data[18] + c.data[21] + c.data[22] +
+                             c.data[23];
+    sum_3x3_subclusters[8] = c.data[12] + c.data[13] + c.data[14] + c.data[17] +
+                             c.data[18] + c.data[19] + c.data[22] + c.data[23] +
+                             c.data[24];
+
+    auto index = std::max_element(sum_3x3_subclusters.begin(),
+                                  sum_3x3_subclusters.end()) -
+                 sum_3x3_subclusters.begin();
+
+    switch (index) {
+    case 0:
+        result.x = c.x - 1;
+        result.y = c.y + 1;
+        result.data = {c.data[0], c.data[1],  c.data[2],  c.data[5], c.data[6],
+                       c.data[7], c.data[10], c.data[11], c.data[12]};
+        break;
+    case 1:
+        result.x = c.x;
+        result.y = c.y + 1;
+        result.data = {c.data[1], c.data[2],  c.data[3],  c.data[6], c.data[7],
+                       c.data[8], c.data[11], c.data[12], c.data[13]};
+        break;
+    case 2:
+        result.x = c.x + 1;
+        result.y = c.y + 1;
+        result.data = {c.data[2], c.data[3],  c.data[4],  c.data[7], c.data[8],
+                       c.data[9], c.data[12], c.data[13], c.data[14]};
+        break;
+    case 3:
+        result.x = c.x - 1;
+        result.y = c.y;
+        result.data = {c.data[5],  c.data[6],  c.data[7],
+                       c.data[10], c.data[11], c.data[12],
+                       c.data[15], c.data[16], c.data[17]};
+        break;
+    case 4:
+        result.x = c.x + 1;
+        result.y = c.y;
+        result.data = {c.data[6],  c.data[7],  c.data[8],
+                       c.data[11], c.data[12], c.data[13],
+                       c.data[16], c.data[17], c.data[18]};
+        break;
+    case 5:
+        result.x = c.x + 1;
+        result.y = c.y;
+        result.data = {c.data[7],  c.data[8],  c.data[9],
+                       c.data[12], c.data[13], c.data[14],
+                       c.data[17], c.data[18], c.data[19]};
+        break;
+    case 6:
+        result.x = c.x + 1;
+        result.y = c.y - 1;
+        result.data = {c.data[10], c.data[11], c.data[12],
+                       c.data[15], c.data[16], c.data[17],
+                       c.data[20], c.data[21], c.data[22]};
+        break;
+    case 7:
+        result.x = c.x + 1;
+        result.y = c.y - 1;
+        result.data = {c.data[11], c.data[12], c.data[13],
+                       c.data[16], c.data[17], c.data[18],
+                       c.data[21], c.data[22], c.data[23]};
+        break;
+    case 8:
+        result.x = c.x + 1;
+        result.y = c.y - 1;
+        result.data = {c.data[12], c.data[13], c.data[14],
+                       c.data[17], c.data[18], c.data[19],
+                       c.data[22], c.data[23], c.data[24]};
+        break;
+    }
+    return result;
+}
+
+BENCHMARK_F(ClustersForReduceFixture, Reduce2x2)(benchmark::State &st) {
+    for (auto _ : st) {
+        // This code gets timed
+        benchmark::DoNotOptimize(reduce_to_2x2<int, 3, 3, uint16_t>(
+            cluster_3x3)); // make sure compiler evaluates the expression
+    }
+}
+
+BENCHMARK_F(ClustersForReduceFixture, SpecificReduce2x2)(benchmark::State &st) {
+    for (auto _ : st) {
+        // This code gets timed
+        benchmark::DoNotOptimize(reduce_to_2x2<int>(cluster_3x3));
+    }
+}
+
+BENCHMARK_F(ClustersForReduceFixture, Reduce3x3)(benchmark::State &st) {
+    for (auto _ : st) {
+        // This code gets timed
+        benchmark::DoNotOptimize(
+            reduce_to_3x3<int, 5, 5, uint16_t>(cluster_5x5));
+    }
+}
+
+BENCHMARK_F(ClustersForReduceFixture, SpecificReduce3x3)(benchmark::State &st) {
+    for (auto _ : st) {
+        // This code gets timed
+        benchmark::DoNotOptimize(reduce_to_3x3<int>(cluster_5x5));
+    }
+}

conda-recipe/conda_build_config.yaml

@@ -3,3 +3,14 @@ python:
     - 3.12
     - 3.13
 
+c_compiler:
+  - gcc                        # [linux]
+
+c_stdlib:
+  - sysroot                    # [linux]
+
+cxx_compiler:
+  - gxx                        # [linux]
+
+c_stdlib_version:              # [linux]
+  - 2.17                       # [linux]

conda-recipe/meta.yaml

@@ -16,6 +16,8 @@ build:
 
 requirements:
   build:
+    - {{ compiler('c') }}
+    - {{ stdlib("c") }}
     - {{ compiler('cxx') }}
     - cmake
     - ninja
@@ -47,4 +49,5 @@ test:
     - 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)
 
@@ -27,6 +28,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/ClusterVector.rst

@@ -12,4 +12,11 @@ ClusterVector
    :members:
    :undoc-members:
    :private-members:
+   
+
+**Free Functions:** 
+
+.. doxygenfunction:: aare::reduce_to_3x3(const ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>&)
+   
+.. doxygenfunction:: aare::reduce_to_2x2(const ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>&)
    

docs/src/Interpolation.rst

@@ -0,0 +1,102 @@
+Interpolation
+==============
+
+Interpolation class for :math:`\eta` Interpolation. 
+
+The Interpolator class provides methods to interpolate the positions of photons based on their :math:`\eta` values. 
+
+.. warning:: 
+   The interpolation might lead to erroneous photon positions for clusters at the boarders of a frame. Make sure to filter out such cases. 
+
+: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: 
+
+
+.. 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*}
+
+
+.. doxygenfunction:: aare::calculate_eta2(const ClusterVector<ClusterType>&)
+
+.. doxygenfunction:: aare::calculate_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>&) 
+
+.. image:: ../figures/Eta2x2Full.png
+    :target: ../figures/Eta2x2Full.png
+    :width: 650px
+    :align: center
+    :alt: Eta2x2 Full 
+
+.. math:: 
+    \begin{equation*}
+        {\color{blue}{\eta_x}}  = \frac{Q_{0,1} + Q_{1,1}}{\sum_i^{2}\sum_j^{2}Q_{i,j}} \quad \quad
+        {\textcolor{green}{\eta_y}} = \frac{Q_{1,0} + Q_{1,1}}{\sum_i^{2}\sum_j^{2}Q_{i,j}}
+    \end{equation*} 
+
+
+.. doxygenfunction:: aare::calculate_full_eta2(const ClusterVector<ClusterType>&)
+
+.. doxygenfunction:: aare::calculate_full_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>&)
+
+.. image:: ../figures/Eta3x3.png
+    :target: ../figures/Eta3x3.png
+    :width: 650px
+    :align: center
+    :alt: Eta3x3
+
+.. math::
+    \begin{equation*}
+         {\color{blue}{\eta_x}} = \frac{\sum_{i}^{3} Q_{i,2} - \sum_{i}^{3} Q_{i,0}}{\sum_{i}^{3}\sum_{j}^{3} Q_{i,j}} \quad \quad
+        {\color{green}{\eta_y}} = \frac{\sum_{j}^{3} Q_{2,j} - \sum_{j}^{3} Q_{0,j}}{\sum_{i}^{3}\sum_{j}^{3} Q_{i,j}}
+    \end{equation*}
+
+.. doxygenfunction:: aare::calculate_eta3(const ClusterVector<Cluster<T, 3,3, CoordType>>&)
+
+.. doxygenfunction:: aare::calculate_eta3(const Cluster<T, 3, 3, CoordType>&)
+
+.. image:: ../figures/Eta3x3Cross.png
+    :target: ../figures/Eta3x3Cross.png
+    :width: 650px
+    :align: center
+    :alt: Cross Eta3x3
+
+.. math:: 
+    \begin{equation*}
+       {\color{blue}{\eta_x}} = \frac{Q_{1,2} - Q_{1,0}}{Q_{1,0} +  Q_{1,1} + Q_{1,0}} \quad \quad
+        {\color{green}{\eta_y}} = \frac{Q_{0,2} - Q_{0,1}}{Q_{0,1} +  Q_{1,1} + Q_{1,2}}
+    \end{equation*}
+
+.. doxygenfunction:: aare::calculate_cross_eta3(const ClusterVector<Cluster<T, 3,3, CoordType>>&)
+
+.. doxygenfunction:: aare::calculate_cross_eta3(const Cluster<T, 3, 3, CoordType>&)
+
+Interpolation class: 
+---------------------
+
+.. Warning:: 
+    Make sure to use the same :math:`\eta`-function during interpolation as given by the joint :math:`\eta`-distribution passed to the constructor. 
+
+.. doxygenclass:: aare::Interpolator
+   :members:
+   :undoc-members:
+   :private-members:
+
+

docs/src/index.rst

@@ -29,6 +29,8 @@ AARE
     pyCtbRawFile
     pyClusterFile
     pyClusterVector
+    pyCluster
+    pyInterpolation
     pyJungfrauDataFile
     pyRawFile
     pyRawMasterFile
@@ -47,10 +49,12 @@ AARE
     Frame
     File
     Dtype
+    Cluster
     ClusterFinder
     ClusterFinderMT
     ClusterFile
     ClusterVector
+    Interpolation
     JungfrauDataFile
     Pedestal
     RawFile

docs/src/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/pyClusterVector.rst

@@ -1,3 +1,5 @@
+.. _py_clustervector:
+    
 ClusterVector
 ================
 
@@ -28,9 +30,29 @@ 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:
     :undoc-members:
     :show-inheritance:
-    :inherited-members:
+    :inherited-members:
+
+
+**Free Functions:** 
+
+.. autofunction:: reduce_to_3x3
+   :noindex:
+
+   Reduce a single Cluster to 3x3 by taking the 3x3 subcluster with highest photon energy.
+
+.. autofunction:: reduce_to_2x2
+   :noindex:
+
+   Reduce a single Cluster to 2x2 by taking the 2x2 subcluster with highest photon energy.

docs/src/pyInterpolation.rst

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

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

@@ -13,4 +13,5 @@ dependencies:
   - pybind11
   - numpy
   - matplotlib
+  - nlohmann_json
 

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.
@@ -7,11 +7,13 @@ Script to update VERSION file with semantic versioning if provided as an argumen
 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:
@@ -26,9 +28,9 @@ def get_version():
 
     # Check at least one argument is passed
     if len(sys.argv) < 2:
-        return "0.0.0"
-
-    version = sys.argv[1]
+        version =   datetime.today().strftime('%Y.%-m.%-d')
+    else: 
+        version = sys.argv[1]
 
     try:
         v = Version(version)  # normalize check if version follows PEP 440 specification
@@ -45,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.")
@@ -54,4 +57,4 @@ def write_version_to_file(version):
 if __name__ == "__main__":
 
     version = get_version()
-    write_version_to_file(version)
+    write_version_to_file(version)

include/aare/ArrayExpr.hpp

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

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 {
-    cBottomLeft = 0,
-    cBottomRight = 1,
-    cTopLeft = 2,
-    cTopRight = 3
-};
-
 enum class pixel : int {
     pBottomLeft = 0,
     pBottom = 1,
@@ -25,146 +20,427 @@ 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;
-    T sum;
+    /// @brief eta in x direction
+    double x{};
+    /// @brief eta in y direction
+    double y{};
+    /// @brief index of subcluster 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) {
-    Eta2<T> eta{};
 
-    auto max_sum = cl.max_sum_2x2();
-    eta.sum = max_sum.first;
-    auto c = max_sum.second;
+    static_assert(ClusterSizeX > 1 && ClusterSizeY > 1);
+    Eta2<T> eta{};
 
     size_t cluster_center_index =
         (ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
 
-    size_t index_bottom_left_max_2x2_subcluster =
-        (int(c / (ClusterSizeX - 1))) * ClusterSizeX + c % (ClusterSizeX - 1);
+    auto max_sum = cl.max_sum_2x2();
+    eta.sum = max_sum.sum;
+    corner c = max_sum.index;
 
-    // check that cluster center is in max subcluster
-    if (cluster_center_index != index_bottom_left_max_2x2_subcluster &&
-        cluster_center_index != index_bottom_left_max_2x2_subcluster + 1 &&
-        cluster_center_index !=
-            index_bottom_left_max_2x2_subcluster + ClusterSizeX &&
-        cluster_center_index !=
-            index_bottom_left_max_2x2_subcluster + ClusterSizeX + 1)
-        throw std::runtime_error("Photon center is not in max 2x2_subcluster");
-
-    if ((cluster_center_index - index_bottom_left_max_2x2_subcluster) %
-            ClusterSizeX ==
-        0) {
-        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]));
-    } else {
-        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 - 1] +
-                                         cl.data[cluster_center_index]));
+    // subcluster top right from center
+    switch (c) {
+    case (corner::cTopLeft):
+        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):
+        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):
+        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):
+        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;
     }
-    if ((cluster_center_index - index_bottom_left_max_2x2_subcluster) /
-            ClusterSizeX <
-        1) {
-        assert(cluster_center_index + ClusterSizeX <
-               ClusterSizeX * ClusterSizeY); // suppress warning
-        if ((cl.data[cluster_center_index] +
-             cl.data[cluster_center_index + ClusterSizeX]) != 0)
-            eta.y = static_cast<double>(
+
+    eta.c = c;
+
+    return eta;
+}
+
+/**
+ * @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{};
+
+    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>(
-                        (cl.data[cluster_center_index] +
-                         cl.data[cluster_center_index + ClusterSizeX]));
-    } else {
-        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]));
+                    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; // TODO only supported for 2x2 and 3x3 clusters -> at least no
-               // underyling enum class
+    eta.c = c;
+
     return eta;
 }
 
-// TODO! Look up eta2 calculation - photon center should be top right corner
 template <typename T>
-Eta2<T> calculate_eta2(const Cluster<T, 2, 2, int16_t> &cl) {
+Eta2<T> calculate_eta2(const Cluster<T, 2, 2, uint16_t> &cl) {
     Eta2<T> eta{};
 
-    if ((cl.data[0] + cl.data[1]) != 0)
-        eta.x = static_cast<double>(cl.data[1]) / (cl.data[0] + cl.data[1]);
-    if ((cl.data[0] + cl.data[2]) != 0)
-        eta.y = static_cast<double>(cl.data[2]) / (cl.data[0] + cl.data[2]);
+    // 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();
-    eta.c = static_cast<int>(corner::cBottomLeft); // TODO! This is not correct,
-                                                   // but need to put something
+
     return eta;
 }
 
-// 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) {
+template <typename T>
+Eta2<T> calculate_full_eta2(const Cluster<T, 2, 2, uint16_t> &cl) {
 
     Eta2<T> eta{};
 
-    T sum = 0;
+    eta.sum = cl.sum();
 
-    std::for_each(std::begin(cl.data), std::end(cl.data),
-                  [&sum](T x) { sum += x; });
+    const uint8_t photon_hit_index =
+        std::max_element(cl.data.begin(), cl.data.end()) - cl.data.begin();
 
-    eta.sum = sum;
+    eta.c = static_cast<corner>(3 - photon_hit_index);
 
-    eta.c = corner::cBottomLeft;
+    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, uint16_t> &cl) {
+    Eta2<T> eta{};
+
+    eta.x = 0;
+    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, uint16_t> &cl) {
+    Eta2<T> eta{};
+
+    eta.x = static_cast<double>(cl.data[1]) / cl.data[0];
+    eta.y = 0;
+    eta.sum = cl.sum();
+}
+
+/**
+ * @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 photon_energy = cl.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]);
+            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,8 +19,12 @@
 namespace aare {
 
 // requires clause c++20 maybe update
+
+/**
+ * @brief Cluster struct
+ */
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType = int16_t>
+          typename CoordType = uint16_t>
 struct Cluster {
 
     static_assert(std::is_arithmetic_v<T>, "T needs to be an arithmetic type");
@@ -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,9 +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{}); }
 
-    std::pair<T, int> max_sum_2x2() const {
+    // TODO: handle 1 dimensional clusters
+    /**
+     * @brief sum of 2x2 subcluster with highest energy
+     * @return photon energy of subcluster, 2x2 subcluster index relative to
+     * cluster center
+     */
+    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;
@@ -49,31 +67,166 @@ 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 num_2x2_subclusters =
-                (ClusterSizeX - 1) * (ClusterSizeY - 1);
+            constexpr size_t cluster_center_index =
+                (ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
 
-            std::array<T, num_2x2_subclusters> sum_2x2_subcluster;
-            for (size_t i = 0; i < ClusterSizeY - 1; ++i) {
-                for (size_t j = 0; j < ClusterSizeX - 1; ++j)
-                    sum_2x2_subcluster[i * (ClusterSizeX - 1) + j] =
-                        data[i * ClusterSizeX + j] +
-                        data[i * ClusterSizeX + j + 1] +
-                        data[(i + 1) * ClusterSizeX + j] +
-                        data[(i + 1) * ClusterSizeX + j + 1];
+            std::array<T, 4> sum_2x2_subcluster{0};
+            // subcluster top left from center
+            sum_2x2_subcluster[0] =
+                data[cluster_center_index] + data[cluster_center_index - 1] +
+                data[cluster_center_index - ClusterSizeX] +
+                data[cluster_center_index - 1 - ClusterSizeX];
+            // subcluster top right from center
+            if (ClusterSizeX > 2) {
+                sum_2x2_subcluster[1] =
+                    data[cluster_center_index] +
+                    data[cluster_center_index + 1] +
+                    data[cluster_center_index - ClusterSizeX] +
+                    data[cluster_center_index - ClusterSizeX + 1];
+            }
+            // subcluster bottom left from center
+            if (ClusterSizeY > 2) {
+                sum_2x2_subcluster[2] =
+                    data[cluster_center_index] +
+                    data[cluster_center_index - 1] +
+                    data[cluster_center_index + ClusterSizeX] +
+                    data[cluster_center_index + ClusterSizeX - 1];
+            }
+            // subcluster bottom right from center
+            if (ClusterSizeX > 2 && ClusterSizeY > 2) {
+                sum_2x2_subcluster[3] =
+                    data[cluster_center_index] +
+                    data[cluster_center_index + 1] +
+                    data[cluster_center_index + ClusterSizeX] +
+                    data[cluster_center_index + ClusterSizeX + 1];
             }
 
             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}};
         }
     }
 };
 
+/**
+ * @brief Reduce a cluster to a 2x2 cluster by selecting the 2x2 block with the
+ * 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 = 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");
+
+    Cluster<T, 2, 2, CoordType> result{};
+
+    auto [sum, index] = c.max_sum_2x2();
+
+    constexpr int16_t cluster_center_index =
+        (ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
+
+    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;
+    result.y = c.y;
+
+    result.data = {
+        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, 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 corner::cTopLeft:
+        result.data = {c.data[0], c.data[1], c.data[3], c.data[4]};
+        break;
+    case corner::cTopRight:
+        result.data = {c.data[1], c.data[2], c.data[4], c.data[5]};
+        break;
+    case corner::cBottomLeft:
+        result.data = {c.data[3], c.data[4], c.data[6], c.data[7]};
+        break;
+    case corner::cBottomRight:
+        result.data = {c.data[4], c.data[5], c.data[7], c.data[8]};
+        break;
+    }
+
+    return result;
+}
+
+/**
+ * @brief Reduce a cluster to a 3x3 cluster
+ * @param c Cluster to reduce
+ * @return reduced cluster
+ */
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = int16_t>
+Cluster<T, 3, 3, CoordType>
+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{};
+
+    int16_t cluster_center_index =
+        (ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
+
+    result.x = c.x;
+    result.y = c.y;
+
+    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]};
+
+    return result;
+}
+
 // Type Traits for is_cluster_type
 template <typename T>
 struct is_cluster : std::false_type {}; // Default case: Not a Cluster

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;
@@ -19,11 +28,9 @@ class ClusterFinder {
     const PEDESTAL_TYPE c3;
     Pedestal<PEDESTAL_TYPE> m_pedestal;
     ClusterVector<ClusterType> m_clusters;
-    const uint32_t ClusterSizeX;
-    const uint32_t ClusterSizeY;
 
-    static const uint8_t SavedClusterSizeX = ClusterType::cluster_size_x;
-    static const uint8_t SavedClusterSizeY = ClusterType::cluster_size_y;
+    static const uint8_t ClusterSizeX = ClusterType::cluster_size_x;
+    static const uint8_t ClusterSizeY = ClusterType::cluster_size_y;
     using CT = typename ClusterType::value_type;
 
   public:
@@ -36,12 +43,10 @@ class ClusterFinder {
      *
      */
     ClusterFinder(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0,
-                  size_t capacity = 1000000,
-                  uint32_t cluster_size_x = 3, uint32_t cluster_size_y = 3)
+                  size_t capacity = 1000000)
         : m_image_size(image_size), m_nSigma(nSigma),
-          c2(sqrt((cluster_size_y + 1) / 2 * (cluster_size_x + 1) / 2)),
-          c3(sqrt(cluster_size_x * cluster_size_y)),
-          ClusterSizeX(cluster_size_x), ClusterSizeY(cluster_size_y),
+          c2(sqrt((ClusterSizeY + 1) / 2 * (ClusterSizeX + 1) / 2)),
+          c3(sqrt(ClusterSizeX * ClusterSizeY)),
           m_pedestal(image_size[0], image_size[1]), m_clusters(capacity) {
         LOG(logDEBUG) << "ClusterFinder: "
                       << "image_size: " << image_size[0] << "x" << image_size[1]
@@ -78,9 +83,6 @@ class ClusterFinder {
         // //  4,4 -> +/- 2
         int dy = ClusterSizeY / 2;
         int dx = ClusterSizeX / 2;
-        int dy2 = SavedClusterSizeY / 2;
-        int dx2 = SavedClusterSizeX / 2;
-
         int has_center_pixel_x =
             ClusterSizeX %
             2; // for even sized clusters there is no proper cluster center and
@@ -142,14 +144,16 @@ class ClusterFinder {
                     // It's worth redoing the look since most of the time we
                     // don't have a photon
                     int i = 0;
-                    for (int ir = -dy2; ir < dy2 + has_center_pixel_y; ir++) {
-                        for (int ic = -dx2; ic < dx2 + has_center_pixel_y; ic++) {
+                    for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) {
+                        for (int ic = -dx; ic < dx + has_center_pixel_x; ic++) {
                             if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
                                 iy + ir >= 0 && iy + ir < frame.shape(0)) {
-                                    
-                                CT tmp = static_cast<CT>(frame(iy + ir, ix + ic)) - static_cast<CT>(m_pedestal.mean(iy + ir, ix + ic));
-                                cluster.data[i] = tmp; // Watch for out of bounds access
-                                
+                                CT tmp =
+                                    static_cast<CT>(frame(iy + ir, ix + ic)) -
+                                    static_cast<CT>(
+                                        m_pedestal.mean(iy + ir, ix + ic));
+                                cluster.data[i] =
+                                    tmp; // Watch for out of bounds access
                             }
                             i++;
                         }

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:
@@ -121,8 +123,7 @@ class ClusterFinderMT {
      * @param n_threads number of threads to use
      */
     ClusterFinderMT(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0,
-                    size_t capacity = 2000, size_t n_threads = 3,
-                    uint32_t cluster_size_x = 3, uint32_t cluster_size_y = 3)
+                    size_t capacity = 2000, size_t n_threads = 3)
         : m_n_threads(n_threads) {
 
         LOG(logDEBUG1) << "ClusterFinderMT: "
@@ -135,7 +136,7 @@ class ClusterFinderMT {
             m_cluster_finders.push_back(
                 std::make_unique<
                     ClusterFinder<ClusterType, FRAME_TYPE, PEDESTAL_TYPE>>(
-                    image_size, nSigma, capacity, cluster_size_x, cluster_size_y));
+                    image_size, nSigma, capacity));
         }
         for (size_t i = 0; i < n_threads; i++) {
             m_input_queues.emplace_back(std::make_unique<InputQueue>(200));

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,12 +29,11 @@ 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>
-class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> 
-{
+class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
 
     std::vector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> m_data{};
     int32_t m_frame_number{0}; // TODO! Check frame number size and type
@@ -87,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;
     }
@@ -173,4 +172,42 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>
     }
 };
 
+/**
+ * @brief Reduce a cluster to a 2x2 cluster by selecting the 2x2 block with the
+ * 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>
+ClusterVector<Cluster<T, 2, 2, CoordType>> reduce_to_2x2(
+    const ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>
+        &cv) {
+    ClusterVector<Cluster<T, 2, 2, CoordType>> result;
+    for (const auto &c : cv) {
+        result.push_back(reduce_to_2x2(c));
+    }
+    return result;
+}
+
+/**
+ * @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>
+ClusterVector<Cluster<T, 3, 3, CoordType>> reduce_to_3x3(
+    const ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>
+        &cv) {
+    ClusterVector<Cluster<T, 3, 3, CoordType>> result;
+    for (const auto &c : cv) {
+        result.push_back(reduce_to_3x3(c));
+    }
+    return result;
+}
+
 } // namespace aare

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"

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"
@@ -105,7 +106,7 @@ class Frame {
      * @tparam T type of the pixels
      * @return NDView<T, 2>
      */
-    template <typename T> NDView<T, 2> view() {
+    template <typename T> NDView<T, 2> view() & {
         std::array<ssize_t, 2> shape = {static_cast<ssize_t>(m_rows),
                                         static_cast<ssize_t>(m_cols)};
         T *data = reinterpret_cast<T *>(m_data);

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"
@@ -17,7 +18,10 @@ struct Photon {
 };
 
 class Interpolator {
+    // marginal CDF of eta_x (if rosenblatt applied), conditional
+    // CDF of eta_x conditioned on eta_y
     NDArray<double, 3> m_ietax;
+    // conditional CDF of eta_y conditioned on eta_x
     NDArray<double, 3> m_ietay;
 
     NDArray<double, 1> m_etabinsx;
@@ -25,106 +29,210 @@ 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
+     * @param xbins bin edges for etaX
+     * @param ybins bin edges for etaY
+     * @param ebins bin edges for photon energy
+     * @note note first dimension is etaX, second etaY, third photon energy
+     */
     Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
                  NDView<double, 1> ybins, NDView<double, 1> ebins);
+
+    /**
+     * @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
+     * @note note 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)
+     */
+    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);
+
+  private:
+    /**
+     * @brief implements underlying interpolation logic based on EtaFunction
+     * Type
+     * @tparam EtaFunction Function object that calculates desired eta default:
+     * @param u: transformed photon position in x between [0,1]
+     * @param v: transformed photon position in y between [0,1]
+     * @param c: corner of eta
+     */
+    template <auto EtaFunction, typename ClusterType>
+    void interpolation_logic(Photon &photon, const double u, const double v,
+                             const corner c = corner::cTopLeft);
+
+    /**
+     *  @brief bilinear interpolation of the transformed eta values
+     * @param ix index of etaX bin
+     * @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);
 };
 
-// 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) {
+    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 <auto EtaFunction, typename ClusterType, typename Enable>
 std::vector<Photon>
 Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) {
     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);
+        // std::cout << "eta.x: " << eta.x << " eta.y: " << eta.y << std::endl;
 
-            // fmt::print("ex: {}, ix: {}, iy: {}\n", ie, ix, iy);
+        // 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);
 
-            double dX, dY;
-            // cBottomLeft = 0,
-            // cBottomRight = 1,
-            // cTopLeft = 2,
-            // cTopRight = 3
-            switch (static_cast<corner>(eta.c)) {
-            case corner::cTopLeft:
-                dX = -1.;
-                dY = 0;
-                break;
-            case corner::cTopRight:;
-                dX = 0;
-                dY = 0;
-                break;
-            case corner::cBottomLeft:
-                dX = -1.;
-                dY = -1.;
-                break;
-            case corner::cBottomRight:
-                dX = 0.;
-                dY = -1.;
-                break;
-            }
-            photon.x += m_ietax(ix, iy, ie) * 2 + dX;
-            photon.y += m_ietay(ix, iy, ie) * 2 + 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);
+        // std::cout << "ix: " << ix << " iy: " << iy << std::endl;
 
-            Photon photon;
-            photon.x = cluster.x;
-            photon.y = cluster.y;
-            photon.energy = static_cast<decltype(photon.energy)>(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);
 
-            // 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);
+        double ietax_interpolated = m_ietax(ix, iy, ie);
+        double ietay_interpolated = m_ietay(ix, iy, ie);
 
-            photon.x += m_ietax(ix, iy, ie) *
-                        2; // eta goes between 0 and 1 but we could move the hit
-                           // anywhere in the 2x2
-            photon.y += m_ietay(ix, iy, ie) * 2;
-            photons.push_back(photon);
-        }
+        interpolation_logic<EtaFunction, ClusterType>(
+            photon, ietax_interpolated, ietay_interpolated, eta.c);
 
-    } else {
-        throw std::runtime_error(
-            "Only 3x3 and 2x2 clusters are supported for interpolation");
+        photons.push_back(photon);
     }
 
     return photons;
 }
 
+template <auto EtaFunction, typename ClusterType>
+void Interpolator::interpolation_logic(Photon &photon, const double u,
+                                       const double v, const corner c) {
+
+    // std::cout << "u: " << u << " v: " << v << std::endl;
+
+    // TODO: try to call this with std::is_same_v and have it constexpr if
+    // possible
+    if (EtaFunction == &calculate_eta2<typename ClusterType::value_type,
+                                       ClusterType::cluster_size_x,
+                                       ClusterType::cluster_size_y,
+                                       typename ClusterType::coord_type> ||
+        EtaFunction == &calculate_full_eta2<typename ClusterType::value_type,
+                                            ClusterType::cluster_size_x,
+                                            ClusterType::cluster_size_y,
+                                            typename ClusterType::coord_type>) {
+        double dX{}, dY{};
+
+        // TODO: could also chaneg the sign of the eta calculation
+        switch (c) {
+        case corner::cTopLeft:
+            dX = -1.0;
+            dY = -1.0;
+            break;
+        case corner::cTopRight:;
+            dX = 0.0;
+            dY = -1.0;
+            break;
+        case corner::cBottomLeft:
+            dX = -1.0;
+            dY = 0.0;
+            break;
+        case corner::cBottomRight:
+            dX = 0.0;
+            dY = 0.0;
+            break;
+        }
+        photon.x = photon.x + 0.5 + u + dX; // use pixel center + 0.5
+        photon.y = photon.y + 0.5 + v +
+                   dY; // eta2 calculates the ratio between bottom and sum of
+                       // bottom and top  shift by 1 add eta value correctly
+    } else {
+        photon.x += u;
+        photon.y += v;
+    }
+}
+
 } // namespace aare

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,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 /*
 Container holding image data, or a time series of image data in contigious

include/aare/NDView.hpp

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

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>
@@ -42,14 +43,16 @@ class RawFileNameComponents {
 
 class ScanParameters {
     bool m_enabled = false;
-    std::string m_dac;
+    DACIndex m_dac{};
     int m_start = 0;
     int m_stop = 0;
     int m_step = 0;
-    // TODO! add settleTime, requires string to time conversion
+    int64_t m_settleTime = 0; // [ns]
 
   public:
     ScanParameters(const std::string &par);
+    ScanParameters(const bool enabled, const DACIndex dac, const int start,
+                   const int stop, const int step, const int64_t settleTime);
     ScanParameters() = default;
     ScanParameters(const ScanParameters &) = default;
     ScanParameters &operator=(const ScanParameters &) = default;
@@ -57,8 +60,9 @@ class ScanParameters {
     int start() const;
     int stop() const;
     int step() const;
-    const std::string &dac() const;
+    DACIndex dac() const;
     bool enabled() const;
+    int64_t settleTime() const;
     void increment_stop();
 };
 

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>

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,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 
 #include <aare/NDView.hpp>

include/aare/defs.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 
 #include "aare/Dtype.hpp"
@@ -215,6 +216,135 @@ enum class DetectorType {
     Unknown
 };
 
+/**
+ * @brief Enum class to define the Digital to Analog converter
+ * The values are the same as in slsDetectorPackage
+ */
+enum DACIndex {
+    DAC_0,
+    DAC_1,
+    DAC_2,
+    DAC_3,
+    DAC_4,
+    DAC_5,
+    DAC_6,
+    DAC_7,
+    DAC_8,
+    DAC_9,
+    DAC_10,
+    DAC_11,
+    DAC_12,
+    DAC_13,
+    DAC_14,
+    DAC_15,
+    DAC_16,
+    DAC_17,
+    VSVP,
+    VTRIM,
+    VRPREAMP,
+    VRSHAPER,
+    VSVN,
+    VTGSTV,
+    VCMP_LL,
+    VCMP_LR,
+    VCAL,
+    VCMP_RL,
+    RXB_RB,
+    RXB_LB,
+    VCMP_RR,
+    VCP,
+    VCN,
+    VISHAPER,
+    VTHRESHOLD,
+    IO_DELAY,
+    VREF_DS,
+    VOUT_CM,
+    VIN_CM,
+    VREF_COMP,
+    VB_COMP,
+    VDD_PROT,
+    VIN_COM,
+    VREF_PRECH,
+    VB_PIXBUF,
+    VB_DS,
+    VREF_H_ADC,
+    VB_COMP_FE,
+    VB_COMP_ADC,
+    VCOM_CDS,
+    VREF_RSTORE,
+    VB_OPA_1ST,
+    VREF_COMP_FE,
+    VCOM_ADC1,
+    VREF_L_ADC,
+    VREF_CDS,
+    VB_CS,
+    VB_OPA_FD,
+    VCOM_ADC2,
+    VCASSH,
+    VTH2,
+    VRSHAPER_N,
+    VIPRE_OUT,
+    VTH3,
+    VTH1,
+    VICIN,
+    VCAS,
+    VCAL_N,
+    VIPRE,
+    VCAL_P,
+    VDCSH,
+    VBP_COLBUF,
+    VB_SDA,
+    VCASC_SFP,
+    VIPRE_CDS,
+    IBIAS_SFP,
+    ADC_VPP,
+    HIGH_VOLTAGE,
+    TEMPERATURE_ADC,
+    TEMPERATURE_FPGA,
+    TEMPERATURE_FPGAEXT,
+    TEMPERATURE_10GE,
+    TEMPERATURE_DCDC,
+    TEMPERATURE_SODL,
+    TEMPERATURE_SODR,
+    TEMPERATURE_FPGA2,
+    TEMPERATURE_FPGA3,
+    TRIMBIT_SCAN,
+    V_POWER_A = 100,
+    V_POWER_B = 101,
+    V_POWER_C = 102,
+    V_POWER_D = 103,
+    V_POWER_IO = 104,
+    V_POWER_CHIP = 105,
+    I_POWER_A = 106,
+    I_POWER_B = 107,
+    I_POWER_C = 108,
+    I_POWER_D = 109,
+    I_POWER_IO = 110,
+    V_LIMIT = 111,
+    SLOW_ADC0 = 1000,
+    SLOW_ADC1,
+    SLOW_ADC2,
+    SLOW_ADC3,
+    SLOW_ADC4,
+    SLOW_ADC5,
+    SLOW_ADC6,
+    SLOW_ADC7,
+    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 };
 
@@ -231,6 +361,15 @@ 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
+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);
 
 } // namespace aare

include/aare/logger.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 /*Utility to log to console*/
 

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,34 +29,33 @@ def _get_class(name, cluster_size, dtype):
 
 
 
-def ClusterFinder(image_size, saved_cluster_size, checked_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++.
     """
-    cls = _get_class("ClusterFinder", saved_cluster_size, dtype)
-    return cls(image_size, n_sigma=n_sigma, capacity=capacity, cluster_size_x=checked_cluster_size[0], cluster_size_y=checked_cluster_size[1])
+    cls = _get_class("ClusterFinder", cluster_size, dtype)
+    return cls(image_size, n_sigma=n_sigma, capacity=capacity)
 
 
 
-def ClusterFinderMT(image_size, saved_cluster_size = (3,3), checked_cluster_size = (3,3), dtype=np.int32, n_sigma=5, capacity = 1024, n_threads = 3): 
+def ClusterFinderMT(image_size, cluster_size = (3,3), dtype=np.int32, n_sigma=5, capacity = 1024, n_threads = 3): 
     """ 
     Factory function to create a ClusterFinderMT object. Provides a cleaner syntax for 
     the templated ClusterFinderMT in C++.
     """
 
-    cls = _get_class("ClusterFinderMT", saved_cluster_size, dtype)
-    return cls(image_size, n_sigma=n_sigma, capacity=capacity, n_threads=n_threads, cluster_size_x=checked_cluster_size[0], cluster_size_y=checked_cluster_size[1])
+    cls = _get_class("ClusterFinderMT", cluster_size, dtype)
+    return cls(image_size, n_sigma=n_sigma, capacity=capacity, n_threads=n_threads)
 
 
-
-def ClusterCollector(clusterfindermt, cluster_size = (3,3), dtype=np.int32): 
+def ClusterCollector(clusterfindermt, dtype=np.int32): 
     """ 
     Factory function to create a ClusterCollector object. Provides a cleaner syntax for 
     the templated ClusterCollector in C++.
     """
-
-    cls = _get_class("ClusterCollector", cluster_size, dtype)
+    
+    cls = _get_class("ClusterCollector", clusterfindermt.cluster_size, dtype)
     return cls(clusterfindermt)
 
 def ClusterFileSink(clusterfindermt, cluster_file, dtype=np.int32): 
@@ -66,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++.
@@ -84,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,17 +8,19 @@ 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
 

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>
@@ -24,7 +25,8 @@ void define_Cluster(py::module &m, const std::string &typestr) {
     py::class_<Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>>(
         m, class_name.c_str(), py::buffer_protocol())
 
-        .def(py::init([](uint8_t x, uint8_t y, py::array_t<Type> data) {
+        .def(py::init([](CoordType x, CoordType y,
+                         py::array_t<Type, py::array::forcecast> data) {
             py::buffer_info buf_info = data.request();
             Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType> cluster;
             cluster.x = x;
@@ -34,31 +36,73 @@ void define_Cluster(py::module &m, const std::string &typestr) {
                 cluster.data[i] = r(i);
             }
             return cluster;
-        }));
+        }))
 
-    /*
-    //TODO! Review if to keep or not
-    .def_property(
-        "data",
-        [](ClusterType &c) -> py::array {
-            return py::array(py::buffer_info(
-                c.data, sizeof(Type),
-                py::format_descriptor<Type>::format(), // Type
-                                                       // format
-                1, // Number of dimensions
-                {static_cast<ssize_t>(ClusterSizeX *
-                                      ClusterSizeY)}, // Shape (flattened)
-                {sizeof(Type)} // Stride (step size between elements)
-                ));
+        // TODO! Review if to keep or not
+        .def_property_readonly(
+            "data",
+            [](Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType> &c)
+                -> py::array {
+                return py::array(py::buffer_info(
+                    c.data.data(), sizeof(Type),
+                    py::format_descriptor<Type>::format(), // Type
+                                                           // format
+                    2, // Number of dimensions
+                    {static_cast<ssize_t>(ClusterSizeX),
+                     static_cast<ssize_t>(ClusterSizeY)}, // Shape (flattened)
+                    {sizeof(Type) * ClusterSizeY, sizeof(Type)}
+                    // Stride (step size between elements)
+                    ));
+            })
+
+        .def_readonly("x",
+                      &Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>::x)
+
+        .def_readonly("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,
+          typename CoordType = int16_t>
+void reduce_to_3x3(py::module &m) {
+
+    m.def(
+        "reduce_to_3x3",
+        [](const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
+            return reduce_to_3x3(cl);
         },
-        [](ClusterType &c, py::array_t<Type> arr) {
-            py::buffer_info buf_info = arr.request();
-            Type *ptr = static_cast<Type *>(buf_info.ptr);
-            std::copy(ptr, ptr + ClusterSizeX * ClusterSizeY,
-                      c.data); // TODO dont iterate over centers!!!
+        py::return_value_policy::move, R"(Reduce cluster to 3x3 subcluster)");
+}
 
-        });
-    */
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = int16_t>
+void reduce_to_2x2(py::module &m) {
+
+    m.def(
+        "reduce_to_2x2",
+        [](const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
+            return reduce_to_2x2(cl);
+        },
+        py::return_value_policy::move,
+        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"
@@ -44,14 +45,15 @@ void define_ClusterFile(py::module &m, const std::string &typestr) {
                  auto v = new ClusterVector<ClusterType>(self.read_frame());
                  return v;
              })
-        .def("set_roi", &ClusterFile<ClusterType>::set_roi,
-             py::arg("roi"))
+        .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) {
                 auto view = make_view_2d(noise_map);
                 self.set_noise_map(view);
-            }, py::arg("noise_map"))
+            },
+            py::arg("noise_map"))
 
         .def("set_gain_map",
              [](ClusterFile<ClusterType> &self, py::array_t<double> gain_map) {
@@ -80,15 +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);
-          });
-}
-
 #pragma GCC diagnostic pop

python/src/bind_ClusterFileSink.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_ClusterFinder.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #include "aare/ClusterCollector.hpp"
 #include "aare/ClusterFileSink.hpp"
 #include "aare/ClusterFinder.hpp"
@@ -30,9 +31,8 @@ void define_ClusterFinder(py::module &m, const std::string &typestr) {
 
     py::class_<ClusterFinder<ClusterType, uint16_t, pd_type>>(
         m, class_name.c_str())
-        .def(py::init<Shape<2>, pd_type, size_t, uint32_t, uint32_t>(), py::arg("image_size"),
-             py::arg("n_sigma") = 5.0, py::arg("capacity") = 1'000'000,
-             py::arg("cluster_size_x") = 3, py::arg("cluster_size_y") = 3)
+        .def(py::init<Shape<2>, pd_type, size_t>(), py::arg("image_size"),
+             py::arg("n_sigma") = 5.0, py::arg("capacity") = 1'000'000)
         .def("push_pedestal_frame",
              [](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
                 py::array_t<uint16_t> frame) {

python/src/bind_ClusterFinderMT.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #include "aare/ClusterCollector.hpp"
 #include "aare/ClusterFileSink.hpp"
 #include "aare/ClusterFinder.hpp"
@@ -30,10 +31,9 @@ void define_ClusterFinderMT(py::module &m, const std::string &typestr) {
 
     py::class_<ClusterFinderMT<ClusterType, uint16_t, pd_type>>(
         m, class_name.c_str())
-        .def(py::init<Shape<2>, pd_type, size_t, size_t, uint32_t, uint32_t>(),
+        .def(py::init<Shape<2>, pd_type, size_t, size_t>(),
              py::arg("image_size"), py::arg("n_sigma") = 5.0,
-             py::arg("capacity") = 2048, py::arg("n_threads") = 3,
-             py::arg("cluster_size_x") = 3, py::arg("cluster_size_y") = 3)
+             py::arg("capacity") = 2048, py::arg("n_threads") = 3)
         .def("push_pedestal_frame",
              [](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
                 py::array_t<uint16_t> frame) {

python/src/bind_ClusterVector.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #include "aare/ClusterCollector.hpp"
 #include "aare/ClusterFileSink.hpp"
 #include "aare/ClusterFinder.hpp"
@@ -44,11 +45,16 @@ void define_ClusterVector(py::module &m, const std::string &typestr) {
                  auto *vec = new std::vector<Type>(self.sum());
                  return return_vector(vec);
              })
-        .def("sum_2x2",
-             [](ClusterVector<ClusterType> &self) {
-                 auto *vec = new std::vector<Type>(self.sum_2x2());
-                 return return_vector(vec);
-             })
+        .def(
+            "sum_2x2",
+            [](ClusterVector<ClusterType> &self) {
+                auto *vec = new std::vector<Sum_index_pair<Type, corner>>(
+                    self.sum_2x2());
+
+                return return_vector(vec);
+            },
+            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
+            )")
         .def_property_readonly("size", &ClusterVector<ClusterType>::size)
         .def("item_size", &ClusterVector<ClusterType>::item_size)
         .def_property_readonly("fmt",
@@ -104,4 +110,47 @@ void define_ClusterVector(py::module &m, const std::string &typestr) {
           });
 }
 
+template <typename Type, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = uint16_t>
+void define_2x2_reduction(py::module &m) {
+    m.def(
+        "reduce_to_2x2",
+        [](const ClusterVector<
+            Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>> &cv) {
+            return new ClusterVector<Cluster<Type, 2, 2, CoordType>>(
+                reduce_to_2x2(cv));
+        },
+        R"(
+        Reduce cluster to 2x2 subcluster by taking the 2x2 subcluster with 
+        the highest photon energy.
+
+        Parameters
+
+        cv : ClusterVector (clusters are 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.)
+
+        )",
+        py::arg("clustervector"));
+}
+
+template <typename Type, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = uint16_t>
+void define_3x3_reduction(py::module &m) {
+
+    m.def(
+        "reduce_to_3x3",
+        [](const ClusterVector<
+            Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>> &cv) {
+            return new ClusterVector<Cluster<Type, 3, 3, CoordType>>(
+                reduce_to_3x3(cv));
+        },
+        R"(
+        Reduce cluster to 3x3 subcluster 
+
+        Parameters
+        
+        cv : ClusterVector   
+        )",
+        py::arg("clustervector"));
+}
+
 #pragma GCC diagnostic pop

python/src/bind_Eta.hpp

@@ -0,0 +1,104 @@
+#include "aare/CalculateEta.hpp"
+
+#include <cstdint>
+// #include <pybind11/native_enum.h> only for version 3
+#include <pybind11/pybind11.h>
+
+namespace py = pybind11;
+using namespace ::aare;
+
+template <typename T>
+void define_eta(py::module &m, const std::string &typestr) {
+    auto class_name = fmt::format("Eta{}", typestr);
+
+    py::class_<Eta2<T>>(m, class_name.c_str())
+        .def(py::init<>())
+        .def_readonly("x", &Eta2<T>::x, "eta x value")
+        .def_readonly("y", &Eta2<T>::y, "eta y value")
+        .def_readonly("c", &Eta2<T>::c,
+                      "eta corner value cTopLeft, cTopRight, "
+                      "cBottomLeft, cBottomRight")
+        .def_readonly("sum", &Eta2<T>::sum, "photon energy of cluster");
+}
+
+void define_corner_enum(py::module &m) {
+    py::enum_<corner>(m, "corner", "enum.Enum")
+        .value("cTopLeft", corner::cTopLeft)
+        .value("cTopRight", corner::cTopRight)
+        .value("cBottomLeft", corner::cBottomLeft)
+        .value("cBottomRight", corner::cBottomRight)
+        .export_values();
+}
+
+template <typename Type, uint8_t CoordSizeX, uint8_t CoordSizeY,
+          typename CoordType = uint16_t>
+void register_calculate_2x2eta(py::module &m) {
+    using ClusterType = Cluster<Type, CoordSizeX, CoordSizeY, CoordType>;
+
+    m.def(
+        "calculate_eta2",
+        [](const aare::ClusterVector<ClusterType> &clusters) {
+            auto eta2 = new std::vector<Eta2<typename ClusterType::value_type>>(
+                calculate_eta2(clusters));
+            return return_vector(eta2);
+        },
+        R"(calculates eta2x2)", py::arg("clusters"));
+
+    m.def(
+        "calculate_eta2",
+        [](const aare::Cluster<Type, CoordSizeX, CoordSizeY, CoordType>
+               &cluster) { return calculate_eta2(cluster); },
+        R"(calculates eta2x2)", py::arg("cluster"));
+
+    m.def(
+        "calculate_full_eta2",
+        [](const aare::Cluster<Type, CoordSizeX, CoordSizeY, CoordType>
+               &cluster) { return calculate_full_eta2(cluster); },
+        R"(calculates full eta2x2)", py::arg("cluster"));
+
+    m.def(
+        "calculate_full_eta2",
+        [](const aare::ClusterVector<ClusterType> &clusters) {
+            auto eta2 = new std::vector<Eta2<typename ClusterType::value_type>>(
+                calculate_full_eta2(clusters));
+            return return_vector(eta2);
+        },
+        R"(calculates full eta2x2)", py::arg("clusters"));
+}
+
+template <typename Type, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = uint16_t>
+void register_calculate_3x3eta(py::module &m) {
+    using ClusterType = Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>;
+
+    m.def(
+        "calculate_eta3",
+        [](const aare::ClusterVector<ClusterType> &clusters) {
+            auto eta = new std::vector<Eta2<Type>>(calculate_eta3(clusters));
+            return return_vector(eta);
+        },
+        R"(calculates eta3x3 using entire cluster)", py::arg("clusters"));
+
+    m.def(
+        "calculate_cross_eta3",
+        [](const aare::ClusterVector<ClusterType> &clusters) {
+            auto eta =
+                new std::vector<Eta2<Type>>(calculate_cross_eta3(clusters));
+            return return_vector(eta);
+        },
+        R"(calculates eta3x3 taking into account cross pixels in cluster)",
+        py::arg("clusters"));
+
+    m.def(
+        "calculate_eta3",
+        [](const ClusterType &cluster) { return calculate_eta3(cluster); },
+        R"(calculates eta3x3 using entire cluster)", py::arg("cluster"));
+
+    m.def(
+        "calculate_cross_eta3",
+        [](const ClusterType &cluster) {
+            return calculate_cross_eta3(cluster);
+        },
+        R"(calculates eta3x3 taking into account cross pixels in cluster)",
+        py::arg("cluster"));
+}

python/src/bind_calibration.hpp

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

python/src/ctb_raw_file.hpp

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

python/src/file.hpp

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

python/src/fit.hpp

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

python/src/interpolation.hpp

@@ -1,3 +1,5 @@
+// SPDX-License-Identifier: MPL-2.0
+#include "aare/CalculateEta.hpp"
 #include "aare/Interpolator.hpp"
 #include "aare/NDArray.hpp"
 #include "aare/NDView.hpp"
@@ -9,19 +11,41 @@
 
 namespace py = pybind11;
 
+#define REGISTER_INTERPOLATOR_ETA2(T, N, M, U)                                 \
+    register_interpolate<T, N, M, U, aare::calculate_full_eta2<T, N, M, U>>(   \
+        interpolator, "_full_eta2", "full eta2");                              \
+    register_interpolate<T, N, M, U, aare::calculate_eta2<T, N, M, U>>(        \
+        interpolator, "", "eta2");
+
+#define REGISTER_INTERPOLATOR_ETA3(T, N, M, U)                                 \
+    register_interpolate<T, N, M, U, aare::calculate_eta3<T, N, M, U>>(        \
+        interpolator, "_eta3", "full eta3");                                   \
+    register_interpolate<T, N, M, U, aare::calculate_cross_eta3<T, N, M, U>>(  \
+        interpolator, "_cross_eta3", "cross eta3");
+
 template <typename Type, uint8_t CoordSizeX, uint8_t CoordSizeY,
-          typename CoordType = uint16_t>
-void register_interpolate(py::class_<aare::Interpolator> &interpolator) {
+          typename CoordType = uint16_t, auto EtaFunction>
+void register_interpolate(py::class_<aare::Interpolator> &interpolator,
+                          const std::string &typestr = "",
+                          const std::string &doc_string_etatype = "eta2x2") {
 
     using ClusterType = Cluster<Type, CoordSizeX, CoordSizeY, CoordType>;
 
-    interpolator.def("interpolate",
-                     [](aare::Interpolator &self,
-                        const ClusterVector<ClusterType> &clusters) {
-                         auto photons = self.interpolate<ClusterType>(clusters);
-                         auto *ptr = new std::vector<Photon>{photons};
-                         return return_vector(ptr);
-                     });
+    const std::string docstring = "interpolation based on " +
+                                  doc_string_etatype +
+                                  "\n\nReturns:\n interpolated photons";
+
+    auto function_name = fmt::format("interpolate{}", typestr);
+
+    interpolator.def(
+        function_name.c_str(),
+        [](aare::Interpolator &self,
+           const ClusterVector<ClusterType> &clusters) {
+            auto photons = self.interpolate<EtaFunction, ClusterType>(clusters);
+            auto *ptr = new std::vector<Photon>{photons};
+            return return_vector(ptr);
+        },
+        docstring.c_str(), py::arg("cluster_vector"));
 }
 
 void define_interpolation_bindings(py::module &m) {
@@ -30,33 +54,91 @@ void define_interpolation_bindings(py::module &m) {
 
     auto interpolator =
         py::class_<aare::Interpolator>(m, "Interpolator")
-            .def(py::init([](py::array_t<double, py::array::c_style |
-                                                     py::array::forcecast>
-                                 etacube,
-                             py::array_t<double> xbins,
-                             py::array_t<double> ybins,
-                             py::array_t<double> ebins) {
-                return Interpolator(make_view_3d(etacube), make_view_1d(xbins),
-                                    make_view_1d(ybins), make_view_1d(ebins));
-            }))
+            .def(py::init(
+                [](py::array_t<double,
+                               py::array::c_style | py::array::forcecast>
+                       etacube,
+                   py::array_t<double> xbins, py::array_t<double> ybins,
+                   py::array_t<double> ebins) {
+                    return Interpolator(
+                        make_view_3d(etacube), make_view_1d(xbins),
+                        make_view_1d(ybins), make_view_1d(ebins));
+                }), 
+                R"doc(
+                Constructor 
+
+                Args:
+                
+                etacube: 
+                    joint distribution of eta_x, eta_y and photon energy (**Note:** for the joint distribution first dimension is eta_x, second: eta_y, third: energy bins.)
+                xbins: 
+                    bin edges of etax
+                ybins: 
+                    bin edges of etay
+                ebins: 
+                    bin edges of photon energy
+                )doc",
+                py::arg("etacube"),
+                py::arg("xbins"), py::arg("ybins"),
+                py::arg("ebins"))
+
+            .def(py::init(
+                [](py::array_t<double> xbins, py::array_t<double> ybins,
+                   py::array_t<double> ebins) {
+                    return Interpolator(make_view_1d(xbins),
+                                        make_view_1d(ybins),
+                                        make_view_1d(ebins));
+                }),
+                R"(
+                Constructor
+
+                Args: 
+                
+                xbins: 
+                    bin edges of etax
+                ybins: 
+                    bin edges of etay
+                ebins: 
+                    bin edges of photon energy
+                )", py::arg("xbins"),
+                py::arg("ybins"), py::arg("ebins"))
+            .def(
+                "rosenblatttransform",
+                [](Interpolator &self,
+                   py::array_t<double,
+                               py::array::c_style | py::array::forcecast>
+                       etacube) {
+                    return self.rosenblatttransform(make_view_3d(etacube));
+                },
+                R"(
+                calculated the rosenblatttransform for the given distribution
+                
+                etacube: 
+                    joint distribution of eta_x, eta_y and photon energy (**Note:** for the joint distribution first dimension is eta_x, second: eta_y, third: energy bins.)
+                )",
+                py::arg("etacube"))
             .def("get_ietax",
                  [](Interpolator &self) {
                      auto *ptr = new NDArray<double, 3>{};
                      *ptr = self.get_ietax();
                      return return_image_data(ptr);
-                 })
+                 }, R"(conditional CDF of etax conditioned on etay, marginal CDF of etax (if rosenblatt transform applied))")
             .def("get_ietay", [](Interpolator &self) {
                 auto *ptr = new NDArray<double, 3>{};
                 *ptr = self.get_ietay();
                 return return_image_data(ptr);
-            });
+            }, R"(conditional CDF of etay conditioned on etax)");
 
-    register_interpolate<int, 3, 3, uint16_t>(interpolator);
-    register_interpolate<float, 3, 3, uint16_t>(interpolator);
-    register_interpolate<double, 3, 3, uint16_t>(interpolator);
-    register_interpolate<int, 2, 2, uint16_t>(interpolator);
-    register_interpolate<float, 2, 2, uint16_t>(interpolator);
-    register_interpolate<double, 2, 2, uint16_t>(interpolator);
+    REGISTER_INTERPOLATOR_ETA3(int, 3, 3, uint16_t);
+    REGISTER_INTERPOLATOR_ETA3(float, 3, 3, uint16_t);
+    REGISTER_INTERPOLATOR_ETA3(double, 3, 3, uint16_t);
+
+    REGISTER_INTERPOLATOR_ETA2(int, 3, 3, uint16_t);
+    REGISTER_INTERPOLATOR_ETA2(float, 3, 3, uint16_t);
+    REGISTER_INTERPOLATOR_ETA2(double, 3, 3, uint16_t);
+    REGISTER_INTERPOLATOR_ETA2(int, 2, 2, uint16_t);
+    REGISTER_INTERPOLATOR_ETA2(float, 2, 2, uint16_t);
+    REGISTER_INTERPOLATOR_ETA2(double, 2, 2, uint16_t);
 
     // TODO! Evaluate without converting to double
     m.def(

python/src/jungfrau_data_file.hpp

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

python/src/module.cpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 // Files with bindings to the different classes
 
 // New style file naming
@@ -8,6 +9,7 @@
 #include "bind_ClusterFinder.hpp"
 #include "bind_ClusterFinderMT.hpp"
 #include "bind_ClusterVector.hpp"
+#include "bind_Eta.hpp"
 #include "bind_calibration.hpp"
 
 // TODO! migrate the other names
@@ -42,12 +44,16 @@ double, 'f' for float)
 #define DEFINE_CLUSTER_BINDINGS(T, N, M, U, TYPE_CODE)                         \
     define_ClusterFile<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE);         \
     define_ClusterVector<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE);       \
+    define_Cluster<T, N, M, U>(m, #N "x" #M #TYPE_CODE);                       \
+    register_calculate_2x2eta<T, N, M, U>(m);                                  \
+    define_2x2_reduction<T, N, M, U>(m);                                       \
+    reduce_to_2x2<T, N, M, U>(m);
+
+#define DEFINE_BINDINGS_CLUSTERFINDER(T, N, M, U, TYPE_CODE)                   \
     define_ClusterFinder<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE);       \
     define_ClusterFinderMT<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE);     \
     define_ClusterFileSink<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE);     \
-    define_ClusterCollector<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE);    \
-    define_Cluster<T, N, M, U>(m, #N "x" #M #TYPE_CODE);                       \
-    register_calculate_eta<T, N, M, U>(m);
+    define_ClusterCollector<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE);
 
 PYBIND11_MODULE(_aare, m) {
     define_file_io_bindings(m);
@@ -85,8 +91,74 @@ PYBIND11_MODULE(_aare, m) {
     DEFINE_CLUSTER_BINDINGS(double, 9, 9, uint16_t, d);
     DEFINE_CLUSTER_BINDINGS(float, 9, 9, uint16_t, f);
 
-    DEFINE_CLUSTER_BINDINGS(int, 21, 21, uint16_t, i);
-    DEFINE_CLUSTER_BINDINGS(double, 21, 21, uint16_t, d);
-    DEFINE_CLUSTER_BINDINGS(float, 21, 21, uint16_t, f);
-    
+    DEFINE_CLUSTER_BINDINGS(int16_t, 3, 3, uint16_t, i16);
+
+    DEFINE_BINDINGS_CLUSTERFINDER(int, 3, 3, uint16_t, i);
+    DEFINE_BINDINGS_CLUSTERFINDER(double, 3, 3, uint16_t, d);
+    DEFINE_BINDINGS_CLUSTERFINDER(float, 3, 3, uint16_t, f);
+
+    DEFINE_BINDINGS_CLUSTERFINDER(int, 5, 5, uint16_t, i);
+    DEFINE_BINDINGS_CLUSTERFINDER(double, 5, 5, uint16_t, d);
+    DEFINE_BINDINGS_CLUSTERFINDER(float, 5, 5, uint16_t, f);
+
+    DEFINE_BINDINGS_CLUSTERFINDER(int, 7, 7, uint16_t, i);
+    DEFINE_BINDINGS_CLUSTERFINDER(double, 7, 7, uint16_t, d);
+    DEFINE_BINDINGS_CLUSTERFINDER(float, 7, 7, uint16_t, f);
+
+    DEFINE_BINDINGS_CLUSTERFINDER(int, 9, 9, uint16_t, i);
+    DEFINE_BINDINGS_CLUSTERFINDER(double, 9, 9, uint16_t, d);
+    DEFINE_BINDINGS_CLUSTERFINDER(float, 9, 9, uint16_t, f);
+
+    define_3x3_reduction<int, 3, 3, uint16_t>(m);
+    define_3x3_reduction<double, 3, 3, uint16_t>(m);
+    define_3x3_reduction<float, 3, 3, uint16_t>(m);
+    define_3x3_reduction<int, 5, 5, uint16_t>(m);
+    define_3x3_reduction<double, 5, 5, uint16_t>(m);
+    define_3x3_reduction<float, 5, 5, uint16_t>(m);
+    define_3x3_reduction<int, 7, 7, uint16_t>(m);
+    define_3x3_reduction<double, 7, 7, uint16_t>(m);
+    define_3x3_reduction<float, 7, 7, uint16_t>(m);
+    define_3x3_reduction<int, 9, 9, uint16_t>(m);
+    define_3x3_reduction<double, 9, 9, uint16_t>(m);
+    define_3x3_reduction<float, 9, 9, uint16_t>(m);
+
+    reduce_to_3x3<int, 3, 3, uint16_t>(m);
+    reduce_to_3x3<double, 3, 3, uint16_t>(m);
+    reduce_to_3x3<float, 3, 3, uint16_t>(m);
+    reduce_to_3x3<int, 5, 5, uint16_t>(m);
+    reduce_to_3x3<double, 5, 5, uint16_t>(m);
+    reduce_to_3x3<float, 5, 5, uint16_t>(m);
+    reduce_to_3x3<int, 7, 7, uint16_t>(m);
+    reduce_to_3x3<double, 7, 7, uint16_t>(m);
+    reduce_to_3x3<float, 7, 7, uint16_t>(m);
+    reduce_to_3x3<int, 9, 9, uint16_t>(m);
+    reduce_to_3x3<double, 9, 9, uint16_t>(m);
+    reduce_to_3x3<float, 9, 9, uint16_t>(m);
+
+    register_calculate_3x3eta<int, 3, 3, uint16_t>(m);
+    register_calculate_3x3eta<double, 3, 3, uint16_t>(m);
+    register_calculate_3x3eta<float, 3, 3, uint16_t>(m);
+    register_calculate_3x3eta<int16_t, 3, 3, uint16_t>(m);
+
+    using Sum_index_pair_d = Sum_index_pair<double, corner>;
+    PYBIND11_NUMPY_DTYPE(Sum_index_pair_d, sum, index);
+    using Sum_index_pair_f = Sum_index_pair<float, corner>;
+    PYBIND11_NUMPY_DTYPE(Sum_index_pair_f, sum, index);
+    using Sum_index_pair_i = Sum_index_pair<int, corner>;
+    PYBIND11_NUMPY_DTYPE(Sum_index_pair_i, sum, index);
+
+    using eta_d = Eta2<double>;
+    PYBIND11_NUMPY_DTYPE(eta_d, x, y, c, sum);
+    using eta_i = Eta2<int>;
+    PYBIND11_NUMPY_DTYPE(eta_i, x, y, c, sum);
+    using eta_f = Eta2<float>;
+    PYBIND11_NUMPY_DTYPE(eta_f, x, y, c, sum);
+    using eta_i16 = Eta2<int16_t>;
+    PYBIND11_NUMPY_DTYPE(eta_i16, x, y, c, sum);
+
+    define_corner_enum(m);
+    define_eta<float>(m, "f");
+    define_eta<double>(m, "d");
+    define_eta<int>(m, "i");
+    define_eta<int16_t>(m, "i16");
 }