Compare commits

..

35 Commits

Author SHA1 Message Date
Erik Fröjdh
80a2b02345 Dev/decode my302 (#254)
Some checks failed
Build on RHEL8 / build (push) Failing after 0s
Build on RHEL9 / build (push) Failing after 0s
This PR adds support for decoding digital data from the my320 test chip.
- Added BitOffset (strong type)
- Expand 24 to 32 bit 
- Python bindings for decoding my302
- Improved docs
2025-12-09 18:27:02 +01:00
e795310b16 fixed tests (#252)
All checks were successful
Build on RHEL8 / build (push) Successful in 3m11s
Build on RHEL9 / build (push) Successful in 3m46s
- fixed failed tests 
- removed import of pickle, scipy 
- still requires boost_histogram, pytest_check

Co-authored-by: Erik Fröjdh <erik.frojdh@psi.ch>
2025-11-28 11:28:13 +01:00
1d8b68bf75 Run catch2 tests in CI (#253)
Some checks failed
Build on RHEL9 / build (push) Failing after 0s
Build on RHEL8 / build (push) Successful in 3m28s
- Build and run tests in CI
- Added macOS builds (and tests)
- Renamed workflow to build_with_docs.yml
2025-11-27 08:58:24 +01:00
Erik Fröjdh
cf53922bd3 restrict upload artifact 2025-11-26 20:21:57 +01:00
Erik Fröjdh
13cffb1ea8 WIP 2025-11-26 20:11:38 +01:00
Erik Fröjdh
bea373a112 with file 2025-11-26 20:09:05 +01:00
Erik Fröjdh
a66ce15a6c renamed workflow and added macos-latest 2025-11-26 20:07:32 +01:00
Erik Fröjdh
44fd015cfe added catch2 to dev env 2025-11-26 20:00:03 +01:00
Erik Fröjdh
6fe822d5dd WIP 2025-11-26 19:56:53 +01:00
Erik Fröjdh
8ed874a679 using tests 2025-11-26 19:54:43 +01:00
dd5ed138cf Dev/print filepath in error (#251)
Some checks failed
Build on RHEL8 / build (push) Failing after 0s
Build on RHEL9 / build (push) Successful in 3m15s
2025-11-25 11:25:44 +01:00
8201c5e999 Fix/mythenfilereading (#250)
Some checks failed
Build on RHEL8 / build (push) Failing after 0s
Build on RHEL9 / build (push) Successful in 3m23s
Add counter mask as member of RawFile
BugFix: temporarily handle -1 for ROI in mythenfile
2025-11-24 12:29:08 +01:00
03af5927ad Release 2025.11.21 (#249)
All checks were successful
Build on RHEL9 / build (push) Successful in 3m19s
Build on RHEL8 / build (push) Successful in 3m27s
- Updated VERSION and script 
- Updated release notes
2025-11-21 16:34:57 +01:00
Erik Fröjdh
452cfcb60f updated release notes 2025-11-21 16:17:04 +01:00
Erik Fröjdh
e8402d9d36 updated version and version script 2025-11-21 16:14:05 +01:00
a9de336817 added license (#247)
- Added LICENSE file
- Added SPX identifier to source files
2025-11-21 15:11:30 +01:00
6f7cb4ae30 Merge branch 'main' into dev/license 2025-11-21 14:52:54 +01:00
267ca87ab0 Dev/rosenblatttransform (#241)
- added rosenblatttransform 
- added 3x3 eta methods 
- interpolation can be used with various eta functions
- added documentation for interpolation, eta calculation 
- exposed full eta struct in python 
- disable ClusterFinder for 2x2 clusters 
- factory function for ClusterVector

---------

Co-authored-by: Dhanya Thattil <dhanya.thattil@psi.ch>
Co-authored-by: Erik Fröjdh <erik.frojdh@psi.ch>
2025-11-21 14:48:46 +01:00
Erik Fröjdh
5dbb969bcc Merge branch 'dev/license' of github.com:slsdetectorgroup/aare into dev/license 2025-11-21 14:45:12 +01:00
Erik Fröjdh
d58d8ea82f added comment in readme 2025-11-21 14:44:54 +01:00
f61f76ccf7 changed License in update_version.py added to etc
Some checks failed
Build on RHEL9 / build (push) Failing after 3m26s
Build on RHEL8 / build (push) Failing after 3m36s
2025-11-21 10:29:12 +01:00
Erik Fröjdh
200ae91622 also hpp 2025-11-21 10:14:14 +01:00
Erik Fröjdh
53aed8d8c6 added license 2025-11-20 09:01:28 +01:00
7fb500c44c Dev/expose sum 2x2 to python (#238)
Some checks failed
Build on RHEL8 / build (push) Failing after 3m15s
Build on RHEL9 / build (push) Failing after 3m16s
Saverio requested that max_sum_2x2 exposes index information in  python 
- max_sum_2x2 returns a corner as index
- replaced eta corner with corner enum class
- max_sum_2x2 now returns index as well in python 
- added link to Documenation in README

Note: Some Tests fail in EtaCalculation due to previous PR about
updating Eta 2x2 will fix in other PR
2025-10-27 20:04:16 +01:00
8989d2eb4a Merge branch 'main' into dev/expose_sum_2x2_to_python 2025-10-27 19:47:09 +01:00
c8c681faa8 updated release notes 2025-10-27 19:30:43 +01:00
0faaf2bbc7 updated release notes
Some checks failed
Build on RHEL8 / build (push) Failing after 3m9s
Build on RHEL9 / build (push) Failing after 3m16s
2025-10-27 18:45:23 +01:00
Erik Fröjdh
ac83eeff9b added tell and better error checking to cluster file (#239)
Some checks failed
Build on RHEL8 / build (push) Failing after 3m8s
Build on RHEL9 / build (push) Failing after 3m18s
- Check feof and ferror when reading frames
- added tell member function to ClusterFile
2025-10-27 15:46:31 +01:00
df7b9be5a5 added docstrings wrap struct into tuple
Some checks failed
Build on RHEL8 / build (push) Failing after 3m42s
Build on RHEL9 / build (push) Failing after 3m41s
2025-10-23 19:16:33 +02:00
dbffea15c0 fix: included deleted file 2025-10-23 17:50:17 +02:00
6e38c3259b added documentation link in README 2025-10-23 17:40:55 +02:00
73e8fd31c9 vector class no longer needed 2025-10-23 17:36:29 +02:00
b28abb2668 updated tests 2025-10-23 17:35:16 +02:00
01fa61cf47 index now returns enum type 2025-10-23 17:34:54 +02:00
790dd63ba3 make max_sum_2x2 properly accessible from python 2025-10-23 15:00:52 +02:00
179 changed files with 4465 additions and 1366 deletions

View File

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

View File

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

View File

@@ -1,3 +1,4 @@
# SPDX-License-Identifier: MPL-2.0
cmake_minimum_required(VERSION 3.15)
project(aare
@@ -442,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

373
LICENSE Normal file
View File

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

View File

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

View File

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

View File

@@ -1 +1 @@
2025.8.22
2025.11.21

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

BIN
docs/figures/Eta2x2.pdf Normal file

Binary file not shown.

BIN
docs/figures/Eta2x2.png Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 6.7 KiB

BIN
docs/figures/Eta2x2Full.pdf Normal file

Binary file not shown.

BIN
docs/figures/Eta2x2Full.png Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 10 KiB

BIN
docs/figures/Eta3x3.pdf Normal file

Binary file not shown.

BIN
docs/figures/Eta3x3.png Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 13 KiB

Binary file not shown.

Binary file not shown.

After

Width:  |  Height:  |  Size: 9.5 KiB

15
docs/src/Cluster.rst Normal file
View File

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

102
docs/src/Interpolation.rst Normal file
View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

103
etc/add_license.py Normal file
View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

@@ -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,108 +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
// TODO: could also chaneg the sign of the eta calculation
switch (static_cast<corner>(eta.c)) {
case corner::cTopLeft:
dX = 0.0;
dY = 0.0;
break;
case corner::cTopRight:;
dX = 1.0;
dY = 0.0;
break;
case corner::cBottomLeft:
dX = 0.0;
dY = 1.0;
break;
case corner::cBottomRight:
dX = 1.0;
dY = 1.0;
break;
}
photon.x -= m_ietax(ix, iy, ie) - dX;
photon.y -= m_ietay(ix, iy, ie) - dY;
photons.push_back(photon);
}
} else if (clusters.cluster_size_x() == 2 ||
clusters.cluster_size_y() == 2) {
for (const ClusterType &cluster : clusters) {
auto eta = calculate_eta2(cluster);
// 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);
// TODO: why 2?
photon.x -=
m_ietax(ix, iy, ie); // eta goes between 0 and 1 but we could
// move the hit anywhere in the 2x2
photon.y -= m_ietay(ix, iy, ie);
photons.push_back(photon);
}
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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/defs.hpp"
#include <algorithm>
@@ -102,6 +103,7 @@ class RawMasterFile {
std::optional<size_t> m_digital_samples;
std::optional<size_t> m_transceiver_samples;
std::optional<size_t> m_number_of_rows;
std::optional<uint8_t> m_counter_mask;
std::optional<ROI> m_roi;
@@ -132,6 +134,7 @@ class RawMasterFile {
std::optional<size_t> digital_samples() const;
std::optional<size_t> transceiver_samples() const;
std::optional<size_t> number_of_rows() const;
std::optional<uint8_t> counter_mask() const;
std::optional<ROI> roi() const;

View File

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

View File

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

View File

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

View File

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

View File

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

View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/Dtype.hpp"
@@ -331,6 +332,19 @@ enum DACIndex {
SLOW_ADC_TEMP
};
// helper pair class to easily expose in python
template <typename T1, typename T2> struct Sum_index_pair {
T1 sum;
T2 index;
};
enum class corner : int {
cTopLeft = 0,
cTopRight = 1,
cBottomLeft = 2,
cBottomRight = 3
};
enum class TimingMode { Auto, Trigger };
enum class FrameDiscardPolicy { NoDiscard, Discard, DiscardPartial };
@@ -358,4 +372,15 @@ constexpr uint16_t ADC_MASK =
*/
template <> DACIndex StringTo(const std::string &arg);
class BitOffset{
uint8_t m_offset{};
public:
BitOffset() = default;
explicit BitOffset(uint32_t offset);
uint8_t value() const {return m_offset;}
bool operator==(const BitOffset& other) const;
bool operator<(const BitOffset& other) const;
};
} // namespace aare

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

24
python/aare/Cluster.py Normal file
View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

Some files were not shown because too many files have changed in this diff Show More