14 Commits

Author SHA1 Message Date
c35f4a7746 some refactoring
Some checks failed
Build on RHEL9 / build (push) Failing after 57s
Build on RHEL8 / build (push) Failing after 1m11s
2025-06-12 11:20:32 +02:00
ba8778cf44 added some more tests and debugged 2025-06-11 11:18:34 +02:00
6bc8b0c4a7 some more tests and refactoring 2025-06-04 09:11:36 +02:00
1369bc780e added some more tests
Some checks failed
Build on RHEL9 / build (push) Failing after 54s
Build on RHEL8 / build (push) Failing after 1m9s
2025-05-30 14:35:37 +02:00
9cfe1ac5e6 function to write to file 2025-05-30 10:54:13 +02:00
6f4cc219b7 MythenDetectorSpecifications cleanup 2025-05-30 10:29:25 +02:00
0d5c6fed61 some file restructuring 2025-05-30 10:00:41 +02:00
54f76100c2 take pow srqt from cmath 2025-05-30 09:18:53 +02:00
e04bf6be30 added minimum mythen file reader
Some checks failed
Build on RHEL9 / build (push) Failing after 59s
Build on RHEL8 / build (push) Failing after 1m10s
2025-05-29 23:17:11 +02:00
bd0ff3d7da function to redistribute histogram to fixed angle bin sizes 2025-05-29 18:19:36 +02:00
df1335529c implementation of FlatField and MythenDetector class 2025-05-28 10:41:33 +02:00
b94be4cbe8 Hdf5Reader supporting reading into frame, ndarray and generic bytestream 2025-05-26 09:06:47 +02:00
6328369ce9 added parameter conversion 2025-05-19 18:30:29 +02:00
67b94eefb0 reading file for initial calibration 2025-05-19 16:46:19 +02:00
117 changed files with 3472 additions and 3000 deletions

View File

@ -1,9 +1,9 @@
name: Build pkgs and deploy if on main
on:
release:
types:
- published
push:
branches:
- main
jobs:
build:

View File

@ -2,10 +2,7 @@ name: Build the package using cmake then documentation
on:
workflow_dispatch:
pull_request:
release:
types:
- published
push:
permissions:
@ -43,7 +40,7 @@ jobs:
run: |
mkdir build
cd build
cmake .. -DAARE_SYSTEM_LIBRARIES=ON -DAARE_PYTHON_BINDINGS=ON -DAARE_DOCS=ON
cmake .. -DAARE_SYSTEM_LIBRARIES=ON -DAARE_DOCS=ON
make -j 2
make docs
@ -58,7 +55,7 @@ jobs:
url: ${{ steps.deployment.outputs.page_url }}
runs-on: ubuntu-latest
needs: build
if: (github.event_name == 'release' && github.event.action == 'published') || (github.event_name == 'workflow_dispatch' )
if: github.ref == 'refs/heads/main'
steps:
- name: Deploy to GitHub Pages
id: deployment

View File

@ -62,6 +62,7 @@ option(AARE_FETCH_CATCH "Use FetchContent to download catch2" ON)
option(AARE_FETCH_JSON "Use FetchContent to download nlohmann::json" ON)
option(AARE_FETCH_ZMQ "Use FetchContent to download libzmq" ON)
option(AARE_FETCH_LMFIT "Use FetchContent to download lmfit" ON)
option(AARE_FETCH_HDF5 "Use FetchContent to download hdf5-devel" OFF)
#Convenience option to use system libraries only (no FetchContent)
@ -73,15 +74,13 @@ if(AARE_SYSTEM_LIBRARIES)
set(AARE_FETCH_CATCH OFF CACHE BOOL "Disabled FetchContent for catch2" FORCE)
set(AARE_FETCH_JSON OFF CACHE BOOL "Disabled FetchContent for nlohmann::json" FORCE)
set(AARE_FETCH_ZMQ OFF CACHE BOOL "Disabled FetchContent for libzmq" FORCE)
set(AARE_FETCH_HDF5 OFF CACHE BOOL "Disabled FetchContent for hdf5" FORCE)
# Still fetch lmfit when setting AARE_SYSTEM_LIBRARIES since this is not available
# on conda-forge
endif()
if(AARE_VERBOSE)
add_compile_definitions(AARE_VERBOSE)
add_compile_definitions(AARE_LOG_LEVEL=aare::logDEBUG5)
else()
add_compile_definitions(AARE_LOG_LEVEL=aare::logERROR)
endif()
if(AARE_CUSTOM_ASSERT)
@ -93,7 +92,6 @@ if(AARE_BENCHMARKS)
endif()
set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
if(AARE_FETCH_LMFIT)
@ -223,6 +221,23 @@ else()
find_package(nlohmann_json 3.11.3 REQUIRED)
endif()
if(AARE_FETCH_HDF5)
message(FATAL_ERROR "Fetching HDF5 via FetchContent is not supported here. Please install it via your system.
For Ubuntu: sudo apt install libhdf5-dev
For Red Hat: sudo dnf install hdf5-devel
For MacOS: brew install hdf5")
else()
find_package(HDF5 QUIET COMPONENTS CXX)
if (HDF5_FOUND)
message(STATUS "Found HDF5: ${HDF5_INCLUDE_DIRS}")
else()
message(FATAL_ERROR "HDF5 was NOT found! Please install it via your system
For Ubuntu: sudo apt install libhdf5-dev
For Red Hat: sudo dnf install hdf5-devel
For MacOS: brew install hdf5")
endif()
endif()
include(GNUInstallDirs)
# If conda build, always set lib dir to 'lib'
@ -335,10 +350,6 @@ if(AARE_ASAN)
)
endif()
if(AARE_TESTS)
enable_testing()
add_subdirectory(tests)
@ -348,6 +359,7 @@ endif()
###------------------------------------------------------------------------------------------
set(PUBLICHEADERS
include/aare/AngleCalibration.hpp
include/aare/ArrayExpr.hpp
include/aare/CalculateEta.hpp
include/aare/Cluster.hpp
@ -362,10 +374,14 @@ set(PUBLICHEADERS
include/aare/Fit.hpp
include/aare/FileInterface.hpp
include/aare/FilePtr.hpp
include/aare/FlatField.hpp
include/aare/Frame.hpp
include/aare/GainMap.hpp
include/aare/geo_helpers.hpp
include/aare/Hdf5FileReader.hpp
include/aare/JungfrauDataFile.hpp
include/aare/MythenDetectorSpecifications.hpp
include/aare/MythenFileReader.hpp
include/aare/NDArray.hpp
include/aare/NDView.hpp
include/aare/NumpyFile.hpp
@ -381,6 +397,7 @@ set(PUBLICHEADERS
set(SourceFiles
${CMAKE_CURRENT_SOURCE_DIR}/src/AngleCalibration.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/CtbRawFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/defs.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.cpp
@ -398,6 +415,7 @@ set(SourceFiles
${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/utils/task.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/utils/ifstream_helpers.cpp
)
@ -416,6 +434,7 @@ target_link_libraries(
PUBLIC
fmt::fmt
nlohmann_json::nlohmann_json
HDF5::HDF5
${STD_FS_LIB} # from helpers.cmake
PRIVATE
aare_compiler_flags
@ -436,12 +455,14 @@ endif()
if(AARE_TESTS)
set(TestSources
${CMAKE_CURRENT_SOURCE_DIR}/src/algorithm.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/AngleCalibration.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/defs.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/decode.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/geo_helpers.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Hdf5FileReader.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NDArray.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NDView.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinder.test.cpp
@ -452,10 +473,10 @@ if(AARE_TESTS)
${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinderMT.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Pedestal.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/MythenFileReader.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/utils/task.test.cpp
)

View File

@ -1,22 +0,0 @@
# Release notes
### head
Features:
- Cluster finder now works with 5x5, 7x7 and 9x9 clusters
### 2025.05.22
Features:
- Added scurve fitting
Bugfixes:
- Fixed crash when opening raw files with large number of data files

View File

@ -1 +1 @@
2025.5.22
0.0.0

View File

@ -41,8 +41,8 @@ BENCHMARK_F(ClusterFixture, Calculate2x2Eta)(benchmark::State &st) {
}
// almost takes double the time
BENCHMARK_F(ClusterFixture, CalculateGeneralEtaFor2x2Cluster)
(benchmark::State &st) {
BENCHMARK_F(ClusterFixture,
CalculateGeneralEtaFor2x2Cluster)(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
Eta2 eta = calculate_eta2<int, 2, 2>(cluster_2x2);
@ -59,8 +59,8 @@ BENCHMARK_F(ClusterFixture, Calculate3x3Eta)(benchmark::State &st) {
}
// almost takes double the time
BENCHMARK_F(ClusterFixture, CalculateGeneralEtaFor3x3Cluster)
(benchmark::State &st) {
BENCHMARK_F(ClusterFixture,
CalculateGeneralEtaFor3x3Cluster)(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
Eta2 eta = calculate_eta2<int, 3, 3>(cluster_3x3);

View File

@ -1,132 +1,136 @@
#include "aare/NDArray.hpp"
#include <benchmark/benchmark.h>
#include "aare/NDArray.hpp"
using aare::NDArray;
constexpr ssize_t size = 1024;
class TwoArrays : public benchmark::Fixture {
public:
NDArray<int, 2> a{{size, size}, 0};
NDArray<int, 2> b{{size, size}, 0};
void SetUp(::benchmark::State &state) {
for (uint32_t i = 0; i < size; i++) {
for (uint32_t j = 0; j < size; j++) {
a(i, j) = i * j + 1;
b(i, j) = i * j + 1;
}
}
public:
NDArray<int,2> a{{size,size},0};
NDArray<int,2> b{{size,size},0};
void SetUp(::benchmark::State& state) {
for(uint32_t i = 0; i < size; i++){
for(uint32_t j = 0; j < size; j++){
a(i, j)= i*j+1;
b(i, j)= i*j+1;
}
}
}
// void TearDown(::benchmark::State& state) {
// }
// void TearDown(::benchmark::State& state) {
// }
};
BENCHMARK_F(TwoArrays, AddWithOperator)(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
NDArray<int, 2> res = a + b;
benchmark::DoNotOptimize(res);
}
BENCHMARK_F(TwoArrays, AddWithOperator)(benchmark::State& st) {
for (auto _ : st) {
// This code gets timed
NDArray<int,2> res = a+b;
benchmark::DoNotOptimize(res);
}
}
BENCHMARK_F(TwoArrays, AddWithIndex)(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
NDArray<int, 2> res(a.shape());
for (uint32_t i = 0; i < a.size(); i++) {
res(i) = a(i) + b(i);
}
benchmark::DoNotOptimize(res);
BENCHMARK_F(TwoArrays, AddWithIndex)(benchmark::State& st) {
for (auto _ : st) {
// This code gets timed
NDArray<int,2> res(a.shape());
for (uint32_t i = 0; i < a.size(); i++) {
res(i) = a(i) + b(i);
}
benchmark::DoNotOptimize(res);
}
}
BENCHMARK_F(TwoArrays, SubtractWithOperator)(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
NDArray<int, 2> res = a - b;
benchmark::DoNotOptimize(res);
}
BENCHMARK_F(TwoArrays, SubtractWithOperator)(benchmark::State& st) {
for (auto _ : st) {
// This code gets timed
NDArray<int,2> res = a-b;
benchmark::DoNotOptimize(res);
}
}
BENCHMARK_F(TwoArrays, SubtractWithIndex)(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
NDArray<int, 2> res(a.shape());
for (uint32_t i = 0; i < a.size(); i++) {
res(i) = a(i) - b(i);
}
benchmark::DoNotOptimize(res);
BENCHMARK_F(TwoArrays, SubtractWithIndex)(benchmark::State& st) {
for (auto _ : st) {
// This code gets timed
NDArray<int,2> res(a.shape());
for (uint32_t i = 0; i < a.size(); i++) {
res(i) = a(i) - b(i);
}
benchmark::DoNotOptimize(res);
}
}
BENCHMARK_F(TwoArrays, MultiplyWithOperator)(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
NDArray<int, 2> res = a * b;
benchmark::DoNotOptimize(res);
}
BENCHMARK_F(TwoArrays, MultiplyWithOperator)(benchmark::State& st) {
for (auto _ : st) {
// This code gets timed
NDArray<int,2> res = a*b;
benchmark::DoNotOptimize(res);
}
}
BENCHMARK_F(TwoArrays, MultiplyWithIndex)(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
NDArray<int, 2> res(a.shape());
for (uint32_t i = 0; i < a.size(); i++) {
res(i) = a(i) * b(i);
}
benchmark::DoNotOptimize(res);
BENCHMARK_F(TwoArrays, MultiplyWithIndex)(benchmark::State& st) {
for (auto _ : st) {
// This code gets timed
NDArray<int,2> res(a.shape());
for (uint32_t i = 0; i < a.size(); i++) {
res(i) = a(i) * b(i);
}
benchmark::DoNotOptimize(res);
}
}
BENCHMARK_F(TwoArrays, DivideWithOperator)(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
NDArray<int, 2> res = a / b;
benchmark::DoNotOptimize(res);
}
BENCHMARK_F(TwoArrays, DivideWithOperator)(benchmark::State& st) {
for (auto _ : st) {
// This code gets timed
NDArray<int,2> res = a/b;
benchmark::DoNotOptimize(res);
}
}
BENCHMARK_F(TwoArrays, DivideWithIndex)(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
NDArray<int, 2> res(a.shape());
for (uint32_t i = 0; i < a.size(); i++) {
res(i) = a(i) / b(i);
}
benchmark::DoNotOptimize(res);
BENCHMARK_F(TwoArrays, DivideWithIndex)(benchmark::State& st) {
for (auto _ : st) {
// This code gets timed
NDArray<int,2> res(a.shape());
for (uint32_t i = 0; i < a.size(); i++) {
res(i) = a(i) / b(i);
}
benchmark::DoNotOptimize(res);
}
}
BENCHMARK_F(TwoArrays, FourAddWithOperator)(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
NDArray<int, 2> res = a + b + a + b;
benchmark::DoNotOptimize(res);
}
BENCHMARK_F(TwoArrays, FourAddWithOperator)(benchmark::State& st) {
for (auto _ : st) {
// This code gets timed
NDArray<int,2> res = a+b+a+b;
benchmark::DoNotOptimize(res);
}
}
BENCHMARK_F(TwoArrays, FourAddWithIndex)(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
NDArray<int, 2> res(a.shape());
for (uint32_t i = 0; i < a.size(); i++) {
res(i) = a(i) + b(i) + a(i) + b(i);
}
benchmark::DoNotOptimize(res);
BENCHMARK_F(TwoArrays, FourAddWithIndex)(benchmark::State& st) {
for (auto _ : st) {
// This code gets timed
NDArray<int,2> res(a.shape());
for (uint32_t i = 0; i < a.size(); i++) {
res(i) = a(i) + b(i) + a(i) + b(i);
}
benchmark::DoNotOptimize(res);
}
}
BENCHMARK_F(TwoArrays, MultiplyAddDivideWithOperator)(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
NDArray<int, 2> res = a * a + b / a;
benchmark::DoNotOptimize(res);
}
BENCHMARK_F(TwoArrays, MultiplyAddDivideWithOperator)(benchmark::State& st) {
for (auto _ : st) {
// This code gets timed
NDArray<int,2> res = a*a+b/a;
benchmark::DoNotOptimize(res);
}
}
BENCHMARK_F(TwoArrays, MultiplyAddDivideWithIndex)(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
NDArray<int, 2> res(a.shape());
for (uint32_t i = 0; i < a.size(); i++) {
res(i) = a(i) * a(i) + b(i) / a(i);
}
benchmark::DoNotOptimize(res);
BENCHMARK_F(TwoArrays, MultiplyAddDivideWithIndex)(benchmark::State& st) {
for (auto _ : st) {
// This code gets timed
NDArray<int,2> res(a.shape());
for (uint32_t i = 0; i < a.size(); i++) {
res(i) = a(i) * a(i) + b(i) / a(i);
}
benchmark::DoNotOptimize(res);
}
}
BENCHMARK_MAIN();

View File

@ -4,5 +4,4 @@ ClusterFile
.. doxygenclass:: aare::ClusterFile
:members:
:undoc-members:
:private-members:
:private-members:

View File

@ -1,47 +0,0 @@
****************
Philosophy
****************
Fast code with a simple interface
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Aare should be fast and efficient, but also easy to use. We strive to keep a simple interface that feels intuitive.
Internally we use C++ for performance and the ability to integrate the library in other programs, but we see most
users using the Python interface.
Live at head
~~~~~~~~~~~~~~~~~~
As a user of the library you should be able to, and is expected to, use the latest version. Bug fixes will rarely be backported
to older releases. By upgrading frequently you will benefit from the latest features and minimize the effort to maintain your scripts/code
by doing several small upgrades instead of one big upgrade.
API
~~~~~~~~~~~~~~~~~~
We aim to keep the API stable and only break it for good reasons. But specially now in the early stages of development
the API will change. On those occasions it will be clearly stated in the release notes. However, the norm should be a
backward compatible API.
Documentation
~~~~~~~~~~~~~~~~~~
Being a library it is important to have a well documented API. We use Doxygen to generate the C++ documentation
and Sphinx for the Python part. Breathe is used to integrate the two into one Sphinx html site. The documentation is built
automatically on release by the CI and published to GitHub pages. In addition to the generated API documentation,
certain classes might need more descriptions of the usage. This is then placed in the .rst files in the docs/src directory.
.. attention::
The code should be well documented, but using descriptive names is more important. In the same spirit
if a function is called `getNumberOfFrames()` you don't need to write a comment saying that it gets the
number of frames.
Dependencies
~~~~~~~~~~~~~~~~~~
Deployment in the scientific community is often tricky. Either due to old OS versions or the lack of package managers.
We strive to keep the dependencies to a minimum and will vendor some libraries to simplify deployment even though it comes
at a cost of build time.

View File

@ -2,21 +2,18 @@ Requirements
==============================================
- C++17 compiler (gcc 8/clang 7)
- CMake 3.15+
- CMake 3.14+
**Internally used libraries**
.. note ::
To save compile time some of the dependencies can also be picked up from the system/conda environment by specifying:
These can also be picked up from the system/conda environment by specifying:
-DAARE_SYSTEM_LIBRARIES=ON during the cmake configuration.
To simplify deployment we build and statically link a few libraries.
- fmt
- lmfit - https://jugit.fz-juelich.de/mlz/lmfit
- nlohmann_json
- pybind11
- fmt
- nlohmann_json
- ZeroMQ
**Extra dependencies for building documentation**

View File

@ -1,86 +0,0 @@
****************
Workflow
****************
This page describes how we develop aare.
GitHub centric
~~~~~~~~~~~~~~~~~~
We use GitHub for all development. Issues and pull requests provide a platform for collaboration as well
as a record of the development process. Even if we discuss things in person, we record the outcome in an issue.
If a particular implementation is chosen over another, the reason should be recorded in the pull request.
Branches
~~~~~~~~~~~~~~~~~~
We aim for an as lightweight branching strategy as possible. Short-lived feature branches are merged back into main.
The main branch is expected to always be in a releasable state. A release is simply a tag on main which provides a
reference and triggers the CI to build the release artifacts (conda, pypi etc.). For large features consider merging
smaller chunks into main as they are completed, rather than waiting for the entire feature to be finished. Worst case
make sure your feature branch merges with main regularly to avoid large merge conflicts later on.
.. note::
The main branch is expected to always work. Feel free to pull from main instead of sticking to a
release
Releases
~~~~~~~~~~~~~~~~~~
Release early, release often. As soon as "enough" new features have been implemented, a release is created.
A release should not be a big thing, rather a routine part of development that does not require any special person or
unfamiliar steps.
Checklists for deployment
~~~~~~~~~~~~~~~~~~
**Feature:**
#. Create a new issue for the feature (label feature)
#. Create a new branch from main.
#. Implement the feature including test and documentation
#. Add the feature to RELEASE.md under head
#. Create a pull request linked to the issue
#. Code is reviewed by at least one other person
#. Once approved, the branch is merged into main
**BugFix:**
Essentially the same as for a feature, if possible start with
a failing test that demonstrates the bug.
#. Create a new issue for the bug (label bug)
#. Create a new branch from main.
#. **Write a test that fails for the bug**
#. Implement the fix
#. **Run the test to ensure it passes**
#. Add the bugfix to RELEASE.md under head
#. Create a pull request linked to the issue.
#. Code is reviewed by at least one other person
#. Once approved, the branch is merged into main
**Release:**
#. Once "enough" new features have been implemented, a release is created
#. Update RELEASE.md with the tag of the release and verify that it is complete
#. Create the release in GitHub describing the new features and bug fixes
#. CI makes magic
**Update documentation only:**
.. attention::
It's possible to update the documentation without changing the code, but take
care since the docs will reflect the code in main and not the latest release.
#. Create a PR to main with the documentation changes
#. Create a pull request linked to the issue.
#. Code is reviewed by at least one other person
#. Once merged you can manually trigger the CI workflow for documentation

View File

@ -63,6 +63,4 @@ AARE
:caption: Developer
:maxdepth: 3
Philosophy
Workflow
Tests

View File

@ -2,24 +2,9 @@
ClusterFile
============
The :class:`ClusterFile` class is the main interface to read and write clusters in aare. Unfortunately the
old file format does not include metadata like the cluster size and the data type. This means that the
user has to know this information from other sources. Specifying the wrong cluster size or data type
will lead to garbage data being read.
.. py:currentmodule:: aare
.. autoclass:: ClusterFile
:members:
:undoc-members:
:inherited-members:
Below is the API of the ClusterFile_Cluster3x3i but all variants share the same API.
.. autoclass:: aare._aare.ClusterFile_Cluster3x3i
:special-members: __init__
:members:
:undoc-members:
:show-inheritance:

View File

@ -2,10 +2,8 @@ ClusterVector
================
The ClusterVector, holds clusters from the ClusterFinder. Since it is templated
in C++ we use a suffix indicating the type of cluster it holds. The suffix follows
the same pattern as for ClusterFile i.e. ``ClusterVector_Cluster3x3i``
for a vector holding 3x3 integer clusters.
in C++ we use a suffix indicating the data type in python. The suffix is
``_i`` for integer, ``_f`` for float, and ``_d`` for double.
At the moment the functionality from python is limited and it is not supported
to push_back clusters to the vector. The intended use case is to pass it to
@ -28,8 +26,7 @@ C++ functions that support the ClusterVector or to view it as a numpy array.
.. py:currentmodule:: aare
.. autoclass:: aare._aare.ClusterVector_Cluster3x3i
:special-members: __init__
.. autoclass:: ClusterVector_i
:members:
:undoc-members:
:show-inheritance:

View File

@ -5,12 +5,9 @@ dependencies:
- anaconda-client
- conda-build
- doxygen
- sphinx
- sphinx=7.1.2
- breathe
- sphinx_rtd_theme
- furo
- zeromq
- pybind11
- numpy
- matplotlib

View File

@ -0,0 +1,164 @@
#pragma once
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <memory>
#include <sstream>
#include <string>
#include <vector>
#include "FlatField.hpp"
#include "MythenDetectorSpecifications.hpp"
#include "MythenFileReader.hpp"
#include "NDArray.hpp"
namespace aare {
using parameters =
std::tuple<std::vector<double>, std::vector<double>, std::vector<double>>;
class AngleCalibration {
public:
AngleCalibration(
std::shared_ptr<MythenDetectorSpecifications> mythen_detector_,
std::shared_ptr<FlatField> flat_field_,
std::shared_ptr<MythenFileReader> mythen_file_reader_);
/** set the histogram bin width [degrees] */
void set_histogram_bin_width(double bin_width);
double get_histogram_bin_width() const;
ssize_t get_new_num_bins() const;
/** reads the historical Detector Group (DG) parameters from file **/
void read_initial_calibration_from_file(const std::string &filename);
std::vector<double> get_centers() const;
std::vector<double> get_conversions() const;
std::vector<double> get_offsets() const;
NDView<double, 1> get_new_photon_counts() const;
NDView<double, 1> get_new_statistical_errors() const;
/** converts DG parameters to easy EE parameters e.g.geometric
* parameters */
parameters convert_to_EE_parameters() const;
std::tuple<double, double, double>
convert_to_EE_parameters(const size_t module_index) const;
std::tuple<double, double, double>
convert_to_EE_parameters(const double center, const double conversion,
const double offset) const;
/** converts DG parameters to BC parameters e.g. best computing
* parameters */
parameters convert_to_BC_parameters() const;
/**
* calculates new histogram with fixed sized angle bins
* for several acquisitions at different detector angles for given frame
* indices
* @param start_frame_index, end_frame_index gives range of frames
*/
void
calculate_fixed_bin_angle_width_histogram(const size_t start_frame_index,
const size_t end_frame_index);
void write_to_file(const std::string &filename,
const bool store_nonzero_bins = false,
const std::filesystem::path &filepath =
std::filesystem::current_path()) const;
/** calculates diffraction angle from EE module parameters (used in
* Beer's Law)
* @param strip_index local strip index of module
*/
double diffraction_angle_from_EE_parameters(
const double module_center_distance, const double normal_distance,
const double angle, const size_t strip_index,
const double distance_to_strip = 0) const;
/** calculates diffraction angle from EE module parameters (used in
* Beer's Law)
* @param center module center
* @param conversion module conversion
* @param offset module offset
* @param strip_index local strip index of module
* @param distance_to_strip distance to strip given by strip_index and
* module -> note needs to be small enough to be in the respective module
*/
double diffraction_angle_from_DG_parameters(
const double center, const double conversion, const double offset,
const size_t strip_index, const double distance_to_strip = 0) const;
/** calculated the strip width expressed as angle [degrees]
* @param strip_index local strip index of module
*/
double angular_strip_width_from_DG_parameters(
const double center, const double conversion, const double offset,
const size_t local_strip_index) const;
double angular_strip_width_from_EE_parameters(
const double module_center_distance, const double normal_distance,
const double angle, const size_t local_strip_index) const;
protected:
/** converts global strip index to local strip index of that module */
size_t global_to_local_strip_index_conversion(
const size_t global_strip_index) const;
/**
* redistributes photon counts with of histogram using one bin per strip
* to histogram with fixed size angle bins
* @param frame MythenFrame storing data from image
* @param bin_counts accumulate new photon counts
* @param new_statistical_weights accumulate new statistical weights
* @param new_errors accumulate new_errors
*/
void redistribute_photon_counts_to_fixed_angle_bins(
const MythenFrame &frame, NDView<double, 1> bin_counts,
NDView<double, 1> new_statistical_weights, NDView<double, 1> new_errors,
NDView<double, 1> inverse_nromalized_flatfield) const;
private:
// TODO: Design maybe have a struct with three vectors, store all three
// sets of parameters as member variables
// TODO: check if interpretation and units are correct
// historical DG parameters
// TODO change to NDArray
std::vector<double> centers; // orthogonal projection of sample onto
// detector (given in strip number) [mm]
// D/pitch
std::vector<double> conversions; // pitch/(normal distance from sample
// to detector (R)) [mm]
// //used for easy conversion
std::vector<double>
offsets; // position of strip zero relative to sample [degrees] phi
// - 180/pi*D/R TODO: expected an arcsin(D/R)?
std::shared_ptr<MythenDetectorSpecifications> mythen_detector;
std::shared_ptr<FlatField> flat_field;
NDArray<double, 1> new_photon_counts;
NDArray<double, 1> new_photon_count_errors;
double histogram_bin_width = 0.0036; // [degrees]
ssize_t num_bins{};
std::shared_ptr<MythenFileReader>
mythen_file_reader; // TODO replace by FileInterface ptr
};
} // namespace aare

View File

@ -1,9 +1,10 @@
#pragma once
#include "aare/defs.hpp"
#include <array>
#include <cstdint>
#include <cstddef>
#include <array>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include "aare/defs.hpp"
namespace aare {
@ -14,9 +15,7 @@ template <typename E, ssize_t Ndim> class ArrayExpr {
auto operator[](size_t i) const { return static_cast<E const &>(*this)[i]; }
auto operator()(size_t i) const { return static_cast<E const &>(*this)[i]; }
auto size() const { return static_cast<E const &>(*this).size(); }
std::array<ssize_t, Ndim> shape() const {
return static_cast<E const &>(*this).shape();
}
std::array<ssize_t, Ndim> shape() const { return static_cast<E const &>(*this).shape(); }
};
template <typename A, typename B, ssize_t Ndim>
@ -48,7 +47,7 @@ class ArraySub : public ArrayExpr<ArraySub<A, B, Ndim>, Ndim> {
};
template <typename A, typename B, ssize_t Ndim>
class ArrayMul : public ArrayExpr<ArrayMul<A, B, Ndim>, Ndim> {
class ArrayMul : public ArrayExpr<ArrayMul<A, B, Ndim>,Ndim> {
const A &arr1_;
const B &arr2_;
@ -75,13 +74,15 @@ class ArrayDiv : public ArrayExpr<ArrayDiv<A, B, Ndim>, Ndim> {
std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); }
};
template <typename A, typename B, ssize_t Ndim>
auto operator+(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
return ArrayAdd<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
}
template <typename A, typename B, ssize_t Ndim>
auto operator-(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
auto operator-(const ArrayExpr<A,Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
return ArraySub<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
}
@ -95,4 +96,6 @@ auto operator/(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
return ArrayDiv<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
}
} // namespace aare

View File

@ -17,8 +17,7 @@ template <class ItemType> class CircularFifo {
public:
CircularFifo() : CircularFifo(100){};
CircularFifo(uint32_t size)
: fifo_size(size), free_slots(size + 1), filled_slots(size + 1) {
CircularFifo(uint32_t size) : fifo_size(size), free_slots(size + 1), filled_slots(size + 1) {
// TODO! how do we deal with alignment for writing? alignas???
// Do we give the user a chance to provide memory locations?
@ -56,8 +55,7 @@ template <class ItemType> class CircularFifo {
bool try_pop_free(ItemType &v) { return free_slots.read(v); }
ItemType pop_value(std::chrono::nanoseconds wait,
std::atomic<bool> &stopped) {
ItemType pop_value(std::chrono::nanoseconds wait, std::atomic<bool> &stopped) {
ItemType v;
while (!filled_slots.read(v) && !stopped) {
std::this_thread::sleep_for(wait);

View File

@ -5,8 +5,6 @@
#include "aare/GainMap.hpp"
#include "aare/NDArray.hpp"
#include "aare/defs.hpp"
#include "aare/logger.hpp"
#include <filesystem>
#include <fstream>
#include <optional>
@ -371,14 +369,10 @@ ClusterFile<ClusterType, Enable>::read_frame_without_cut() {
"Could not read number of clusters");
}
LOG(logDEBUG1) << "Reading " << n_clusters << " clusters from frame "
<< frame_number;
ClusterVector<ClusterType> clusters(n_clusters);
clusters.set_frame_number(frame_number);
clusters.resize(n_clusters);
LOG(logDEBUG1) << "clusters.item_size(): " << clusters.item_size();
clusters.resize(n_clusters);
if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) !=
static_cast<size_t>(n_clusters)) {

View File

@ -21,7 +21,7 @@ class ClusterFileSink {
void process() {
m_stopped = false;
LOG(logDEBUG) << "ClusterFileSink started";
fmt::print("ClusterFileSink started\n");
while (!m_stop_requested || !m_source->isEmpty()) {
if (ClusterVector<ClusterType> *clusters = m_source->frontPtr();
clusters != nullptr) {
@ -41,16 +41,13 @@ class ClusterFileSink {
std::this_thread::sleep_for(m_default_wait);
}
}
LOG(logDEBUG) << "ClusterFileSink stopped";
fmt::print("ClusterFileSink stopped\n");
m_stopped = true;
}
public:
ClusterFileSink(ClusterFinderMT<ClusterType, uint16_t, double> *source,
const std::filesystem::path &fname) {
LOG(logDEBUG) << "ClusterFileSink: "
<< "source: " << source->sink()
<< ", file: " << fname.string();
m_source = source->sink();
m_thread = std::thread(&ClusterFileSink::process, this);
m_file.open(fname, std::ios::binary);

View File

@ -38,11 +38,7 @@ class ClusterFinder {
: m_image_size(image_size), m_nSigma(nSigma),
c2(sqrt((ClusterSizeY + 1) / 2 * (ClusterSizeX + 1) / 2)),
c3(sqrt(ClusterSizeX * ClusterSizeY)),
m_pedestal(image_size[0], image_size[1]), m_clusters(capacity) {
LOG(logDEBUG) << "ClusterFinder: "
<< "image_size: " << image_size[0] << "x" << image_size[1]
<< ", nSigma: " << nSigma << ", capacity: " << capacity;
}
m_pedestal(image_size[0], image_size[1]), m_clusters(capacity) {};
void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
m_pedestal.push(frame);

View File

@ -8,7 +8,6 @@
#include "aare/ClusterFinder.hpp"
#include "aare/NDArray.hpp"
#include "aare/ProducerConsumerQueue.hpp"
#include "aare/logger.hpp"
namespace aare {
@ -124,12 +123,6 @@ class ClusterFinderMT {
size_t capacity = 2000, size_t n_threads = 3)
: m_n_threads(n_threads) {
LOG(logDEBUG1) << "ClusterFinderMT: "
<< "image_size: " << image_size[0] << "x"
<< image_size[1] << ", nSigma: " << nSigma
<< ", capacity: " << capacity
<< ", n_threads: " << n_threads;
for (size_t i = 0; i < n_threads; i++) {
m_cluster_finders.push_back(
std::make_unique<

View File

@ -133,9 +133,9 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
*/
size_t capacity() const { return m_data.capacity(); }
auto begin() const { return m_data.begin(); }
const auto begin() const { return m_data.begin(); }
auto end() const { return m_data.end(); }
const auto end() const { return m_data.end(); }
/**
* @brief Return the size in bytes of a single cluster

View File

@ -1,27 +1,27 @@
#pragma once
#include "aare/FileInterface.hpp"
#include "aare/Frame.hpp"
#include "aare/RawMasterFile.hpp"
#include "aare/Frame.hpp"
#include <filesystem>
#include <fstream>
namespace aare {
namespace aare{
class CtbRawFile {
class CtbRawFile{
RawMasterFile m_master;
std::ifstream m_file;
size_t m_current_frame{0};
size_t m_current_subfile{0};
size_t m_num_subfiles{0};
public:
public:
CtbRawFile(const std::filesystem::path &fname);
void read_into(std::byte *image_buf, DetectorHeader *header = nullptr);
void seek(size_t frame_index); //!< seek to the given frame index
size_t tell() const; //!< get the frame index of the file pointer
void read_into(std::byte *image_buf, DetectorHeader* header = nullptr);
void seek(size_t frame_index); //!< seek to the given frame index
size_t tell() const; //!< get the frame index of the file pointer
// in the specific class we can expose more functionality
@ -29,13 +29,13 @@ class CtbRawFile {
size_t frames_in_file() const;
RawMasterFile master() const;
private:
private:
void find_subfiles();
size_t sub_file_index(size_t frame_index) const {
return frame_index / m_master.max_frames_per_file();
}
void open_data_file(size_t subfile_index);
};
} // namespace aare
}

View File

@ -6,37 +6,31 @@
namespace aare {
// The format descriptor is a single character that specifies the type of the
// data
// The format descriptor is a single character that specifies the type of the data
// - python documentation: https://docs.python.org/3/c-api/arg.html#numbers
// - py::format_descriptor<T>::format() (in pybind11) does not return the same
// format as
// - py::format_descriptor<T>::format() (in pybind11) does not return the same format as
// written in python.org documentation.
// - numpy also doesn't use the same format. and also numpy associates the
// format
// with variable bitdepth types. (e.g. long is int64 on linux64 and int32 on
// win64) https://numpy.org/doc/stable/reference/arrays.scalars.html
// - numpy also doesn't use the same format. and also numpy associates the format
// with variable bitdepth types. (e.g. long is int64 on linux64 and int32 on win64)
// https://numpy.org/doc/stable/reference/arrays.scalars.html
//
// github issue discussing this:
// https://github.com/pybind/pybind11/issues/1908#issuecomment-658358767
//
// [IN LINUX] the difference is for int64 (long) and uint64 (unsigned long). The
// format descriptor is 'q' and 'Q' respectively and in the documentation it is
// 'l' and 'k'.
// [IN LINUX] the difference is for int64 (long) and uint64 (unsigned long). The format
// descriptor is 'q' and 'Q' respectively and in the documentation it is 'l' and 'k'.
// in practice numpy doesn't seem to care when reading buffer info: the library
// interprets 'q' or 'l' as int64 and 'Q' or 'L' as uint64.
// for this reason we decided to use the same format descriptor as pybind to
// avoid any further discrepancies.
// for this reason we decided to use the same format descriptor as pybind to avoid
// any further discrepancies.
// in the following order:
// int8, uint8, int16, uint16, int32, uint32, int64, uint64, float, double
const char DTYPE_FORMAT_DSC[] = {'b', 'B', 'h', 'H', 'i',
'I', 'q', 'Q', 'f', 'd'};
const char DTYPE_FORMAT_DSC[] = {'b', 'B', 'h', 'H', 'i', 'I', 'q', 'Q', 'f', 'd'};
// on linux64 & apple
const char NUMPY_FORMAT_DSC[] = {'b', 'B', 'h', 'H', 'i',
'I', 'l', 'L', 'f', 'd'};
const char NUMPY_FORMAT_DSC[] = {'b', 'B', 'h', 'H', 'i', 'I', 'l', 'L', 'f', 'd'};
/**
* @brief enum class to define the endianess of the system
*/
@ -58,29 +52,12 @@ enum class endian {
*/
class Dtype {
public:
enum TypeIndex {
INT8,
UINT8,
INT16,
UINT16,
INT32,
UINT32,
INT64,
UINT64,
FLOAT,
DOUBLE,
ERROR,
NONE
};
enum TypeIndex { INT8, UINT8, INT16, UINT16, INT32, UINT32, INT64, UINT64, FLOAT, DOUBLE, ERROR, NONE };
uint8_t bitdepth() const;
size_t bytes() const;
std::string format_descr() const {
return std::string(1, DTYPE_FORMAT_DSC[static_cast<int>(m_type)]);
}
std::string numpy_descr() const {
return std::string(1, NUMPY_FORMAT_DSC[static_cast<int>(m_type)]);
}
std::string format_descr() const { return std::string(1, DTYPE_FORMAT_DSC[static_cast<int>(m_type)]); }
std::string numpy_descr() const { return std::string(1, NUMPY_FORMAT_DSC[static_cast<int>(m_type)]); }
explicit Dtype(const std::type_info &t);
explicit Dtype(std::string_view sv);

View File

@ -5,12 +5,12 @@
namespace aare {
/**
* @brief RAII File class for reading, and in the future potentially writing
* image files in various formats. Minimal generic interface. For specail
* fuctions plase use the RawFile or NumpyFile classes directly. Wraps
* FileInterface to abstract the underlying file format
* @note **frame_number** refers the the frame number sent by the detector while
* **frame_index** is the position of the frame in the file
* @brief RAII File class for reading, and in the future potentially writing
* image files in various formats. Minimal generic interface. For specail fuctions
* plase use the RawFile or NumpyFile classes directly.
* Wraps FileInterface to abstract the underlying file format
* @note **frame_number** refers the the frame number sent by the detector while **frame_index**
* is the position of the frame in the file
*/
class File {
std::unique_ptr<FileInterface> file_impl;
@ -25,46 +25,42 @@ class File {
* @throws std::invalid_argument if the file mode is not supported
*
*/
File(const std::filesystem::path &fname, const std::string &mode = "r",
const FileConfig &cfg = {});
/**Since the object is responsible for managing the file we disable copy
* construction */
File(File const &other) = delete;
File(const std::filesystem::path &fname, const std::string &mode="r", const FileConfig &cfg = {});
/**Since the object is responsible for managing the file we disable copy construction */
File(File const &other) = delete;
/**The same goes for copy assignment */
File &operator=(File const &other) = delete;
File& operator=(File const &other) = delete;
File(File &&other) noexcept;
File &operator=(File &&other) noexcept;
File& operator=(File &&other) noexcept;
~File() = default;
// void close(); //!< close the file
Frame
read_frame(); //!< read one frame from the file at the current position
Frame read_frame(size_t frame_index); //!< read one frame at the position
//!< given by frame number
std::vector<Frame> read_n(size_t n_frames); //!< read n_frames from the file
//!< at the current position
Frame read_frame(); //!< read one frame from the file at the current position
Frame read_frame(size_t frame_index); //!< read one frame at the position given by frame number
std::vector<Frame> read_n(size_t n_frames); //!< read n_frames from the file at the current position
void read_into(std::byte *image_buf);
void read_into(std::byte *image_buf, size_t n_frames);
size_t frame_number(); //!< get the frame number at the current position
size_t frame_number(
size_t frame_index); //!< get the frame number at the given frame index
size_t bytes_per_frame() const;
size_t pixels_per_frame() const;
size_t bytes_per_pixel() const;
size_t frame_number(); //!< get the frame number at the current position
size_t frame_number(size_t frame_index); //!< get the frame number at the given frame index
size_t bytes_per_frame() const;
size_t pixels_per_frame() const;
size_t bytes_per_pixel() const;
size_t bitdepth() const;
void seek(size_t frame_index); //!< seek to the given frame index
size_t tell() const; //!< get the frame index of the file pointer
void seek(size_t frame_index); //!< seek to the given frame index
size_t tell() const; //!< get the frame index of the file pointer
size_t total_frames() const;
size_t rows() const;
size_t cols() const;
DetectorType detector_type() const;
};
} // namespace aare

View File

@ -20,10 +20,8 @@ struct FileConfig {
uint64_t rows{};
uint64_t cols{};
bool operator==(const FileConfig &other) const {
return dtype == other.dtype && rows == other.rows &&
cols == other.cols && geometry == other.geometry &&
detector_type == other.detector_type &&
max_frames_per_file == other.max_frames_per_file;
return dtype == other.dtype && rows == other.rows && cols == other.cols && geometry == other.geometry &&
detector_type == other.detector_type && max_frames_per_file == other.max_frames_per_file;
}
bool operator!=(const FileConfig &other) const { return !(*this == other); }
@ -34,11 +32,8 @@ struct FileConfig {
int max_frames_per_file{};
size_t total_frames{};
std::string to_string() const {
return "{ dtype: " + dtype.to_string() +
", rows: " + std::to_string(rows) +
", cols: " + std::to_string(cols) +
", geometry: " + geometry.to_string() +
", detector_type: " + ToString(detector_type) +
return "{ dtype: " + dtype.to_string() + ", rows: " + std::to_string(rows) + ", cols: " + std::to_string(cols) +
", geometry: " + geometry.to_string() + ", detector_type: " + ToString(detector_type) +
", max_frames_per_file: " + std::to_string(max_frames_per_file) +
", total_frames: " + std::to_string(total_frames) + " }";
}
@ -47,8 +42,7 @@ struct FileConfig {
/**
* @brief FileInterface class to define the interface for file operations
* @note parent class for NumpyFile and RawFile
* @note all functions are pure virtual and must be implemented by the derived
* classes
* @note all functions are pure virtual and must be implemented by the derived classes
*/
class FileInterface {
public:
@ -70,20 +64,17 @@ class FileInterface {
* @param n_frames number of frames to read
* @return vector of frames
*/
virtual std::vector<Frame>
read_n(size_t n_frames) = 0; // Is this the right interface?
virtual std::vector<Frame> read_n(size_t n_frames) = 0; // Is this the right interface?
/**
* @brief read one frame from the file at the current position and store it
* in the provided buffer
* @brief read one frame from the file at the current position and store it in the provided buffer
* @param image_buf buffer to store the frame
* @return void
*/
virtual void read_into(std::byte *image_buf) = 0;
/**
* @brief read n_frames from the file at the current position and store them
* in the provided buffer
* @brief read n_frames from the file at the current position and store them in the provided buffer
* @param image_buf buffer to store the frames
* @param n_frames number of frames to read
* @return void
@ -143,6 +134,7 @@ class FileInterface {
*/
virtual size_t bitdepth() const = 0;
virtual DetectorType detector_type() const = 0;
// function to query the data type of the file

View File

@ -12,7 +12,7 @@ class FilePtr {
public:
FilePtr() = default;
FilePtr(const std::filesystem::path &fname, const std::string &mode);
FilePtr(const std::filesystem::path& fname, const std::string& mode);
FilePtr(const FilePtr &) = delete; // we don't want a copy
FilePtr &operator=(const FilePtr &) = delete; // since we handle a resource
FilePtr(FilePtr &&other);

View File

@ -23,19 +23,16 @@ NDArray<double, 1> scurve2(NDView<double, 1> x, NDView<double, 1> par);
} // namespace func
/**
* @brief Estimate the initial parameters for a Gaussian fit
*/
std::array<double, 3> gaus_init_par(const NDView<double, 1> x,
const NDView<double, 1> y);
std::array<double, 3> gaus_init_par(const NDView<double, 1> x, const NDView<double, 1> y);
std::array<double, 2> pol1_init_par(const NDView<double, 1> x,
const NDView<double, 1> y);
std::array<double, 2> pol1_init_par(const NDView<double, 1> x, const NDView<double, 1> y);
std::array<double, 6> scurve_init_par(const NDView<double, 1> x,
const NDView<double, 1> y);
std::array<double, 6> scurve2_init_par(const NDView<double, 1> x,
const NDView<double, 1> y);
std::array<double, 6> scurve_init_par(const NDView<double, 1> x, const NDView<double, 1> y);
std::array<double, 6> scurve2_init_par(const NDView<double, 1> x, const NDView<double, 1> y);
static constexpr int DEFAULT_NUM_THREADS = 4;
@ -46,6 +43,7 @@ static constexpr int DEFAULT_NUM_THREADS = 4;
*/
NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y);
/**
* @brief Fit a 1D Gaussian to each pixel. Data layout [row, col, values]
* @param x x values
@ -56,6 +54,9 @@ NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y);
NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
int n_threads = DEFAULT_NUM_THREADS);
/**
* @brief Fit a 1D Gaussian with error estimates
* @param x x values
@ -66,7 +67,7 @@ NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
*/
void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> par_out, NDView<double, 1> par_err_out,
double &chi2);
double& chi2);
/**
* @brief Fit a 1D Gaussian to each pixel with error estimates. Data layout
@ -79,8 +80,9 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
* @param n_threads number of threads to use
*/
void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out,
NDView<double, 2> chi2_out, int n_threads = DEFAULT_NUM_THREADS);
NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads = DEFAULT_NUM_THREADS
);
NDArray<double, 1> fit_pol1(NDView<double, 1> x, NDView<double, 1> y);
@ -88,33 +90,26 @@ NDArray<double, 3> fit_pol1(NDView<double, 1> x, NDView<double, 3> y,
int n_threads = DEFAULT_NUM_THREADS);
void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> par_out, NDView<double, 1> par_err_out,
double &chi2);
NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2);
// TODO! not sure we need to offer the different version in C++
void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out,
NDView<double, 2> chi2_out, int n_threads = DEFAULT_NUM_THREADS);
NDView<double, 3> par_out, NDView<double, 3> par_err_out,NDView<double, 2> chi2_out,
int n_threads = DEFAULT_NUM_THREADS);
NDArray<double, 1> fit_scurve(NDView<double, 1> x, NDView<double, 1> y);
NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
int n_threads);
void fit_scurve(NDView<double, 1> x, NDView<double, 1> y,
NDView<double, 1> y_err, NDView<double, 1> par_out,
NDView<double, 1> par_err_out, double &chi2);
void fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
NDView<double, 3> y_err, NDView<double, 3> par_out,
NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads);
NDArray<double, 1> fit_scurve(NDView<double, 1> x, NDView<double, 1> y);
NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y, int n_threads);
void fit_scurve(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2);
void fit_scurve(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads);
NDArray<double, 1> fit_scurve2(NDView<double, 1> x, NDView<double, 1> y);
NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
int n_threads);
void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y,
NDView<double, 1> y_err, NDView<double, 1> par_out,
NDView<double, 1> par_err_out, double &chi2);
void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
NDView<double, 3> y_err, NDView<double, 3> par_out,
NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads);
NDArray<double, 1> fit_scurve2(NDView<double, 1> x, NDView<double, 1> y);
NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, int n_threads);
void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2);
void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads);
} // namespace aare

112
include/aare/FlatField.hpp Normal file
View File

@ -0,0 +1,112 @@
/**
* stores flatfield for angle calibration
*/
#pragma once
#include <cmath>
#include <cstdint>
#include <fstream>
#include <iostream>
#include <memory>
#include <sstream>
#include <string>
#include "MythenDetectorSpecifications.hpp"
#include "NDArray.hpp"
namespace aare {
// TODO maybe template now its uint32
class FlatField {
public:
FlatField(std::shared_ptr<MythenDetectorSpecifications> mythen_detector_)
: mythen_detector(mythen_detector_) {
flat_field = NDArray<uint32_t, 1>(
std::array<ssize_t, 1>{mythen_detector->num_strips()}, 0);
}
void read_flatfield_from_file(const std::string &filename) {
std::string word;
uint32_t strip_number{};
try {
std::ifstream file(filename, std::ios_base::in);
if (!file.good()) {
throw std::logic_error("file does not exist");
}
std::stringstream file_buffer;
file_buffer << file.rdbuf();
while (file_buffer >> word) {
strip_number = std::stoi(word);
file_buffer >> word;
if (!mythen_detector->get_bad_channels()[strip_number])
flat_field[strip_number] = std::stod(word);
}
file.close();
} catch (const std::exception &e) {
std::cerr << "Error: " << e.what() << std::endl;
}
}
NDView<uint32_t, 1> get_flatfield() const { return flat_field.view(); }
double calculate_mean(double tolerance = 0.001) const {
auto [sum, count] = std::accumulate(
flat_field.begin(), flat_field.end(),
std::make_pair<double, ssize_t>(0.0, 0),
[&tolerance](std::pair<double, ssize_t> acc, const auto &element) {
return element == 0 ? acc
: std::make_pair(acc.first + element,
acc.second + 1);
});
return sum / count;
}
NDArray<double, 1>
inverse_normalized_flatfield(double tolerance = 0.001) const {
double mean = calculate_mean(tolerance);
NDArray<double, 1> inverse_normalized_flatfield(flat_field.shape());
for (ssize_t i = 0; i < flat_field.size(); ++i) {
inverse_normalized_flatfield[i] =
(flat_field[i] <= tolerance ? 0.0 : mean / flat_field[i]);
if (inverse_normalized_flatfield[i] < tolerance)
mythen_detector->get_bad_channels()[i] = true;
}
return inverse_normalized_flatfield; // TODO: better to have a copy in
// this context but unneccessary
// for angle calibration code
// maybe provide inplace and copy option
// maybe store as member variable access with view
}
NDArray<double, 1> normalized_flatfield(double tolerance = 0.001) const {
double mean = calculate_mean(tolerance);
NDArray<double, 1> normalized_flatfield(flat_field.shape());
for (ssize_t i = 0; i < flat_field.size(); ++i) {
normalized_flatfield[i] = (flat_field[i] == flat_field[i] / mean);
if (normalized_flatfield[i] < tolerance)
mythen_detector->get_bad_channels()[i] = true;
}
return normalized_flatfield;
}
private:
NDArray<uint32_t, 1> flat_field; // TODO: should be 2d
std::shared_ptr<MythenDetectorSpecifications> mythen_detector;
};
} // namespace aare

View File

@ -19,7 +19,7 @@ class Frame {
uint32_t m_cols;
Dtype m_dtype;
std::byte *m_data;
// TODO! Add frame number?
//TODO! Add frame number?
public:
/**
@ -39,7 +39,7 @@ class Frame {
* @param dtype data type of the pixels
*/
Frame(const std::byte *bytes, uint32_t rows, uint32_t cols, Dtype dtype);
~Frame() { delete[] m_data; };
~Frame(){ delete[] m_data; };
/** @warning Copy is disabled to ensure performance when passing
* frames around. Can discuss enabling it.
@ -52,6 +52,7 @@ class Frame {
Frame &operator=(Frame &&other) noexcept;
Frame(Frame &&other) noexcept;
Frame clone() const; //<- Explicit copy
uint32_t rows() const;
@ -92,7 +93,7 @@ class Frame {
if (row >= m_rows || col >= m_cols) {
throw std::out_of_range("Invalid row or column index");
}
// TODO! add tests then reimplement using pixel_ptr
//TODO! add tests then reimplement using pixel_ptr
T data;
std::memcpy(&data, m_data + (row * m_cols + col) * m_dtype.bytes(),
m_dtype.bytes());
@ -101,9 +102,9 @@ class Frame {
/**
* @brief Return an NDView of the frame. This is the preferred way to access
* data in the frame.
*
*
* @tparam T type of the pixels
* @return NDView<T, 2>
* @return NDView<T, 2>
*/
template <typename T> NDView<T, 2> view() {
std::array<ssize_t, 2> shape = {static_cast<ssize_t>(m_rows),
@ -112,7 +113,7 @@ class Frame {
return NDView<T, 2>(data, shape);
}
/**
/**
* @brief Copy the frame data into a new NDArray. This is a deep copy.
*/
template <typename T> NDArray<T> image() {

View File

@ -0,0 +1,212 @@
/************************************************
* @file HD5FFileReader.hpp
* @short HDF5FileReader based on H5File object
***********************************************/
#pragma once
#include "Frame.hpp"
#include "NDArray.hpp"
#include <H5Cpp.h>
#include <array>
#include <cxxabi.h>
#include <optional>
namespace aare {
// return std::type_info
inline const std::type_info &deduce_cpp_type(const H5::DataType datatype) {
if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_UINT8.getId())) {
return typeid(uint8_t);
} else if (H5Tequal(datatype.getId(),
H5::PredType::NATIVE_UINT16.getId())) {
return typeid(uint16_t);
} else if (H5Tequal(datatype.getId(),
H5::PredType::NATIVE_UINT32.getId())) {
return typeid(uint32_t);
} else if (H5Tequal(datatype.getId(),
H5::PredType::NATIVE_UINT64.getId())) {
return typeid(uint64_t);
} else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT8.getId())) {
return typeid(int8_t);
} else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT16.getId())) {
return typeid(int16_t);
} else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT32.getId())) {
return typeid(int32_t);
} else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT64.getId())) {
return typeid(int64_t);
} else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT.getId())) {
return typeid(int);
} else if (H5Tequal(datatype.getId(), H5::PredType::IEEE_F64LE.getId())) {
return typeid(double);
} else if (H5Tequal(datatype.getId(), H5::PredType::IEEE_F32LE.getId())) {
return typeid(float);
} else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_FLOAT.getId())) {
return typeid(float);
} else if (H5Tequal(datatype.getId(),
H5::PredType::NATIVE_DOUBLE.getId())) {
return typeid(float);
} else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_CHAR.getId()) &&
datatype.getId() == H5::PredType::NATIVE_CHAR.getId()) {
return typeid(char);
} else {
throw std::runtime_error("c++ type cannot be deduced");
}
}
struct Subset {
std::vector<hsize_t> shape;
std::vector<hsize_t> offset; // index where data subset should start
};
class HDF5Dataset {
public:
HDF5Dataset(const std::string &datasetname_, const H5::DataSet dataset_)
: datasetname(datasetname_), dataset(dataset_) {
datatype = dataset.getDataType();
cpp_type = &deduce_cpp_type(datatype);
dataspace = dataset.getSpace();
rank = dataspace.getSimpleExtentNdims(); // number of dimensions
shape.resize(rank);
dataspace.getSimpleExtentDims(shape.data(), nullptr);
}
hsize_t get_shape(ssize_t index) const { return shape[index]; }
std::vector<hsize_t> get_shape() const { return shape; }
H5::DataType get_datatype() const { return datatype; }
const std::type_info *get_cpp_type() const { return cpp_type; }
/**
* Reads subset of dataset into the buffer
* e.g. to read one 2d frame pass Subset({shape[1], shape[2]}, {frame_index,
* 0,0})
*/
void
read_into_buffer(std::byte *buffer,
std::optional<const Subset> subset = std::nullopt) const {
if (subset) {
// TODO treat scalar cases
if (static_cast<ssize_t>(subset->offset.size()) != rank) {
throw std::runtime_error("provide an offset for" +
std::to_string(rank) + "dimensions");
}
for (ssize_t i = 0; i < rank; ++i) {
hsize_t size =
i < rank - static_cast<ssize_t>(subset->shape.size())
? 0
: subset->shape[i - (rank - subset->shape.size())];
if ((size + subset->offset[i]) > shape[i]) {
throw std::runtime_error(
"subset is too large or offset is out of bounds");
}
}
H5::DataSpace memspace(static_cast<int>(subset->shape.size()),
subset->shape.data());
dataspace.selectHyperslab(H5S_SELECT_SET, subset->shape.data(),
subset->offset.data());
dataset.read(buffer, datatype, memspace, dataspace);
} else {
dataset.read(buffer, datatype);
}
}
Frame store_as_frame() const {
uint32_t rows{}, cols{};
if (rank == 1) {
rows = 1;
// TODO overflow
cols = static_cast<uint32_t>(shape[0]);
} else if (rank == 2) {
rows = static_cast<uint32_t>(shape[0]);
cols = static_cast<uint32_t>(shape[1]);
} else {
throw std::runtime_error("Frame only supports 2d images");
}
Frame frame(rows, cols, Dtype(*cpp_type));
read_into_buffer(frame.data());
return frame;
}
template <typename T, ssize_t NDim>
NDArray<T, NDim> store_as_ndarray() const {
if (NDim != rank) {
std::cout
<< "Warning: dataset dimension and array dimension mismatch"
<< std::endl; // TODO: replace with log - still valid if we
// want subset
}
if (typeid(T) != *cpp_type) {
std::cout << "Warning: dataset and array type mismatch"
<< std::endl;
}
std::array<ssize_t, NDim> array_shape{};
std::transform(
shape.begin(), shape.end(), array_shape.begin(),
[](const auto dim) { return static_cast<ssize_t>(dim); });
aare::NDArray<T, NDim> dataset_array(array_shape);
read_into_buffer(reinterpret_cast<std::byte *>(dataset_array.data()));
return dataset_array;
}
// getMemDataSize()
private:
std::string datasetname{};
H5::DataSet dataset;
H5::DataSpace dataspace;
H5::DataType datatype;
const std::type_info *cpp_type;
ssize_t rank{};
std::vector<hsize_t> shape{};
};
class HDF5FileReader {
public:
HDF5FileReader() = default;
void open_file(const std::string &filename_) {
filename = filename_;
try {
file = H5::H5File(filename, H5F_ACC_RDONLY);
} catch (H5::Exception &e) {
std::cerr << "Error: " << e.getDetailMsg() << std::endl;
}
}
void close_file() { file.close(); }
HDF5Dataset get_dataset(const std::string &dataset_name) const {
H5::DataSet dataset;
try {
dataset = file.openDataSet(dataset_name);
} catch (H5::Exception &e) {
std::cerr << "Error: " << e.getDetailMsg() << std::endl;
}
// TODO use optional to handle error
return HDF5Dataset(dataset_name, dataset);
}
private:
std::string filename{};
H5::H5File file;
};
} // namespace aare

View File

@ -51,7 +51,7 @@ Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) {
Photon photon;
photon.x = cluster.x;
photon.y = cluster.y;
photon.energy = static_cast<decltype(photon.energy)>(eta.sum);
photon.energy = eta.sum;
// auto ie = nearest_index(m_energy_bins, photon.energy)-1;
// auto ix = nearest_index(m_etabinsx, eta.x)-1;
@ -99,7 +99,7 @@ Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) {
Photon photon;
photon.x = cluster.x;
photon.y = cluster.y;
photon.energy = static_cast<decltype(photon.energy)>(eta.sum);
photon.energy = eta.sum;
// Now do some actual interpolation.
// Find which energy bin the cluster is in

View File

@ -3,113 +3,104 @@
#include <filesystem>
#include <vector>
#include "aare/FileInterface.hpp"
#include "aare/FilePtr.hpp"
#include "aare/NDArray.hpp"
#include "aare/defs.hpp"
#include "aare/NDArray.hpp"
#include "aare/FileInterface.hpp"
namespace aare {
struct JungfrauDataHeader {
struct JungfrauDataHeader{
uint64_t framenum;
uint64_t bunchid;
};
class JungfrauDataFile : public FileInterface {
size_t m_rows{}; //!< number of rows in the image, from find_frame_size();
size_t
m_cols{}; //!< number of columns in the image, from find_frame_size();
size_t m_rows{}; //!< number of rows in the image, from find_frame_size();
size_t m_cols{}; //!< number of columns in the image, from find_frame_size();
size_t m_bytes_per_frame{}; //!< number of bytes per frame excluding header
size_t m_total_frames{}; //!< total number of frames in the series of files
size_t m_offset{}; //!< file index of the first file, allow starting at non
//!< zero file
size_t m_current_file_index{}; //!< The index of the open file
size_t m_current_frame_index{}; //!< The index of the current frame (with
//!< reference to all files)
size_t m_total_frames{}; //!< total number of frames in the series of files
size_t m_offset{}; //!< file index of the first file, allow starting at non zero file
size_t m_current_file_index{}; //!< The index of the open file
size_t m_current_frame_index{}; //!< The index of the current frame (with reference to all files)
std::vector<size_t>
m_last_frame_in_file{}; //!< Used for seeking to the correct file
std::vector<size_t> m_last_frame_in_file{}; //!< Used for seeking to the correct file
std::filesystem::path m_path; //!< path to the files
std::string m_base_name; //!< base name used for formatting file names
FilePtr m_fp; //!< RAII wrapper for a FILE*
FilePtr m_fp; //!< RAII wrapper for a FILE*
using pixel_type = uint16_t;
static constexpr size_t header_size = sizeof(JungfrauDataHeader);
static constexpr size_t n_digits_in_file_index =
6; //!< to format file names
static constexpr size_t header_size = sizeof(JungfrauDataHeader);
static constexpr size_t n_digits_in_file_index = 6; //!< to format file names
public:
JungfrauDataFile(const std::filesystem::path &fname);
std::string base_name()
const; //!< get the base name of the file (without path and extension)
size_t bytes_per_frame() override;
size_t pixels_per_frame() override;
size_t bytes_per_pixel() const;
std::string base_name() const; //!< get the base name of the file (without path and extension)
size_t bytes_per_frame() override;
size_t pixels_per_frame() override;
size_t bytes_per_pixel() const;
size_t bitdepth() const override;
void seek(size_t frame_index)
override; //!< seek to the given frame index (note not byte offset)
size_t tell() override; //!< get the frame index of the file pointer
void seek(size_t frame_index) override; //!< seek to the given frame index (note not byte offset)
size_t tell() override; //!< get the frame index of the file pointer
size_t total_frames() const override;
size_t rows() const override;
size_t cols() const override;
std::array<ssize_t, 2> shape() const;
size_t n_files() const; //!< get the number of files in the series.
std::array<ssize_t,2> shape() const;
size_t n_files() const; //!< get the number of files in the series.
// Extra functions needed for FileInterface
Frame read_frame() override;
Frame read_frame(size_t frame_number) override;
std::vector<Frame> read_n(size_t n_frames = 0) override;
std::vector<Frame> read_n(size_t n_frames=0) override;
void read_into(std::byte *image_buf) override;
void read_into(std::byte *image_buf, size_t n_frames) override;
size_t frame_number(size_t frame_index) override;
DetectorType detector_type() const override;
/**
* @brief Read a single frame from the file into the given buffer.
* @param image_buf buffer to read the frame into. (Note the caller is
* responsible for allocating the buffer)
* @brief Read a single frame from the file into the given buffer.
* @param image_buf buffer to read the frame into. (Note the caller is responsible for allocating the buffer)
* @param header pointer to a JungfrauDataHeader or nullptr to skip header)
*/
void read_into(std::byte *image_buf, JungfrauDataHeader *header = nullptr);
/**
* @brief Read a multiple frames from the file into the given buffer.
* @param image_buf buffer to read the frame into. (Note the caller is
* responsible for allocating the buffer)
* @brief Read a multiple frames from the file into the given buffer.
* @param image_buf buffer to read the frame into. (Note the caller is responsible for allocating the buffer)
* @param n_frames number of frames to read
* @param header pointer to a JungfrauDataHeader or nullptr to skip header)
*/
void read_into(std::byte *image_buf, size_t n_frames,
JungfrauDataHeader *header = nullptr);
/**
void read_into(std::byte *image_buf, size_t n_frames, JungfrauDataHeader *header = nullptr);
/**
* @brief Read a single frame from the file into the given NDArray
* @param image NDArray to read the frame into.
*/
void read_into(NDArray<uint16_t> *image,
JungfrauDataHeader *header = nullptr);
void read_into(NDArray<uint16_t>* image, JungfrauDataHeader* header = nullptr);
JungfrauDataHeader read_header();
std::filesystem::path current_file() const {
return fpath(m_current_file_index + m_offset);
}
std::filesystem::path current_file() const { return fpath(m_current_file_index+m_offset); }
private:
private:
/**
* @brief Find the size of the frame in the file. (256x256, 256x1024,
* 512x1024)
* @brief Find the size of the frame in the file. (256x256, 256x1024, 512x1024)
* @param fname path to the file
* @throws std::runtime_error if the file is empty or the size cannot be
* determined
* @throws std::runtime_error if the file is empty or the size cannot be determined
*/
void find_frame_size(const std::filesystem::path &fname);
void find_frame_size(const std::filesystem::path &fname);
void parse_fname(const std::filesystem::path &fname);
void scan_files();
void open_file(size_t file_index);
std::filesystem::path fpath(size_t frame_index) const;
};
void parse_fname(const std::filesystem::path &fname);
void scan_files();
void open_file(size_t file_index);
std::filesystem::path fpath(size_t frame_index) const;
};
} // namespace aare

View File

@ -0,0 +1,150 @@
#pragma once
#include <cmath>
#include <cstdint>
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>
#include "NDArray.hpp"
namespace aare {
class MythenDetectorSpecifications {
public:
// TODO: constructor that reads from a config file
MythenDetectorSpecifications() {
num_strips_ = max_modules_ * strips_per_module_;
num_connected_modules_ = max_modules_;
bad_channels =
NDArray<bool, 1>(std::array<ssize_t, 1>{num_strips_}, false);
connected_modules = NDArray<bool, 1>(
std::array<ssize_t, 1>{static_cast<ssize_t>(max_modules_)}, true);
}
MythenDetectorSpecifications(const size_t max_modules,
const double exposure_time,
const double bloffset)
: max_modules_(max_modules), exposure_time_(exposure_time),
bloffset_(bloffset) {
num_strips_ = max_modules_ * strips_per_module_;
num_connected_modules_ = max_modules_;
bad_channels =
NDArray<bool, 1>(std::array<ssize_t, 1>{num_strips_}, false);
connected_modules = NDArray<bool, 1>(
std::array<ssize_t, 1>{static_cast<ssize_t>(max_modules_)}, true);
}
void read_bad_channels_from_file(const std::string &filename) {
std::string line;
try {
std::ifstream file(filename, std::ios_base::in);
if (!file.good()) {
throw std::logic_error("file does not exist");
}
while (std::getline(file, line)) {
std::size_t pos = line.find("-");
if (pos == std::string::npos) {
bad_channels(std::stoi(line)) = true;
} else {
size_t line_size = line.size();
for (int i = std::stoi(line.substr(0, pos));
i <= std::stoi(line.substr(pos + 1, line_size - pos));
++i)
bad_channels(i) = true;
}
}
file.close();
} catch (const std::exception &e) {
std::cerr << "Error: " << e.what() << std::endl;
}
}
void read_unconnected_modules_from_file(const std::string &filename) {
std::string line;
try {
std::ifstream file(filename, std::ios_base::in);
if (!file.good()) {
throw std::logic_error("file does not exist");
}
std::stringstream file_buffer;
file_buffer << file.rdbuf();
file_buffer >> line;
num_connected_modules_ -= std::stoi(line);
while (file_buffer >> line) {
size_t module = std::stoi(line);
connected_modules[module] = false;
for (size_t i = module * strips_per_module_;
i < (module + 1) * strips_per_module_; ++i)
bad_channels[i] = true;
}
} catch (const std::exception &e) {
std::cerr << "Error: " << e.what() << std::endl;
}
}
NDView<bool, 1> get_bad_channels() const { return bad_channels.view(); }
NDView<bool, 1> get_connected_modules() const {
return connected_modules.view();
}
static constexpr double pitch() { return pitch_; }
static constexpr size_t strips_per_module() { return strips_per_module_; }
size_t max_modules() const { return max_modules_; }
double exposure_time() const { return exposure_time_; }
double bloffset() const { return bloffset_; }
double dtt0() const { return dtt0_; }
static constexpr double min_angle() { return min_angle_; }
static constexpr double max_angle() { return max_angle_; }
ssize_t num_strips() const { return num_strips_; }
private:
static constexpr size_t strips_per_module_ = 1280;
static constexpr double pitch_ = 0.05; // strip width [mm]
static constexpr double min_angle_ =
-180.0; // maybe shoudnt be static but configurable
static constexpr double max_angle_ = 180.0;
static constexpr double dtt0_ =
0.0; // No idea what this is - probably configurable
size_t max_modules_ = 48;
double exposure_time_ = 5.0; // TODO: could read from acquired file but
// maybe should be configurable
double bloffset_ = 1.532; // what is this? detector offset relative to what?
size_t num_connected_modules_{};
ssize_t num_strips_{};
NDArray<bool, 1> bad_channels;
NDArray<bool, 1> connected_modules; // connected modules
};
} // namespace aare

View File

@ -0,0 +1,82 @@
/************************************************
* @file MythenFileReader.hpp
* @short minimal file reader to read mythen files
***********************************************/
#include <bitset>
#include <filesystem>
#include <string>
#include "Hdf5FileReader.hpp"
#include "NDArray.hpp"
namespace aare {
struct MythenFrame {
NDArray<uint32_t, 1> photon_counts;
double detector_angle{};
// double reference_intensity{}; not needed
std::array<uint8_t, 3> channel_mask{};
};
/** minimal version for a mythen file reader */
class MythenFileReader : public HDF5FileReader {
public:
MythenFileReader(const std::filesystem::path &file_path_,
const std::string &file_prefix_)
: m_base_path(file_path_), file_prefix(file_prefix_) {};
MythenFrame read_frame(ssize_t frame_index) {
// TODO not a good design fixed number of digits in file name for frame
// number -> pad with zeros
// not even sure if files have the same name
std::string current_file_name =
m_base_path / (file_prefix + std::to_string(frame_index) + ".h5");
MythenFrame frame;
open_file(current_file_name);
auto dataset_photon_count =
get_dataset("/entry/instrument/detector/data");
frame.photon_counts =
dataset_photon_count.store_as_ndarray<uint32_t, 1>();
++frame.photon_counts; // Why though?
auto dataset_detector_angle =
get_dataset("/entry/instrument/NDAttributes/DetectorAngle");
dataset_detector_angle.read_into_buffer(
reinterpret_cast<std::byte *>(&frame.detector_angle));
auto dataset_channel_number =
get_dataset("/entry/instrument/NDAttributes/CounterMask");
uint8_t channel_number;
dataset_channel_number.read_into_buffer(
reinterpret_cast<std::byte *>(&channel_number));
std::bitset<3> binary_channel_numbers(channel_number); // 1 0 0
// binary_channel_numbers.flip(); // TODO not sure where most
// significant
// bit is ask Anna again
frame.channel_mask = std::array<uint8_t, 3>{binary_channel_numbers[0],
binary_channel_numbers[1],
binary_channel_numbers[2]};
close_file();
return frame;
}
private:
std::filesystem::path m_base_path{};
std::string file_prefix{};
};
} // namespace aare

View File

@ -33,7 +33,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
* @brief Default constructor. Will construct an empty NDArray.
*
*/
NDArray() : shape_(), strides_(c_strides<Ndim>(shape_)), data_(nullptr){};
NDArray() : shape_(), strides_(c_strides<Ndim>(shape_)), data_(nullptr) {};
/**
* @brief Construct a new NDArray object with a given shape.
@ -139,6 +139,9 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
NDArray<bool, Ndim> operator>(const NDArray &other);
bool equals(const NDArray<T, Ndim> &other,
const T tolerance = std::numeric_limits<T>::epsilon()) const;
bool operator==(const NDArray &other) const;
bool operator!=(const NDArray &other) const;
@ -380,6 +383,7 @@ NDArray<T, Ndim> NDArray<T, Ndim>::operator*(const T &value) {
result *= value;
return result;
}
// template <typename T, ssize_t Ndim> void NDArray<T, Ndim>::Print() {
// if (shape_[0] < 20 && shape_[1] < 20)
// Print_all();
@ -422,7 +426,7 @@ template <typename T, ssize_t Ndim>
void save(NDArray<T, Ndim> &img, std::string &pathname) {
std::ofstream f;
f.open(pathname, std::ios::binary);
f.write(img.buffer(), img.size() * sizeof(T));
f.write(reinterpret_cast<char *>(img.buffer()), img.size() * sizeof(T));
f.close();
}
@ -432,9 +436,38 @@ NDArray<T, Ndim> load(const std::string &pathname,
NDArray<T, Ndim> img{shape};
std::ifstream f;
f.open(pathname, std::ios::binary);
f.read(img.buffer(), img.size() * sizeof(T));
f.read(reinterpret_cast<char *>(img.buffer()), img.size() * sizeof(T));
f.close();
return img;
}
template <typename T, ssize_t Ndim = 1>
NDArray<T, Ndim> load_non_binary_file(const std::string &filename,
const std::array<ssize_t, Ndim> shape) {
std::string word;
NDArray<T, Ndim> array(shape);
try {
std::ifstream file(filename, std::ios_base::in);
if (!file.good()) {
throw std::logic_error("file does not exist");
}
std::stringstream file_buffer;
file_buffer << file.rdbuf();
ssize_t counter = 0;
while (file_buffer >> word && counter < size) {
array[counter] = static_cast<T>(
std::stod(word)); // TODO change for different Types
++counter;
}
file.close();
} catch (const std::exception &e) {
std::cerr << "Error: " << e.what() << std::endl;
}
return array;
}
} // namespace aare

View File

@ -6,6 +6,7 @@
#include <array>
#include <cassert>
#include <cstdint>
#include <fstream>
#include <functional>
#include <iomanip>
#include <iostream>
@ -105,6 +106,20 @@ class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
return true;
}
bool equals(const NDView<T, Ndim> &other, const T tolerance) const {
if (shape_ != other.shape_)
return false;
using SignedT = typename make_signed<T>::type;
for (uint32_t i = 0; i != size_; ++i)
if (std::abs(static_cast<SignedT>(buffer_[i]) -
static_cast<SignedT>(other.buffer_[i])) > tolerance)
return false;
return true;
}
NDView &operator+=(const T val) { return elemenwise(val, std::plus<T>()); }
NDView &operator-=(const T val) { return elemenwise(val, std::minus<T>()); }
NDView &operator*=(const T val) {
@ -206,4 +221,12 @@ template <typename T> NDView<T, 1> make_view(std::vector<T> &vec) {
return NDView<T, 1>(vec.data(), {static_cast<ssize_t>(vec.size())});
}
template <typename T, ssize_t Ndim>
void save(NDView<T, Ndim> img, const std::string &pathname) {
std::ofstream f;
f.open(pathname, std::ios::binary);
f.write(reinterpret_cast<char *>(img.data()), img.size() * sizeof(T));
f.close();
}
} // namespace aare

View File

@ -1,8 +1,9 @@
#pragma once
#include "aare/Dtype.hpp"
#include "aare/defs.hpp"
#include "aare/FileInterface.hpp"
#include "aare/NumpyHelpers.hpp"
#include "aare/defs.hpp"
#include <filesystem>
#include <iostream>
@ -10,12 +11,13 @@
namespace aare {
/**
* @brief NumpyFile class to read and write numpy files
* @note derived from FileInterface
* @note implements all the pure virtual functions from FileInterface
* @note documentation for the functions can also be found in the FileInterface
* class
* @note documentation for the functions can also be found in the FileInterface class
*/
class NumpyFile : public FileInterface {
@ -26,35 +28,26 @@ class NumpyFile : public FileInterface {
* @param mode file mode (r, w)
* @param cfg file configuration
*/
explicit NumpyFile(const std::filesystem::path &fname,
const std::string &mode = "r", FileConfig cfg = {});
explicit NumpyFile(const std::filesystem::path &fname, const std::string &mode = "r", FileConfig cfg = {});
void write(Frame &frame);
Frame read_frame() override { return get_frame(this->current_frame++); }
Frame read_frame(size_t frame_number) override {
return get_frame(frame_number);
}
Frame read_frame(size_t frame_number) override { return get_frame(frame_number); }
std::vector<Frame> read_n(size_t n_frames) override;
void read_into(std::byte *image_buf) override {
return get_frame_into(this->current_frame++, image_buf);
}
void read_into(std::byte *image_buf) override { return get_frame_into(this->current_frame++, image_buf); }
void read_into(std::byte *image_buf, size_t n_frames) override;
size_t frame_number(size_t frame_index) override { return frame_index; };
size_t bytes_per_frame() override;
size_t pixels_per_frame() override;
void seek(size_t frame_number) override {
this->current_frame = frame_number;
}
void seek(size_t frame_number) override { this->current_frame = frame_number; }
size_t tell() override { return this->current_frame; }
size_t total_frames() const override { return m_header.shape[0]; }
size_t rows() const override { return m_header.shape[1]; }
size_t cols() const override { return m_header.shape[2]; }
size_t bitdepth() const override { return m_header.dtype.bitdepth(); }
DetectorType detector_type() const override {
return DetectorType::Unknown;
}
DetectorType detector_type() const override { return DetectorType::Unknown; }
/**
* @brief get the data type of the numpy file
@ -77,8 +70,7 @@ class NumpyFile : public FileInterface {
template <typename T, size_t NDim> NDArray<T, NDim> load() {
NDArray<T, NDim> arr(make_shape<NDim>(m_header.shape));
if (fseek(fp, static_cast<long>(header_size), SEEK_SET)) {
throw std::runtime_error(LOCATION +
"Error seeking to the start of the data");
throw std::runtime_error(LOCATION + "Error seeking to the start of the data");
}
size_t rc = fread(arr.data(), sizeof(T), arr.size(), fp);
if (rc != static_cast<size_t>(arr.size())) {
@ -86,20 +78,16 @@ class NumpyFile : public FileInterface {
}
return arr;
}
template <typename A, typename TYPENAME, A Ndim>
void write(NDView<TYPENAME, Ndim> &frame) {
template <typename A, typename TYPENAME, A Ndim> void write(NDView<TYPENAME, Ndim> &frame) {
write_impl(frame.data(), frame.total_bytes());
}
template <typename A, typename TYPENAME, A Ndim>
void write(NDArray<TYPENAME, Ndim> &frame) {
template <typename A, typename TYPENAME, A Ndim> void write(NDArray<TYPENAME, Ndim> &frame) {
write_impl(frame.data(), frame.total_bytes());
}
template <typename A, typename TYPENAME, A Ndim>
void write(NDView<TYPENAME, Ndim> &&frame) {
template <typename A, typename TYPENAME, A Ndim> void write(NDView<TYPENAME, Ndim> &&frame) {
write_impl(frame.data(), frame.total_bytes());
}
template <typename A, typename TYPENAME, A Ndim>
void write(NDArray<TYPENAME, Ndim> &&frame) {
template <typename A, typename TYPENAME, A Ndim> void write(NDArray<TYPENAME, Ndim> &&frame) {
write_impl(frame.data(), frame.total_bytes());
}

View File

@ -40,18 +40,15 @@ bool parse_bool(const std::string &in);
std::string get_value_from_map(const std::string &mapstr);
std::unordered_map<std::string, std::string>
parse_dict(std::string in, const std::vector<std::string> &keys);
std::unordered_map<std::string, std::string> parse_dict(std::string in, const std::vector<std::string> &keys);
template <typename T, size_t N>
bool in_array(T val, const std::array<T, N> &arr) {
template <typename T, size_t N> bool in_array(T val, const std::array<T, N> &arr) {
return std::find(std::begin(arr), std::end(arr), val) != std::end(arr);
}
bool is_digits(const std::string &str);
aare::Dtype parse_descr(std::string typestring);
size_t write_header(const std::filesystem::path &fname,
const NumpyHeader &header);
size_t write_header(const std::filesystem::path &fname, const NumpyHeader &header);
size_t write_header(std::ostream &out, const NumpyHeader &header);
} // namespace NumpyHelpers

View File

@ -18,15 +18,15 @@ template <typename SUM_TYPE = double> class Pedestal {
uint32_t m_samples;
NDArray<uint32_t, 2> m_cur_samples;
// TODO! in case of int needs to be changed to uint64_t
//TODO! in case of int needs to be changed to uint64_t
NDArray<SUM_TYPE, 2> m_sum;
NDArray<SUM_TYPE, 2> m_sum2;
// Cache mean since it is used over and over in the ClusterFinder
// This optimization is related to the access pattern of the ClusterFinder
// Relies on having more reads than pushes to the pedestal
NDArray<SUM_TYPE, 2> m_mean;
//Cache mean since it is used over and over in the ClusterFinder
//This optimization is related to the access pattern of the ClusterFinder
//Relies on having more reads than pushes to the pedestal
NDArray<SUM_TYPE, 2> m_mean;
public:
Pedestal(uint32_t rows, uint32_t cols, uint32_t n_samples = 1000)
@ -42,7 +42,9 @@ template <typename SUM_TYPE = double> class Pedestal {
}
~Pedestal() = default;
NDArray<SUM_TYPE, 2> mean() { return m_mean; }
NDArray<SUM_TYPE, 2> mean() {
return m_mean;
}
SUM_TYPE mean(const uint32_t row, const uint32_t col) const {
return m_mean(row, col);
@ -69,6 +71,8 @@ template <typename SUM_TYPE = double> class Pedestal {
return variance_array;
}
NDArray<SUM_TYPE, 2> std() {
NDArray<SUM_TYPE, 2> standard_deviation_array({m_rows, m_cols});
for (uint32_t i = 0; i < m_rows * m_cols; i++) {
@ -79,6 +83,8 @@ template <typename SUM_TYPE = double> class Pedestal {
return standard_deviation_array;
}
void clear() {
m_sum = 0;
m_sum2 = 0;
@ -86,12 +92,16 @@ template <typename SUM_TYPE = double> class Pedestal {
m_mean = 0;
}
void clear(const uint32_t row, const uint32_t col) {
m_sum(row, col) = 0;
m_sum2(row, col) = 0;
m_cur_samples(row, col) = 0;
m_mean(row, col) = 0;
}
template <typename T> void push(NDView<T, 2> frame) {
assert(frame.size() == m_rows * m_cols);
@ -112,7 +122,7 @@ template <typename SUM_TYPE = double> class Pedestal {
/**
* Push but don't update the cached mean. Speeds up the process
* when initializing the pedestal.
*
*
*/
template <typename T> void push_no_update(NDView<T, 2> frame) {
assert(frame.size() == m_rows * m_cols);
@ -130,6 +140,9 @@ template <typename SUM_TYPE = double> class Pedestal {
}
}
template <typename T> void push(Frame &frame) {
assert(frame.rows() == static_cast<size_t>(m_rows) &&
frame.cols() == static_cast<size_t>(m_cols));
@ -157,8 +170,7 @@ template <typename SUM_TYPE = double> class Pedestal {
m_sum(row, col) += val - m_sum(row, col) / m_samples;
m_sum2(row, col) += val * val - m_sum2(row, col) / m_samples;
}
// Since we just did a push we know that m_cur_samples(row, col) is at
// least 1
//Since we just did a push we know that m_cur_samples(row, col) is at least 1
m_mean(row, col) = m_sum(row, col) / m_cur_samples(row, col);
}
@ -171,8 +183,7 @@ template <typename SUM_TYPE = double> class Pedestal {
m_cur_samples(row, col)++;
} else {
m_sum(row, col) += val - m_sum(row, col) / m_cur_samples(row, col);
m_sum2(row, col) +=
val * val - m_sum2(row, col) / m_cur_samples(row, col);
m_sum2(row, col) += val * val - m_sum2(row, col) / m_cur_samples(row, col);
}
}
@ -180,16 +191,19 @@ template <typename SUM_TYPE = double> class Pedestal {
* @brief Update the mean of the pedestal. This is used after having done
* push_no_update. It is not necessary to call this function after push.
*/
void update_mean() { m_mean = m_sum / m_cur_samples; }
void update_mean(){
m_mean = m_sum / m_cur_samples;
}
template <typename T>
void push_fast(const uint32_t row, const uint32_t col, const T val_) {
// Assume we reached the steady state where all pixels have
// m_samples samples
template<typename T>
void push_fast(const uint32_t row, const uint32_t col, const T val_){
//Assume we reached the steady state where all pixels have
//m_samples samples
SUM_TYPE val = static_cast<SUM_TYPE>(val_);
m_sum(row, col) += val - m_sum(row, col) / m_samples;
m_sum2(row, col) += val * val - m_sum2(row, col) / m_samples;
m_mean(row, col) = m_sum(row, col) / m_samples;
}
};
} // namespace aare

View File

@ -1,7 +1,7 @@
#pragma once
#include "aare/NDArray.hpp"
#include "aare/defs.hpp"
#include "aare/NDArray.hpp"
namespace aare {
@ -10,11 +10,11 @@ NDArray<ssize_t, 2> GenerateMoench05PixelMap();
NDArray<ssize_t, 2> GenerateMoench05PixelMap1g();
NDArray<ssize_t, 2> GenerateMoench05PixelMapOld();
// Matterhorn02
NDArray<ssize_t, 2> GenerateMH02SingleCounterPixelMap();
//Matterhorn02
NDArray<ssize_t, 2>GenerateMH02SingleCounterPixelMap();
NDArray<ssize_t, 3> GenerateMH02FourCounterPixelMap();
// Eiger
NDArray<ssize_t, 2> GenerateEigerFlipRowsPixelMap();
//Eiger
NDArray<ssize_t, 2>GenerateEigerFlipRowsPixelMap();
} // namespace aare

View File

@ -18,9 +18,9 @@
// @author Jordan DeLong (delong.j@fb.com)
// Changes made by PSD Detector Group:
// Copied: Line 34 constexpr std::size_t hardware_destructive_interference_size
// = 128; from folly/lang/Align.h Changed extension to .hpp Changed namespace to
// aare
// Copied: Line 34 constexpr std::size_t hardware_destructive_interference_size = 128; from folly/lang/Align.h
// Changed extension to .hpp
// Changed namespace to aare
#pragma once
@ -45,14 +45,15 @@ template <class T> struct ProducerConsumerQueue {
ProducerConsumerQueue(const ProducerConsumerQueue &) = delete;
ProducerConsumerQueue &operator=(const ProducerConsumerQueue &) = delete;
ProducerConsumerQueue(ProducerConsumerQueue &&other) {
ProducerConsumerQueue(ProducerConsumerQueue &&other){
size_ = other.size_;
records_ = other.records_;
other.records_ = nullptr;
readIndex_ = other.readIndex_.load(std::memory_order_acquire);
writeIndex_ = other.writeIndex_.load(std::memory_order_acquire);
}
ProducerConsumerQueue &operator=(ProducerConsumerQueue &&other) {
ProducerConsumerQueue &operator=(ProducerConsumerQueue &&other){
size_ = other.size_;
records_ = other.records_;
other.records_ = nullptr;
@ -60,17 +61,16 @@ template <class T> struct ProducerConsumerQueue {
writeIndex_ = other.writeIndex_.load(std::memory_order_acquire);
return *this;
}
ProducerConsumerQueue() : ProducerConsumerQueue(2){};
ProducerConsumerQueue():ProducerConsumerQueue(2){};
// size must be >= 2.
//
// Also, note that the number of usable slots in the queue at any
// given time is actually (size-1), so if you start with an empty queue,
// isFull() will return true after size-1 insertions.
explicit ProducerConsumerQueue(uint32_t size)
: size_(size),
records_(static_cast<T *>(std::malloc(sizeof(T) * size))),
readIndex_(0), writeIndex_(0) {
: size_(size), records_(static_cast<T *>(std::malloc(sizeof(T) * size))), readIndex_(0), writeIndex_(0) {
assert(size >= 2);
if (!records_) {
throw std::bad_alloc();
@ -154,8 +154,7 @@ template <class T> struct ProducerConsumerQueue {
}
bool isEmpty() const {
return readIndex_.load(std::memory_order_acquire) ==
writeIndex_.load(std::memory_order_acquire);
return readIndex_.load(std::memory_order_acquire) == writeIndex_.load(std::memory_order_acquire);
}
bool isFull() const {
@ -176,8 +175,7 @@ template <class T> struct ProducerConsumerQueue {
// be removing items concurrently).
// * It is undefined to call this from any other thread.
size_t sizeGuess() const {
int ret = writeIndex_.load(std::memory_order_acquire) -
readIndex_.load(std::memory_order_acquire);
int ret = writeIndex_.load(std::memory_order_acquire) - readIndex_.load(std::memory_order_acquire);
if (ret < 0) {
ret += size_;
}
@ -194,7 +192,7 @@ template <class T> struct ProducerConsumerQueue {
// const uint32_t size_;
uint32_t size_;
// T *const records_;
T *records_;
T* records_;
alignas(hardware_destructive_interference_size) AtomicIndex readIndex_;
alignas(hardware_destructive_interference_size) AtomicIndex writeIndex_;

View File

@ -1,10 +1,11 @@
#pragma once
#include "aare/FileInterface.hpp"
#include "aare/RawMasterFile.hpp"
#include "aare/Frame.hpp"
#include "aare/NDArray.hpp" //for pixel map
#include "aare/RawMasterFile.hpp"
#include "aare/RawSubFile.hpp"
#include <optional>
namespace aare {
@ -29,11 +30,22 @@ struct ModuleConfig {
* Consider using that unless you need raw file specific functionality.
*/
class RawFile : public FileInterface {
std::vector<std::unique_ptr<RawSubFile>> m_subfiles;
size_t n_subfiles{}; //f0,f1...fn
size_t n_subfile_parts{}; // d0,d1...dn
//TODO! move to vector of SubFile instead of pointers
std::vector<std::vector<RawSubFile *>> subfiles; //subfiles[f0,f1...fn][d0,d1...dn]
// std::vector<xy> positions;
ModuleConfig cfg{0, 0};
RawMasterFile m_master;
size_t m_current_frame{};
size_t m_current_subfile{};
// std::vector<ModuleGeometry> m_module_pixel_0;
// size_t m_rows{};
// size_t m_cols{};
DetectorGeometry m_geometry;
public:
@ -44,7 +56,7 @@ class RawFile : public FileInterface {
*/
RawFile(const std::filesystem::path &fname, const std::string &mode = "r");
virtual ~RawFile() override = default;
virtual ~RawFile() override;
Frame read_frame() override;
Frame read_frame(size_t frame_number) override;
@ -52,10 +64,10 @@ class RawFile : public FileInterface {
void read_into(std::byte *image_buf) override;
void read_into(std::byte *image_buf, size_t n_frames) override;
// TODO! do we need to adapt the API?
//TODO! do we need to adapt the API?
void read_into(std::byte *image_buf, DetectorHeader *header);
void read_into(std::byte *image_buf, size_t n_frames,
DetectorHeader *header);
void read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *header);
size_t frame_number(size_t frame_index) override;
size_t bytes_per_frame() override;
@ -68,21 +80,24 @@ class RawFile : public FileInterface {
size_t cols() const override;
size_t bitdepth() const override;
xy geometry();
size_t n_modules() const;
size_t n_mod() const;
RawMasterFile master() const;
DetectorType detector_type() const override;
private:
/**
* @brief read the frame at the given frame index into the image buffer
* @param frame_number frame number to read
* @param image_buf buffer to store the frame
*/
void get_frame_into(size_t frame_index, std::byte *frame_buffer,
DetectorHeader *header = nullptr);
void get_frame_into(size_t frame_index, std::byte *frame_buffer, DetectorHeader *header = nullptr);
/**
* @brief get the frame at the given frame index
@ -91,6 +106,8 @@ class RawFile : public FileInterface {
*/
Frame get_frame(size_t frame_index);
/**
* @brief read the header of the file
* @param fname path to the data subfile
@ -98,8 +115,12 @@ class RawFile : public FileInterface {
*/
static DetectorHeader read_header(const std::filesystem::path &fname);
// void update_geometry_with_roi();
int find_number_of_subfiles();
void open_subfiles();
void find_geometry();
};
} // namespace aare

View File

@ -45,7 +45,7 @@ class ScanParameters {
int m_start = 0;
int m_stop = 0;
int m_step = 0;
// TODO! add settleTime, requires string to time conversion
//TODO! add settleTime, requires string to time conversion
public:
ScanParameters(const std::string &par);
@ -61,6 +61,7 @@ class ScanParameters {
void increment_stop();
};
/**
* @brief Class for parsing a master file either in our .json format or the old
* .raw format
@ -100,6 +101,7 @@ class RawMasterFile {
std::optional<ROI> m_roi;
public:
RawMasterFile(const std::filesystem::path &fpath);
@ -119,7 +121,6 @@ class RawMasterFile {
size_t total_frames_expected() const;
xy geometry() const;
size_t n_modules() const;
std::optional<size_t> analog_samples() const;
std::optional<size_t> digital_samples() const;
@ -127,8 +128,10 @@ class RawMasterFile {
std::optional<size_t> number_of_rows() const;
std::optional<uint8_t> quad() const;
std::optional<ROI> roi() const;
ScanParameters scan_parameters() const;
private:

View File

@ -10,34 +10,23 @@
namespace aare {
/**
* @brief Class to read a singe subfile written in .raw format. Used from
* RawFile to read the entire detector. Can be used directly to read part of the
* image.
* @brief Class to read a singe subfile written in .raw format. Used from RawFile to read
* the entire detector. Can be used directly to read part of the image.
*/
class RawSubFile {
protected:
std::ifstream m_file;
DetectorType m_detector_type;
size_t m_bitdepth;
std::filesystem::path m_path; //!< path to the subfile
std::string m_base_name; //!< base name used for formatting file names
size_t m_offset{}; //!< file index of the first file, allow starting at non
//!< zero file
size_t m_total_frames{}; //!< total number of frames in the series of files
std::filesystem::path m_fname;
size_t m_rows{};
size_t m_cols{};
size_t m_bytes_per_frame{};
int m_module_index{};
size_t m_current_file_index{}; //!< The index of the open file
size_t m_current_frame_index{}; //!< The index of the current frame (with
//!< reference to all files)
std::vector<size_t>
m_last_frame_in_file{}; //!< Used for seeking to the correct file
size_t m_num_frames{};
uint32_t m_pos_row{};
uint32_t m_pos_col{};
std::optional<NDArray<ssize_t, 2>> m_pixel_map;
public:
@ -51,14 +40,12 @@ class RawSubFile {
* @throws std::invalid_argument if the detector,type pair is not supported
*/
RawSubFile(const std::filesystem::path &fname, DetectorType detector,
size_t rows, size_t cols, size_t bitdepth, uint32_t pos_row = 0,
uint32_t pos_col = 0);
size_t rows, size_t cols, size_t bitdepth, uint32_t pos_row = 0, uint32_t pos_col = 0);
~RawSubFile() = default;
/**
* @brief Seek to the given frame number
* @note Puts the file pointer at the start of the header, not the start of
* the data
* @note Puts the file pointer at the start of the header, not the start of the data
* @param frame_index frame position in file to seek to
* @throws std::runtime_error if the frame number is out of range
*/
@ -66,30 +53,26 @@ class RawSubFile {
size_t tell();
void read_into(std::byte *image_buf, DetectorHeader *header = nullptr);
void read_into(std::byte *image_buf, size_t n_frames,
DetectorHeader *header = nullptr);
void read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *header= nullptr);
void get_part(std::byte *buffer, size_t frame_index);
void read_header(DetectorHeader *header);
size_t rows() const;
size_t cols() const;
size_t frame_number(size_t frame_index);
size_t bytes_per_frame() const { return m_bytes_per_frame; }
size_t pixels_per_frame() const { return m_rows * m_cols; }
size_t bytes_per_pixel() const { return m_bitdepth / bits_per_byte; }
size_t frames_in_file() const { return m_total_frames; }
size_t frames_in_file() const { return m_num_frames; }
private:
template <typename T> void read_with_map(std::byte *image_buf);
private:
template <typename T>
void read_with_map(std::byte *image_buf);
void parse_fname(const std::filesystem::path &fname);
void scan_files();
void open_file(size_t file_index);
std::filesystem::path fpath(size_t file_index) const;
};
} // namespace aare

View File

@ -38,13 +38,11 @@ template <typename T> class VarClusterFinder {
bool use_noise_map = false;
int peripheralThresholdFactor_ = 5;
int current_label;
const std::array<int, 4> di{
{0, -1, -1, -1}}; // row ### 8-neighbour by scaning from left to right
const std::array<int, 4> dj{
{-1, -1, 0, 1}}; // col ### 8-neighbour by scaning from top to bottom
const std::array<int, 4> di{{0, -1, -1, -1}}; // row ### 8-neighbour by scaning from left to right
const std::array<int, 4> dj{{-1, -1, 0, 1}}; // col ### 8-neighbour by scaning from top to bottom
const std::array<int, 8> di_{{0, 0, -1, 1, -1, 1, -1, 1}}; // row
const std::array<int, 8> dj_{{-1, 1, 0, 0, 1, -1, -1, 1}}; // col
std::map<int, int> child; // heirachy: key: child; val: parent
std::map<int, int> child; // heirachy: key: child; val: parent
std::unordered_map<int, Hit> h_size;
std::vector<Hit> hits;
// std::vector<std::vector<int16_t>> row
@ -52,8 +50,7 @@ template <typename T> class VarClusterFinder {
public:
VarClusterFinder(Shape<2> shape, T threshold)
: shape_(shape), labeled_(shape, 0), peripheral_labeled_(shape, 0),
binary_(shape), threshold_(threshold) {
: shape_(shape), labeled_(shape, 0), peripheral_labeled_(shape, 0), binary_(shape), threshold_(threshold) {
hits.reserve(2000);
}
@ -63,9 +60,7 @@ template <typename T> class VarClusterFinder {
noiseMap = noise_map;
use_noise_map = true;
}
void set_peripheralThresholdFactor(int factor) {
peripheralThresholdFactor_ = factor;
}
void set_peripheralThresholdFactor(int factor) { peripheralThresholdFactor_ = factor; }
void find_clusters(NDView<T, 2> img);
void find_clusters_X(NDView<T, 2> img);
void rec_FillHit(int clusterIndex, int i, int j);
@ -149,8 +144,7 @@ template <typename T> int VarClusterFinder<T>::check_neighbours(int i, int j) {
}
}
template <typename T>
void VarClusterFinder<T>::find_clusters(NDView<T, 2> img) {
template <typename T> void VarClusterFinder<T>::find_clusters(NDView<T, 2> img) {
original_ = img;
labeled_ = 0;
peripheral_labeled_ = 0;
@ -162,8 +156,7 @@ void VarClusterFinder<T>::find_clusters(NDView<T, 2> img) {
store_clusters();
}
template <typename T>
void VarClusterFinder<T>::find_clusters_X(NDView<T, 2> img) {
template <typename T> void VarClusterFinder<T>::find_clusters_X(NDView<T, 2> img) {
original_ = img;
int clusterIndex = 0;
for (int i = 0; i < shape_[0]; ++i) {
@ -182,8 +175,7 @@ void VarClusterFinder<T>::find_clusters_X(NDView<T, 2> img) {
h_size.clear();
}
template <typename T>
void VarClusterFinder<T>::rec_FillHit(int clusterIndex, int i, int j) {
template <typename T> void VarClusterFinder<T>::rec_FillHit(int clusterIndex, int i, int j) {
// printf("original_(%d, %d)=%f\n", i, j, original_(i,j));
// printf("h_size[%d].size=%d\n", clusterIndex, h_size[clusterIndex].size);
if (h_size[clusterIndex].size < MAX_CLUSTER_SIZE) {
@ -211,15 +203,11 @@ void VarClusterFinder<T>::rec_FillHit(int clusterIndex, int i, int j) {
} else {
// if (h_size[clusterIndex].size < MAX_CLUSTER_SIZE){
// h_size[clusterIndex].size += 1;
// h_size[clusterIndex].rows[h_size[clusterIndex].size] =
// row; h_size[clusterIndex].cols[h_size[clusterIndex].size]
// = col;
// h_size[clusterIndex].enes[h_size[clusterIndex].size] =
// original_(row, col);
// h_size[clusterIndex].rows[h_size[clusterIndex].size] = row;
// h_size[clusterIndex].cols[h_size[clusterIndex].size] = col;
// h_size[clusterIndex].enes[h_size[clusterIndex].size] = original_(row, col);
// }// ? weather to include peripheral pixels
original_(row, col) =
0; // remove peripheral pixels, to avoid potential influence
// for pedestal updating
original_(row, col) = 0; // remove peripheral pixels, to avoid potential influence for pedestal updating
}
}
}
@ -287,8 +275,8 @@ template <typename T> void VarClusterFinder<T>::store_clusters() {
for (int i = 0; i < shape_[0]; ++i) {
for (int j = 0; j < shape_[1]; ++j) {
if (labeled_(i, j) != 0 || false
// (i-1 >= 0 and labeled_(i-1, j) != 0) or // another circle of
// peripheral pixels (j-1 >= 0 and labeled_(i, j-1) != 0) or
// (i-1 >= 0 and labeled_(i-1, j) != 0) or // another circle of peripheral pixels
// (j-1 >= 0 and labeled_(i, j-1) != 0) or
// (i+1 < shape_[0] and labeled_(i+1, j) != 0) or
// (j+1 < shape_[1] and labeled_(i, j+1) != 0)
) {

View File

@ -1,9 +1,9 @@
#pragma once
#include <aare/NDArray.hpp>
#include <algorithm>
#include <array>
#include <vector>
#include <aare/NDArray.hpp>
namespace aare {
/**
@ -15,24 +15,26 @@ namespace aare {
* @param last iterator to the last element
* @param val value to compare
* @return index of the last element that is smaller than val
*
*
*/
template <typename T>
size_t last_smaller(const T *first, const T *last, T val) {
for (auto iter = first + 1; iter != last; ++iter) {
size_t last_smaller(const T* first, const T* last, T val) {
for (auto iter = first+1; iter != last; ++iter) {
if (*iter >= val) {
return std::distance(first, iter - 1);
return std::distance(first, iter-1);
}
}
return std::distance(first, last - 1);
return std::distance(first, last-1);
}
template <typename T> size_t last_smaller(const NDArray<T, 1> &arr, T val) {
template <typename T>
size_t last_smaller(const NDArray<T, 1>& arr, T val) {
return last_smaller(arr.begin(), arr.end(), val);
}
template <typename T> size_t last_smaller(const std::vector<T> &vec, T val) {
return last_smaller(vec.data(), vec.data() + vec.size(), val);
template <typename T>
size_t last_smaller(const std::vector<T>& vec, T val) {
return last_smaller(vec.data(), vec.data()+vec.size(), val);
}
/**
@ -46,67 +48,64 @@ template <typename T> size_t last_smaller(const std::vector<T> &vec, T val) {
* @return index of the first element that is larger than val
*/
template <typename T>
size_t first_larger(const T *first, const T *last, T val) {
size_t first_larger(const T* first, const T* last, T val) {
for (auto iter = first; iter != last; ++iter) {
if (*iter > val) {
return std::distance(first, iter);
}
}
return std::distance(first, last - 1);
return std::distance(first, last-1);
}
template <typename T> size_t first_larger(const NDArray<T, 1> &arr, T val) {
template <typename T>
size_t first_larger(const NDArray<T, 1>& arr, T val) {
return first_larger(arr.begin(), arr.end(), val);
}
template <typename T> size_t first_larger(const std::vector<T> &vec, T val) {
return first_larger(vec.data(), vec.data() + vec.size(), val);
template <typename T>
size_t first_larger(const std::vector<T>& vec, T val) {
return first_larger(vec.data(), vec.data()+vec.size(), val);
}
/**
* @brief Index of the nearest element to val.
* Requires a sorted array. If there is no difference it takes the first
* element.
* Requires a sorted array. If there is no difference it takes the first element.
* @param first iterator to the first element
* @param last iterator to the last element
* @param val value to compare
* @return index of the nearest element
*/
template <typename T>
size_t nearest_index(const T *first, const T *last, T val) {
auto iter = std::min_element(first, last, [val](T a, T b) {
size_t nearest_index(const T* first, const T* last, T val) {
auto iter = std::min_element(first, last,
[val](T a, T b) {
return std::abs(a - val) < std::abs(b - val);
});
return std::distance(first, iter);
}
template <typename T> size_t nearest_index(const NDArray<T, 1> &arr, T val) {
template <typename T>
size_t nearest_index(const NDArray<T, 1>& arr, T val) {
return nearest_index(arr.begin(), arr.end(), val);
}
template <typename T> size_t nearest_index(const std::vector<T> &vec, T val) {
return nearest_index(vec.data(), vec.data() + vec.size(), val);
template <typename T>
size_t nearest_index(const std::vector<T>& vec, T val) {
return nearest_index(vec.data(), vec.data()+vec.size(), val);
}
template <typename T, size_t N>
size_t nearest_index(const std::array<T, N> &arr, T val) {
return nearest_index(arr.data(), arr.data() + arr.size(), val);
size_t nearest_index(const std::array<T,N>& arr, T val) {
return nearest_index(arr.data(), arr.data()+arr.size(), val);
}
template <typename T> std::vector<T> cumsum(const std::vector<T> &vec) {
template <typename T>
std::vector<T> cumsum(const std::vector<T>& vec) {
std::vector<T> result(vec.size());
std::partial_sum(vec.begin(), vec.end(), result.begin());
return result;
}
template <typename Container> bool all_equal(const Container &c) {
if (!c.empty() &&
std::all_of(begin(c), end(c),
[c](const typename Container::value_type &element) {
return element == c.front();
}))
return true;
return false;
}
} // namespace aare

View File

@ -1,27 +1,26 @@
#pragma once
#include <aare/NDView.hpp>
#include <cstdint>
#include <vector>
#include <aare/NDView.hpp>
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,
NDView<uint16_t, 2> output);
void adc_sar_04_decode64to16(NDView<uint64_t, 2> input,
NDView<uint16_t, 2> output);
void adc_sar_05_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> output);
void adc_sar_04_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> output);
/**
* @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.
* @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.
* @throws std::out_of_range if weights.size() < 16
* @param input 16-bit input value
* @param weights vector of weights, size must be less than or equal to 16
*/
double apply_custom_weights(uint16_t input, const NDView<double, 1> weights);
void apply_custom_weights(NDView<uint16_t, 1> input, NDView<double, 1> output,
const NDView<double, 1> weights);
void apply_custom_weights(NDView<uint16_t, 1> input, NDView<double, 1> output, const NDView<double, 1> weights);
} // namespace aare

View File

@ -9,6 +9,7 @@
#include <stdexcept>
#include <string>
#include <string_view>
#include <type_traits>
#include <variant>
#include <vector>
@ -199,8 +200,6 @@ struct DetectorGeometry {
int module_gap_row{};
int module_gap_col{};
std::vector<ModuleGeometry> module_pixel_0;
auto size() const { return module_pixel_0.size(); }
};
struct ROI {
@ -260,4 +259,12 @@ template <> FrameDiscardPolicy StringTo(const std::string & /*mode*/);
using DataTypeVariants = std::variant<uint16_t, uint32_t>;
template <typename T, bool = std::is_integral_v<T>> struct make_signed {
using type = T;
};
template <typename T> struct make_signed<T, true> {
using type = std::make_signed_t<T>;
};
} // namespace aare

View File

@ -1,15 +1,16 @@
#pragma once
#include "aare/RawMasterFile.hpp" //ROI refactor away
#include "aare/defs.hpp"
namespace aare {
#include "aare/RawMasterFile.hpp" //ROI refactor away
namespace aare{
/**
* @brief Update the detector geometry given a region of interest
*
* @param geo
* @param roi
* @return DetectorGeometry
*
* @param geo
* @param roi
* @return DetectorGeometry
*/
DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, ROI roi);
} // namespace aare

View File

@ -1,141 +0,0 @@
#pragma once
/*Utility to log to console*/
#include <iostream>
#include <sstream>
#include <sys/time.h>
namespace aare {
#define RED "\x1b[31m"
#define GREEN "\x1b[32m"
#define YELLOW "\x1b[33m"
#define BLUE "\x1b[34m"
#define MAGENTA "\x1b[35m"
#define CYAN "\x1b[36m"
#define GRAY "\x1b[37m"
#define DARKGRAY "\x1b[30m"
#define BG_BLACK "\x1b[48;5;232m"
#define BG_RED "\x1b[41m"
#define BG_GREEN "\x1b[42m"
#define BG_YELLOW "\x1b[43m"
#define BG_BLUE "\x1b[44m"
#define BG_MAGENTA "\x1b[45m"
#define BG_CYAN "\x1b[46m"
#define RESET "\x1b[0m"
#define BOLD "\x1b[1m"
enum TLogLevel {
logERROR,
logWARNING,
logINFOBLUE,
logINFOGREEN,
logINFORED,
logINFOCYAN,
logINFOMAGENTA,
logINFO,
logDEBUG, // constructors, destructors etc. should still give too much
// output
logDEBUG1,
logDEBUG2,
logDEBUG3,
logDEBUG4,
logDEBUG5
};
// Compiler should optimize away anything below this value
#ifndef AARE_LOG_LEVEL
#define AARE_LOG_LEVEL \
"LOG LEVEL NOT SET IN CMAKE" // This is configured in the main
// CMakeLists.txt
#endif
#define __AT__ \
std::string(__FILE__) + std::string("::") + std::string(__func__) + \
std::string("(): ")
#define __SHORT_FORM_OF_FILE__ \
(strrchr(__FILE__, '/') ? strrchr(__FILE__, '/') + 1 : __FILE__)
#define __SHORT_AT__ \
std::string(__SHORT_FORM_OF_FILE__) + std::string("::") + \
std::string(__func__) + std::string("(): ")
class Logger {
std::ostringstream os;
TLogLevel m_level = AARE_LOG_LEVEL;
public:
Logger() = default;
explicit Logger(TLogLevel level) : m_level(level){};
~Logger() {
// output in the destructor to allow for << syntax
os << RESET << '\n';
std::clog << os.str() << std::flush; // Single write
}
static TLogLevel &
ReportingLevel() { // singelton eeh TODO! Do we need a runtime option?
static TLogLevel reportingLevel = logDEBUG5;
return reportingLevel;
}
// Danger this buffer need as many elements as TLogLevel
static const char *Color(TLogLevel level) noexcept {
static const char *const colors[] = {
RED BOLD, YELLOW BOLD, BLUE, GREEN, RED, CYAN, MAGENTA,
RESET, RESET, RESET, RESET, RESET, RESET, RESET};
// out of bounds
if (level < 0 || level >= sizeof(colors) / sizeof(colors[0])) {
return RESET;
}
return colors[level];
}
// Danger this buffer need as many elements as TLogLevel
static std::string ToString(TLogLevel level) {
static const char *const buffer[] = {
"ERROR", "WARNING", "INFO", "INFO", "INFO",
"INFO", "INFO", "INFO", "DEBUG", "DEBUG1",
"DEBUG2", "DEBUG3", "DEBUG4", "DEBUG5"};
// out of bounds
if (level < 0 || level >= sizeof(buffer) / sizeof(buffer[0])) {
return "UNKNOWN";
}
return buffer[level];
}
std::ostringstream &Get() {
os << Color(m_level) << "- " << Timestamp() << " " << ToString(m_level)
<< ": ";
return os;
}
static std::string Timestamp() {
constexpr size_t buffer_len = 12;
char buffer[buffer_len];
time_t t;
::time(&t);
tm r;
strftime(buffer, buffer_len, "%X", localtime_r(&t, &r));
buffer[buffer_len - 1] = '\0';
struct timeval tv;
gettimeofday(&tv, nullptr);
constexpr size_t result_len = 100;
char result[result_len];
snprintf(result, result_len, "%s.%03ld", buffer,
static_cast<long>(tv.tv_usec) / 1000);
result[result_len - 1] = '\0';
return result;
}
};
// TODO! Do we need to keep the runtime option?
#define LOG(level) \
if (level > AARE_LOG_LEVEL) \
; \
else if (level > aare::Logger::ReportingLevel()) \
; \
else \
aare::Logger(level).Get()
} // namespace aare

View File

@ -4,9 +4,9 @@
#include <string>
namespace aare {
/**
/**
* @brief Get the error message from an ifstream object
*/
*/
std::string ifstream_error_msg(std::ifstream &ifs);
} // namespace aare

View File

@ -1,18 +1,18 @@
#include <thread>
#include <utility>
#include <vector>
#include <utility>
namespace aare {
template <typename F>
void RunInParallel(F func, const std::vector<std::pair<int, int>> &tasks) {
// auto tasks = split_task(0, y.shape(0), n_threads);
std::vector<std::thread> threads;
for (auto &task : tasks) {
threads.push_back(std::thread(func, task.first, task.second));
template<typename F>
void RunInParallel(F func, const std::vector<std::pair<int, int>>& tasks) {
// auto tasks = split_task(0, y.shape(0), n_threads);
std::vector<std::thread> threads;
for (auto &task : tasks) {
threads.push_back(std::thread(func, task.first, task.second));
}
for (auto &thread : threads) {
thread.join();
}
}
for (auto &thread : threads) {
thread.join();
}
}
} // namespace aare

View File

@ -1,39 +1,22 @@
from . import _aare
from ._aare import ClusterFinder_Cluster3x3i, ClusterFinder_Cluster2x2i, ClusterFinderMT_Cluster3x3i, ClusterFinderMT_Cluster2x2i, ClusterCollector_Cluster3x3i, ClusterCollector_Cluster2x2i
from ._aare import ClusterFileSink_Cluster3x3i, ClusterFileSink_Cluster2x2i
import numpy as np
_supported_cluster_sizes = [(2,2), (3,3), (5,5), (7,7), (9,9),]
def _type_to_char(dtype):
if dtype == np.int32:
return 'i'
elif dtype == np.float32:
return 'f'
elif dtype == np.float64:
return 'd'
else:
raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32, np.float32, and np.float64 are supported.")
def _get_class(name, cluster_size, dtype):
"""
Helper function to get the class based on the name, cluster size, and dtype.
"""
try:
class_name = f"{name}_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
def ClusterFinder(image_size, cluster_size, n_sigma=5, dtype = np.int32, capacity = 1024):
"""
Factory function to create a ClusterFinder object. Provides a cleaner syntax for
the templated ClusterFinder in C++.
"""
cls = _get_class("ClusterFinder", cluster_size, dtype)
return cls(image_size, n_sigma=n_sigma, capacity=capacity)
if dtype == np.int32 and cluster_size == (3,3):
return ClusterFinder_Cluster3x3i(image_size, n_sigma = n_sigma, capacity=capacity)
elif dtype == np.int32 and cluster_size == (2,2):
return ClusterFinder_Cluster2x2i(image_size, n_sigma = n_sigma, capacity=capacity)
else:
#TODO! add the other formats
raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.")
def ClusterFinderMT(image_size, cluster_size = (3,3), dtype=np.int32, n_sigma=5, capacity = 1024, n_threads = 3):
@ -42,9 +25,15 @@ def ClusterFinderMT(image_size, cluster_size = (3,3), dtype=np.int32, n_sigma=5,
the templated ClusterFinderMT in C++.
"""
cls = _get_class("ClusterFinderMT", cluster_size, dtype)
return cls(image_size, n_sigma=n_sigma, capacity=capacity, n_threads=n_threads)
if dtype == np.int32 and cluster_size == (3,3):
return ClusterFinderMT_Cluster3x3i(image_size, n_sigma = n_sigma,
capacity = capacity, n_threads = n_threads)
elif dtype == np.int32 and cluster_size == (2,2):
return ClusterFinderMT_Cluster2x2i(image_size, n_sigma = n_sigma,
capacity = capacity, n_threads = n_threads)
else:
#TODO! add the other formats
raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.")
def ClusterCollector(clusterfindermt, cluster_size = (3,3), dtype=np.int32):
@ -53,8 +42,14 @@ def ClusterCollector(clusterfindermt, cluster_size = (3,3), dtype=np.int32):
the templated ClusterCollector in C++.
"""
cls = _get_class("ClusterCollector", cluster_size, dtype)
return cls(clusterfindermt)
if dtype == np.int32 and cluster_size == (3,3):
return ClusterCollector_Cluster3x3i(clusterfindermt)
elif dtype == np.int32 and cluster_size == (2,2):
return ClusterCollector_Cluster2x2i(clusterfindermt)
else:
#TODO! add the other formats
raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.")
def ClusterFileSink(clusterfindermt, cluster_file, dtype=np.int32):
"""
@ -62,26 +57,11 @@ def ClusterFileSink(clusterfindermt, cluster_file, dtype=np.int32):
the templated ClusterCollector in C++.
"""
cls = _get_class("ClusterFileSink", clusterfindermt.cluster_size, dtype)
return cls(clusterfindermt, cluster_file)
def ClusterFile(fname, cluster_size=(3,3), dtype=np.int32, chunk_size = 1000):
"""
Factory function to create a ClusterFile object. Provides a cleaner syntax for
the templated ClusterFile in C++.
.. code-block:: python
from aare import ClusterFile
with ClusterFile("clusters.clust", cluster_size=(3,3), dtype=np.int32) as cf:
# cf is now a ClusterFile_Cluster3x3i object but you don't need to know that.
for clusters in cf:
# Loop over clusters in chunks of 1000
# The type of clusters will be a ClusterVector_Cluster3x3i in this case
"""
cls = _get_class("ClusterFile", cluster_size, dtype)
return cls(fname, chunk_size=chunk_size)
if dtype == np.int32 and clusterfindermt.cluster_size == (3,3):
return ClusterFileSink_Cluster3x3i(clusterfindermt, cluster_file)
elif dtype == np.int32 and clusterfindermt.cluster_size == (2,2):
return ClusterFileSink_Cluster2x2i(clusterfindermt, cluster_file)
else:
#TODO! add the other formats
raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.")

View File

@ -5,12 +5,13 @@ from . import _aare
from ._aare import File, RawMasterFile, RawSubFile, JungfrauDataFile
from ._aare import Pedestal_d, Pedestal_f, ClusterFinder_Cluster3x3i, VarClusterFinder
from ._aare import DetectorType
from ._aare import ClusterFile_Cluster3x3i as ClusterFile
from ._aare import hitmap
from ._aare import ROI
# from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i
from .ClusterFinder import ClusterFinder, ClusterCollector, ClusterFinderMT, ClusterFileSink, ClusterFile
from .ClusterFinder import ClusterFinder, ClusterCollector, ClusterFinderMT, ClusterFileSink
from .ClusterVector import ClusterVector

View File

@ -1,89 +1,79 @@
import sys
sys.path.append('/home/l_msdetect/erik/aare/build')
from aare._aare import ClusterVector_i, Interpolator
from aare import RawSubFile, DetectorType, RawFile
from pathlib import Path
path = Path("/home/l_msdetect/erik/data/aare-test-data/raw/jungfrau/")
f = RawSubFile(path/"jungfrau_single_d0_f0_0.raw", DetectorType.Jungfrau, 512, 1024, 16)
# f = RawFile(path/"jungfrau_single_master_0.json")
# from aare._aare import ClusterVector_i, Interpolator
# import pickle
# import numpy as np
# import matplotlib.pyplot as plt
# import boost_histogram as bh
# import torch
# import math
# import time
import pickle
import numpy as np
import matplotlib.pyplot as plt
import boost_histogram as bh
import torch
import math
import time
# def gaussian_2d(mx, my, sigma = 1, res=100, grid_size = 2):
# """
# Generate a 2D gaussian as position mx, my, with sigma=sigma.
# The gaussian is placed on a 2x2 pixel matrix with resolution
# res in one dimesion.
# """
# x = torch.linspace(0, pixel_size*grid_size, res)
# x,y = torch.meshgrid(x,x, indexing="ij")
# return 1 / (2*math.pi*sigma**2) * \
# torch.exp(-((x - my)**2 / (2*sigma**2) + (y - mx)**2 / (2*sigma**2)))
def gaussian_2d(mx, my, sigma = 1, res=100, grid_size = 2):
"""
Generate a 2D gaussian as position mx, my, with sigma=sigma.
The gaussian is placed on a 2x2 pixel matrix with resolution
res in one dimesion.
"""
x = torch.linspace(0, pixel_size*grid_size, res)
x,y = torch.meshgrid(x,x, indexing="ij")
return 1 / (2*math.pi*sigma**2) * \
torch.exp(-((x - my)**2 / (2*sigma**2) + (y - mx)**2 / (2*sigma**2)))
# scale = 1000 #Scale factor when converting to integer
# pixel_size = 25 #um
# grid = 2
# resolution = 100
# sigma_um = 10
# xa = np.linspace(0,grid*pixel_size,resolution)
# ticks = [0, 25, 50]
scale = 1000 #Scale factor when converting to integer
pixel_size = 25 #um
grid = 2
resolution = 100
sigma_um = 10
xa = np.linspace(0,grid*pixel_size,resolution)
ticks = [0, 25, 50]
# hit = np.array((20,20))
# etahist_fname = "/home/l_msdetect/erik/tmp/test_hist.pkl"
hit = np.array((20,20))
etahist_fname = "/home/l_msdetect/erik/tmp/test_hist.pkl"
# local_resolution = 99
# grid_size = 3
# xaxis = np.linspace(0,grid_size*pixel_size, local_resolution)
# t = gaussian_2d(hit[0],hit[1], grid_size = grid_size, sigma = 10, res = local_resolution)
# pixels = t.reshape(grid_size, t.shape[0] // grid_size, grid_size, t.shape[1] // grid_size).sum(axis = 3).sum(axis = 1)
# pixels = pixels.numpy()
# pixels = (pixels*scale).astype(np.int32)
# v = ClusterVector_i(3,3)
# v.push_back(1,1, pixels)
local_resolution = 99
grid_size = 3
xaxis = np.linspace(0,grid_size*pixel_size, local_resolution)
t = gaussian_2d(hit[0],hit[1], grid_size = grid_size, sigma = 10, res = local_resolution)
pixels = t.reshape(grid_size, t.shape[0] // grid_size, grid_size, t.shape[1] // grid_size).sum(axis = 3).sum(axis = 1)
pixels = pixels.numpy()
pixels = (pixels*scale).astype(np.int32)
v = ClusterVector_i(3,3)
v.push_back(1,1, pixels)
# with open(etahist_fname, "rb") as f:
# hist = pickle.load(f)
# eta = hist.view().copy()
# etabinsx = np.array(hist.axes.edges.T[0].flat)
# etabinsy = np.array(hist.axes.edges.T[1].flat)
# ebins = np.array(hist.axes.edges.T[2].flat)
# p = Interpolator(eta, etabinsx[0:-1], etabinsy[0:-1], ebins[0:-1])
with open(etahist_fname, "rb") as f:
hist = pickle.load(f)
eta = hist.view().copy()
etabinsx = np.array(hist.axes.edges.T[0].flat)
etabinsy = np.array(hist.axes.edges.T[1].flat)
ebins = np.array(hist.axes.edges.T[2].flat)
p = Interpolator(eta, etabinsx[0:-1], etabinsy[0:-1], ebins[0:-1])
# #Generate the hit
#Generate the hit
# tmp = p.interpolate(v)
# print(f'tmp:{tmp}')
# pos = np.array((tmp['x'], tmp['y']))*25
tmp = p.interpolate(v)
print(f'tmp:{tmp}')
pos = np.array((tmp['x'], tmp['y']))*25
# print(pixels)
# fig, ax = plt.subplots(figsize = (7,7))
# ax.pcolormesh(xaxis, xaxis, t)
# ax.plot(*pos, 'o')
# ax.set_xticks([0,25,50,75])
# ax.set_yticks([0,25,50,75])
# ax.set_xlim(0,75)
# ax.set_ylim(0,75)
# ax.grid()
# print(f'{hit=}')
# print(f'{pos=}')
print(pixels)
fig, ax = plt.subplots(figsize = (7,7))
ax.pcolormesh(xaxis, xaxis, t)
ax.plot(*pos, 'o')
ax.set_xticks([0,25,50,75])
ax.set_yticks([0,25,50,75])
ax.set_xlim(0,75)
ax.set_ylim(0,75)
ax.grid()
print(f'{hit=}')
print(f'{pos=}')

View File

@ -1,64 +0,0 @@
#include "aare/Cluster.hpp"
#include <cstdint>
#include <filesystem>
#include <fmt/format.h>
#include <pybind11/numpy.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl_bind.h>
namespace py = pybind11;
using pd_type = double;
using namespace aare;
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
template <typename Type, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
void define_Cluster(py::module &m, const std::string &typestr) {
auto class_name = fmt::format("Cluster{}", typestr);
py::class_<Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>>(
m, class_name.c_str(), py::buffer_protocol())
.def(py::init([](uint8_t x, uint8_t y, py::array_t<Type> data) {
py::buffer_info buf_info = data.request();
Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType> cluster;
cluster.x = x;
cluster.y = y;
auto r = data.template unchecked<1>(); // no bounds checks
for (py::ssize_t i = 0; i < data.size(); ++i) {
cluster.data[i] = r(i);
}
return cluster;
}));
/*
//TODO! Review if to keep or not
.def_property(
"data",
[](ClusterType &c) -> py::array {
return py::array(py::buffer_info(
c.data, sizeof(Type),
py::format_descriptor<Type>::format(), // Type
// format
1, // Number of dimensions
{static_cast<ssize_t>(ClusterSizeX *
ClusterSizeY)}, // Shape (flattened)
{sizeof(Type)} // Stride (step size between elements)
));
},
[](ClusterType &c, py::array_t<Type> arr) {
py::buffer_info buf_info = arr.request();
Type *ptr = static_cast<Type *>(buf_info.ptr);
std::copy(ptr, ptr + ClusterSizeX * ClusterSizeY,
c.data); // TODO dont iterate over centers!!!
});
*/
}
#pragma GCC diagnostic pop

View File

@ -1,44 +0,0 @@
#include "aare/ClusterCollector.hpp"
#include "aare/ClusterFileSink.hpp"
#include "aare/ClusterFinder.hpp"
#include "aare/ClusterFinderMT.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/NDView.hpp"
#include "aare/Pedestal.hpp"
#include "np_helper.hpp"
#include <cstdint>
#include <filesystem>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl_bind.h>
namespace py = pybind11;
using pd_type = double;
using namespace aare;
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
void define_ClusterCollector(py::module &m, const std::string &typestr) {
auto class_name = fmt::format("ClusterCollector_{}", typestr);
using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
py::class_<ClusterCollector<ClusterType>>(m, class_name.c_str())
.def(py::init<ClusterFinderMT<ClusterType, uint16_t, double> *>())
.def("stop", &ClusterCollector<ClusterType>::stop)
.def(
"steal_clusters",
[](ClusterCollector<ClusterType> &self) {
auto v = new std::vector<ClusterVector<ClusterType>>(
self.steal_clusters());
return v; // TODO change!!!
},
py::return_value_policy::take_ownership);
}
#pragma GCC diagnostic pop

View File

@ -1,37 +0,0 @@
#include "aare/ClusterCollector.hpp"
#include "aare/ClusterFileSink.hpp"
#include "aare/ClusterFinder.hpp"
#include "aare/ClusterFinderMT.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/NDView.hpp"
#include "aare/Pedestal.hpp"
#include "np_helper.hpp"
#include <cstdint>
#include <filesystem>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl_bind.h>
namespace py = pybind11;
using pd_type = double;
using namespace aare;
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
void define_ClusterFileSink(py::module &m, const std::string &typestr) {
auto class_name = fmt::format("ClusterFileSink_{}", typestr);
using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
py::class_<ClusterFileSink<ClusterType>>(m, class_name.c_str())
.def(py::init<ClusterFinderMT<ClusterType, uint16_t, double> *,
const std::filesystem::path &>())
.def("stop", &ClusterFileSink<ClusterType>::stop);
}
#pragma GCC diagnostic pop

View File

@ -1,77 +0,0 @@
#include "aare/ClusterCollector.hpp"
#include "aare/ClusterFileSink.hpp"
#include "aare/ClusterFinder.hpp"
#include "aare/ClusterFinderMT.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/NDView.hpp"
#include "aare/Pedestal.hpp"
#include "np_helper.hpp"
#include <cstdint>
#include <filesystem>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl_bind.h>
namespace py = pybind11;
using pd_type = double;
using namespace aare;
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
void define_ClusterFinder(py::module &m, const std::string &typestr) {
auto class_name = fmt::format("ClusterFinder_{}", typestr);
using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
py::class_<ClusterFinder<ClusterType, uint16_t, pd_type>>(
m, class_name.c_str())
.def(py::init<Shape<2>, pd_type, size_t>(), py::arg("image_size"),
py::arg("n_sigma") = 5.0, py::arg("capacity") = 1'000'000)
.def("push_pedestal_frame",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
py::array_t<uint16_t> frame) {
auto view = make_view_2d(frame);
self.push_pedestal_frame(view);
})
.def("clear_pedestal",
&ClusterFinder<ClusterType, uint16_t, pd_type>::clear_pedestal)
.def_property_readonly(
"pedestal",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self) {
auto pd = new NDArray<pd_type, 2>{};
*pd = self.pedestal();
return return_image_data(pd);
})
.def_property_readonly(
"noise",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self) {
auto arr = new NDArray<pd_type, 2>{};
*arr = self.noise();
return return_image_data(arr);
})
.def(
"steal_clusters",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
bool realloc_same_capacity) {
ClusterVector<ClusterType> clusters =
self.steal_clusters(realloc_same_capacity);
return clusters;
},
py::arg("realloc_same_capacity") = false)
.def(
"find_clusters",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
py::array_t<uint16_t> frame, uint64_t frame_number) {
auto view = make_view_2d(frame);
self.find_clusters(view, frame_number);
return;
},
py::arg(), py::arg("frame_number") = 0);
}
#pragma GCC diagnostic pop

View File

@ -1,81 +0,0 @@
#include "aare/ClusterCollector.hpp"
#include "aare/ClusterFileSink.hpp"
#include "aare/ClusterFinder.hpp"
#include "aare/ClusterFinderMT.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/NDView.hpp"
#include "aare/Pedestal.hpp"
#include "np_helper.hpp"
#include <cstdint>
#include <filesystem>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl_bind.h>
namespace py = pybind11;
using pd_type = double;
using namespace aare;
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
void define_ClusterFinderMT(py::module &m, const std::string &typestr) {
auto class_name = fmt::format("ClusterFinderMT_{}", typestr);
using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
py::class_<ClusterFinderMT<ClusterType, uint16_t, pd_type>>(
m, class_name.c_str())
.def(py::init<Shape<2>, pd_type, size_t, size_t>(),
py::arg("image_size"), py::arg("n_sigma") = 5.0,
py::arg("capacity") = 2048, py::arg("n_threads") = 3)
.def("push_pedestal_frame",
[](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
py::array_t<uint16_t> frame) {
auto view = make_view_2d(frame);
self.push_pedestal_frame(view);
})
.def(
"find_clusters",
[](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
py::array_t<uint16_t> frame, uint64_t frame_number) {
auto view = make_view_2d(frame);
self.find_clusters(view, frame_number);
return;
},
py::arg(), py::arg("frame_number") = 0)
.def_property_readonly(
"cluster_size",
[](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self) {
return py::make_tuple(ClusterSizeX, ClusterSizeY);
})
.def("clear_pedestal",
&ClusterFinderMT<ClusterType, uint16_t, pd_type>::clear_pedestal)
.def("sync", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::sync)
.def("stop", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::stop)
.def("start", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::start)
.def(
"pedestal",
[](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
size_t thread_index) {
auto pd = new NDArray<pd_type, 2>{};
*pd = self.pedestal(thread_index);
return return_image_data(pd);
},
py::arg("thread_index") = 0)
.def(
"noise",
[](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
size_t thread_index) {
auto arr = new NDArray<pd_type, 2>{};
*arr = self.noise(thread_index);
return return_image_data(arr);
},
py::arg("thread_index") = 0);
}
#pragma GCC diagnostic pop

View File

@ -44,11 +44,10 @@ void define_ClusterVector(py::module &m, const std::string &typestr) {
auto *vec = new std::vector<Type>(self.sum());
return return_vector(vec);
})
.def("sum_2x2",
[](ClusterVector<ClusterType> &self) {
auto *vec = new std::vector<Type>(self.sum_2x2());
return return_vector(vec);
})
.def("sum_2x2", [](ClusterVector<ClusterType> &self){
auto *vec = new std::vector<Type>(self.sum_2x2());
return return_vector(vec);
})
.def_property_readonly("size", &ClusterVector<ClusterType>::size)
.def("item_size", &ClusterVector<ClusterType>::item_size)
.def_property_readonly("fmt",
@ -102,6 +101,4 @@ void define_ClusterVector(py::module &m, const std::string &typestr) {
return hitmap;
});
}
#pragma GCC diagnostic pop
}

211
python/src/cluster.hpp Normal file
View File

@ -0,0 +1,211 @@
#include "aare/ClusterCollector.hpp"
#include "aare/ClusterFileSink.hpp"
#include "aare/ClusterFinder.hpp"
#include "aare/ClusterFinderMT.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/NDView.hpp"
#include "aare/Pedestal.hpp"
#include "np_helper.hpp"
#include <cstdint>
#include <filesystem>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl_bind.h>
namespace py = pybind11;
using pd_type = double;
using namespace aare;
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
template <typename Type, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
void define_cluster(py::module &m, const std::string &typestr) {
auto class_name = fmt::format("Cluster{}", typestr);
py::class_<Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>>(
m, class_name.c_str(), py::buffer_protocol())
.def(py::init([](uint8_t x, uint8_t y, py::array_t<Type> data) {
py::buffer_info buf_info = data.request();
Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType> cluster;
cluster.x = x;
cluster.y = y;
auto r = data.template unchecked<1>(); // no bounds checks
for (py::ssize_t i = 0; i < data.size(); ++i) {
cluster.data[i] = r(i);
}
return cluster;
}));
/*
.def_property(
"data",
[](ClusterType &c) -> py::array {
return py::array(py::buffer_info(
c.data, sizeof(Type),
py::format_descriptor<Type>::format(), // Type
// format
1, // Number of dimensions
{static_cast<ssize_t>(ClusterSizeX *
ClusterSizeY)}, // Shape (flattened)
{sizeof(Type)} // Stride (step size between elements)
));
},
[](ClusterType &c, py::array_t<Type> arr) {
py::buffer_info buf_info = arr.request();
Type *ptr = static_cast<Type *>(buf_info.ptr);
std::copy(ptr, ptr + ClusterSizeX * ClusterSizeY,
c.data); // TODO dont iterate over centers!!!
});
*/
}
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
void define_cluster_finder_mt_bindings(py::module &m,
const std::string &typestr) {
auto class_name = fmt::format("ClusterFinderMT_{}", typestr);
using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
py::class_<ClusterFinderMT<ClusterType, uint16_t, pd_type>>(
m, class_name.c_str())
.def(py::init<Shape<2>, pd_type, size_t, size_t>(),
py::arg("image_size"), py::arg("n_sigma") = 5.0,
py::arg("capacity") = 2048, py::arg("n_threads") = 3)
.def("push_pedestal_frame",
[](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
py::array_t<uint16_t> frame) {
auto view = make_view_2d(frame);
self.push_pedestal_frame(view);
})
.def(
"find_clusters",
[](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
py::array_t<uint16_t> frame, uint64_t frame_number) {
auto view = make_view_2d(frame);
self.find_clusters(view, frame_number);
return;
},
py::arg(), py::arg("frame_number") = 0)
.def_property_readonly("cluster_size", [](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self){
return py::make_tuple(ClusterSizeX, ClusterSizeY);
})
.def("clear_pedestal",
&ClusterFinderMT<ClusterType, uint16_t, pd_type>::clear_pedestal)
.def("sync", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::sync)
.def("stop", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::stop)
.def("start", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::start)
.def(
"pedestal",
[](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
size_t thread_index) {
auto pd = new NDArray<pd_type, 2>{};
*pd = self.pedestal(thread_index);
return return_image_data(pd);
},
py::arg("thread_index") = 0)
.def(
"noise",
[](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
size_t thread_index) {
auto arr = new NDArray<pd_type, 2>{};
*arr = self.noise(thread_index);
return return_image_data(arr);
},
py::arg("thread_index") = 0);
}
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
void define_cluster_collector_bindings(py::module &m,
const std::string &typestr) {
auto class_name = fmt::format("ClusterCollector_{}", typestr);
using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
py::class_<ClusterCollector<ClusterType>>(m, class_name.c_str())
.def(py::init<ClusterFinderMT<ClusterType, uint16_t, double> *>())
.def("stop", &ClusterCollector<ClusterType>::stop)
.def(
"steal_clusters",
[](ClusterCollector<ClusterType> &self) {
auto v = new std::vector<ClusterVector<ClusterType>>(
self.steal_clusters());
return v; // TODO change!!!
},
py::return_value_policy::take_ownership);
}
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
void define_cluster_file_sink_bindings(py::module &m,
const std::string &typestr) {
auto class_name = fmt::format("ClusterFileSink_{}", typestr);
using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
py::class_<ClusterFileSink<ClusterType>>(m, class_name.c_str())
.def(py::init<ClusterFinderMT<ClusterType, uint16_t, double> *,
const std::filesystem::path &>())
.def("stop", &ClusterFileSink<ClusterType>::stop);
}
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
void define_cluster_finder_bindings(py::module &m, const std::string &typestr) {
auto class_name = fmt::format("ClusterFinder_{}", typestr);
using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
py::class_<ClusterFinder<ClusterType, uint16_t, pd_type>>(
m, class_name.c_str())
.def(py::init<Shape<2>, pd_type, size_t>(), py::arg("image_size"),
py::arg("n_sigma") = 5.0, py::arg("capacity") = 1'000'000)
.def("push_pedestal_frame",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
py::array_t<uint16_t> frame) {
auto view = make_view_2d(frame);
self.push_pedestal_frame(view);
})
.def("clear_pedestal",
&ClusterFinder<ClusterType, uint16_t, pd_type>::clear_pedestal)
.def_property_readonly(
"pedestal",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self) {
auto pd = new NDArray<pd_type, 2>{};
*pd = self.pedestal();
return return_image_data(pd);
})
.def_property_readonly(
"noise",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self) {
auto arr = new NDArray<pd_type, 2>{};
*arr = self.noise();
return return_image_data(arr);
})
.def(
"steal_clusters",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
bool realloc_same_capacity) {
ClusterVector<ClusterType> clusters =
self.steal_clusters(realloc_same_capacity);
return clusters;
},
py::arg("realloc_same_capacity") = false)
.def(
"find_clusters",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
py::array_t<uint16_t> frame, uint64_t frame_number) {
auto view = make_view_2d(frame);
self.find_clusters(view, frame_number);
return;
},
py::arg(), py::arg("frame_number") = 0);
}
#pragma GCC diagnostic pop

View File

@ -21,7 +21,8 @@ using namespace ::aare;
template <typename Type, uint8_t CoordSizeX, uint8_t CoordSizeY,
typename CoordType = uint16_t>
void define_ClusterFile(py::module &m, const std::string &typestr) {
void define_cluster_file_io_bindings(py::module &m,
const std::string &typestr) {
using ClusterType = Cluster<Type, CoordSizeX, CoordSizeY, CoordType>;
@ -38,20 +39,19 @@ void define_ClusterFile(py::module &m, const std::string &typestr) {
self.read_clusters(n_clusters));
return v;
},
py::return_value_policy::take_ownership, py::arg("n_clusters"))
py::return_value_policy::take_ownership)
.def("read_frame",
[](ClusterFile<ClusterType> &self) {
auto v = new ClusterVector<ClusterType>(self.read_frame());
return v;
})
.def("set_roi", &ClusterFile<ClusterType>::set_roi,
py::arg("roi"))
.def("set_roi", &ClusterFile<ClusterType>::set_roi)
.def(
"set_noise_map",
[](ClusterFile<ClusterType> &self, py::array_t<int32_t> noise_map) {
auto view = make_view_2d(noise_map);
self.set_noise_map(view);
}, py::arg("noise_map"))
})
.def("set_gain_map",
[](ClusterFile<ClusterType> &self, py::array_t<double> gain_map) {

View File

@ -6,8 +6,8 @@
#include "aare/RawMasterFile.hpp"
#include "aare/RawSubFile.hpp"
#include "aare/decode.hpp"
#include "aare/defs.hpp"
#include "aare/decode.hpp"
// #include "aare/fClusterFileV2.hpp"
#include "np_helper.hpp"
@ -26,103 +26,95 @@ using namespace ::aare;
void define_ctb_raw_file_io_bindings(py::module &m) {
m.def("adc_sar_05_decode64to16", [](py::array_t<uint8_t> input) {
if (input.ndim() != 2) {
throw std::runtime_error(
"Only 2D arrays are supported at this moment");
}
m.def("adc_sar_05_decode64to16", [](py::array_t<uint8_t> input) {
// Create a 2D output array with the same shape as the input
std::vector<ssize_t> shape{input.shape(0),
input.shape(1) /
static_cast<ssize_t>(bits_per_byte)};
py::array_t<uint16_t> output(shape);
if(input.ndim() != 2){
throw std::runtime_error("Only 2D arrays are supported at this moment");
}
// Create a view of the input and output arrays
NDView<uint64_t, 2> input_view(
reinterpret_cast<uint64_t *>(input.mutable_data()),
{output.shape(0), output.shape(1)});
NDView<uint16_t, 2> output_view(output.mutable_data(),
{output.shape(0), output.shape(1)});
//Create a 2D output array with the same shape as the input
std::vector<ssize_t> shape{input.shape(0), input.shape(1)/static_cast<ssize_t>(bits_per_byte)};
py::array_t<uint16_t> output(shape);
adc_sar_05_decode64to16(input_view, output_view);
//Create a view of the input and output arrays
NDView<uint64_t, 2> input_view(reinterpret_cast<uint64_t*>(input.mutable_data()), {output.shape(0), output.shape(1)});
NDView<uint16_t, 2> output_view(output.mutable_data(), {output.shape(0), output.shape(1)});
adc_sar_05_decode64to16(input_view, output_view);
return output;
});
m.def("adc_sar_04_decode64to16", [](py::array_t<uint8_t> input) {
if(input.ndim() != 2){
throw std::runtime_error("Only 2D arrays are supported at this moment");
}
//Create a 2D output array with the same shape as the input
std::vector<ssize_t> shape{input.shape(0), input.shape(1)/static_cast<ssize_t>(bits_per_byte)};
py::array_t<uint16_t> output(shape);
//Create a view of the input and output arrays
NDView<uint64_t, 2> input_view(reinterpret_cast<uint64_t*>(input.mutable_data()), {output.shape(0), output.shape(1)});
NDView<uint16_t, 2> output_view(output.mutable_data(), {output.shape(0), output.shape(1)});
adc_sar_04_decode64to16(input_view, output_view);
return output;
});
m.def(
"apply_custom_weights",
[](py::array_t<uint16_t, py::array::c_style | py::array::forcecast> &input,
py::array_t<double, py::array::c_style | py::array::forcecast>
&weights) {
// Create new array with same shape as the input array (uninitialized values)
py::buffer_info buf = input.request();
py::array_t<double> output(buf.shape);
// Use NDViews to call into the C++ library
auto weights_view = make_view_1d(weights);
NDView<uint16_t, 1> input_view(input.mutable_data(), {input.size()});
NDView<double, 1> output_view(output.mutable_data(), {output.size()});
apply_custom_weights(input_view, output_view, weights_view);
return output;
});
m.def("adc_sar_04_decode64to16", [](py::array_t<uint8_t> input) {
if (input.ndim() != 2) {
throw std::runtime_error(
"Only 2D arrays are supported at this moment");
}
py::class_<CtbRawFile>(m, "CtbRawFile")
.def(py::init<const std::filesystem::path &>())
.def("read_frame",
[](CtbRawFile &self) {
size_t image_size = self.image_size_in_bytes();
py::array image;
std::vector<ssize_t> shape;
shape.reserve(2);
shape.push_back(1);
shape.push_back(image_size);
// Create a 2D output array with the same shape as the input
std::vector<ssize_t> shape{input.shape(0),
input.shape(1) /
static_cast<ssize_t>(bits_per_byte)};
py::array_t<uint16_t> output(shape);
py::array_t<DetectorHeader> header(1);
// Create a view of the input and output arrays
NDView<uint64_t, 2> input_view(
reinterpret_cast<uint64_t *>(input.mutable_data()),
{output.shape(0), output.shape(1)});
NDView<uint16_t, 2> output_view(output.mutable_data(),
{output.shape(0), output.shape(1)});
// always read bytes
image = py::array_t<uint8_t>(shape);
adc_sar_04_decode64to16(input_view, output_view);
self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()),
header.mutable_data());
return output;
});
return py::make_tuple(header, image);
})
.def("seek", &CtbRawFile::seek)
.def("tell", &CtbRawFile::tell)
.def("master", &CtbRawFile::master)
m.def("apply_custom_weights",
[](py::array_t<uint16_t, py::array::c_style | py::array::forcecast>
&input,
py::array_t<double, py::array::c_style | py::array::forcecast>
&weights) {
// Create new array with same shape as the input array
// (uninitialized values)
py::buffer_info buf = input.request();
py::array_t<double> output(buf.shape);
.def_property_readonly("image_size_in_bytes",
&CtbRawFile::image_size_in_bytes)
// Use NDViews to call into the C++ library
auto weights_view = make_view_1d(weights);
NDView<uint16_t, 1> input_view(input.mutable_data(),
{input.size()});
NDView<double, 1> output_view(output.mutable_data(),
{output.size()});
.def_property_readonly("frames_in_file", &CtbRawFile::frames_in_file);
apply_custom_weights(input_view, output_view, weights_view);
return output;
});
py::class_<CtbRawFile>(m, "CtbRawFile")
.def(py::init<const std::filesystem::path &>())
.def("read_frame",
[](CtbRawFile &self) {
size_t image_size = self.image_size_in_bytes();
py::array image;
std::vector<ssize_t> shape;
shape.reserve(2);
shape.push_back(1);
shape.push_back(image_size);
py::array_t<DetectorHeader> header(1);
// always read bytes
image = py::array_t<uint8_t>(shape);
self.read_into(
reinterpret_cast<std::byte *>(image.mutable_data()),
header.mutable_data());
return py::make_tuple(header, image);
})
.def("seek", &CtbRawFile::seek)
.def("tell", &CtbRawFile::tell)
.def("master", &CtbRawFile::master)
.def_property_readonly("image_size_in_bytes",
&CtbRawFile::image_size_in_bytes)
.def_property_readonly("frames_in_file", &CtbRawFile::frames_in_file);
}

View File

@ -20,13 +20,17 @@
namespace py = pybind11;
using namespace ::aare;
// Disable warnings for unused parameters, as we ignore some
// in the __exit__ method
//Disable warnings for unused parameters, as we ignore some
//in the __exit__ method
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
void define_file_io_bindings(py::module &m) {
py::enum_<DetectorType>(m, "DetectorType")
.value("Jungfrau", DetectorType::Jungfrau)
.value("Eiger", DetectorType::Eiger)
@ -37,10 +41,13 @@ void define_file_io_bindings(py::module &m) {
.value("ChipTestBoard", DetectorType::ChipTestBoard)
.value("Unknown", DetectorType::Unknown);
PYBIND11_NUMPY_DTYPE(DetectorHeader, frameNumber, expLength, packetNumber,
bunchId, timestamp, modId, row, column, reserved,
debug, roundRNumber, detType, version, packetMask);
py::class_<File>(m, "File")
.def(py::init([](const std::filesystem::path &fname) {
return File(fname, "r", {});
@ -105,18 +112,45 @@ void define_file_io_bindings(py::module &m) {
reinterpret_cast<std::byte *>(image.mutable_data()));
return image;
})
.def("read_n",
[](File &self, size_t n_frames) {
// adjust for actual frames left in the file
n_frames =
std::min(n_frames, self.total_frames() - self.tell());
if (n_frames == 0) {
throw std::runtime_error("No frames left in file");
}
std::vector<size_t> shape{n_frames, self.rows(), self.cols()};
.def("read_n", [](File &self, size_t n_frames) {
//adjust for actual frames left in the file
n_frames = std::min(n_frames, self.total_frames()-self.tell());
if(n_frames == 0){
throw std::runtime_error("No frames left in file");
}
std::vector<size_t> shape{n_frames, self.rows(), self.cols()};
py::array image;
const uint8_t item_size = self.bytes_per_pixel();
if (item_size == 1) {
image = py::array_t<uint8_t>(shape);
} else if (item_size == 2) {
image = py::array_t<uint16_t>(shape);
} else if (item_size == 4) {
image = py::array_t<uint32_t>(shape);
}
self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()),
n_frames);
return image;
})
.def("__enter__", [](File &self) { return &self; })
.def("__exit__",
[](File &self,
const std::optional<pybind11::type> &exc_type,
const std::optional<pybind11::object> &exc_value,
const std::optional<pybind11::object> &traceback) {
// self.close();
})
.def("__iter__", [](File &self) { return &self; })
.def("__next__", [](File &self) {
try{
const uint8_t item_size = self.bytes_per_pixel();
py::array image;
const uint8_t item_size = self.bytes_per_pixel();
std::vector<ssize_t> shape;
shape.reserve(2);
shape.push_back(self.rows());
shape.push_back(self.cols());
if (item_size == 1) {
image = py::array_t<uint8_t>(shape);
} else if (item_size == 2) {
@ -125,41 +159,14 @@ void define_file_io_bindings(py::module &m) {
image = py::array_t<uint32_t>(shape);
}
self.read_into(
reinterpret_cast<std::byte *>(image.mutable_data()),
n_frames);
reinterpret_cast<std::byte *>(image.mutable_data()));
return image;
})
.def("__enter__", [](File &self) { return &self; })
.def("__exit__",
[](File &self, const std::optional<pybind11::type> &exc_type,
const std::optional<pybind11::object> &exc_value,
const std::optional<pybind11::object> &traceback) {
// self.close();
})
.def("__iter__", [](File &self) { return &self; })
.def("__next__", [](File &self) {
try {
const uint8_t item_size = self.bytes_per_pixel();
py::array image;
std::vector<ssize_t> shape;
shape.reserve(2);
shape.push_back(self.rows());
shape.push_back(self.cols());
if (item_size == 1) {
image = py::array_t<uint8_t>(shape);
} else if (item_size == 2) {
image = py::array_t<uint16_t>(shape);
} else if (item_size == 4) {
image = py::array_t<uint32_t>(shape);
}
self.read_into(
reinterpret_cast<std::byte *>(image.mutable_data()));
return image;
} catch (std::runtime_error &e) {
}catch(std::runtime_error &e){
throw py::stop_iteration();
}
});
py::class_<FileConfig>(m, "FileConfig")
.def(py::init<>())
.def_readwrite("rows", &FileConfig::rows)
@ -176,6 +183,8 @@ void define_file_io_bindings(py::module &m) {
return "<FileConfig: " + a.to_string() + ">";
});
py::class_<ScanParameters>(m, "ScanParameters")
.def(py::init<const std::string &>())
.def(py::init<const ScanParameters &>())
@ -186,6 +195,7 @@ void define_file_io_bindings(py::module &m) {
.def_property_readonly("stop", &ScanParameters::stop)
.def_property_readonly("step", &ScanParameters::step);
py::class_<ROI>(m, "ROI")
.def(py::init<>())
.def(py::init<ssize_t, ssize_t, ssize_t, ssize_t>(), py::arg("xmin"),
@ -194,21 +204,23 @@ void define_file_io_bindings(py::module &m) {
.def_readwrite("xmax", &ROI::xmax)
.def_readwrite("ymin", &ROI::ymin)
.def_readwrite("ymax", &ROI::ymax)
.def("__str__",
[](const ROI &self) {
return fmt::format("ROI: xmin: {} xmax: {} ymin: {} ymax: {}",
self.xmin, self.xmax, self.ymin, self.ymax);
})
.def("__repr__",
[](const ROI &self) {
return fmt::format(
"<ROI: xmin: {} xmax: {} ymin: {} ymax: {}>", self.xmin,
self.xmax, self.ymin, self.ymax);
})
.def("__str__", [](const ROI& self){
return fmt::format("ROI: xmin: {} xmax: {} ymin: {} ymax: {}", self.xmin, self.xmax, self.ymin, self.ymax);
})
.def("__repr__", [](const ROI& self){
return fmt::format("<ROI: xmin: {} xmax: {} ymin: {} ymax: {}>", self.xmin, self.xmax, self.ymin, self.ymax);
})
.def("__iter__", [](const ROI &self) {
return py::make_iterator(&self.xmin, &self.ymax + 1); // NOLINT
return py::make_iterator(&self.xmin, &self.ymax+1); //NOLINT
});
#pragma GCC diagnostic pop
// py::class_<ClusterHeader>(m, "ClusterHeader")
// .def(py::init<>())

View File

@ -9,6 +9,7 @@
namespace py = pybind11;
using namespace pybind11::literals;
void define_fit_bindings(py::module &m) {
// TODO! Evaluate without converting to double
@ -60,8 +61,7 @@ void define_fit_bindings(py::module &m) {
py::array_t<double, py::array::c_style | py::array::forcecast> par) {
auto x_view = make_view_1d(x);
auto par_view = make_view_1d(par);
auto y =
new NDArray<double, 1>{aare::func::scurve(x_view, par_view)};
auto y = new NDArray<double, 1>{aare::func::scurve(x_view, par_view)};
return return_image_data(y);
},
R"(
@ -82,8 +82,7 @@ void define_fit_bindings(py::module &m) {
py::array_t<double, py::array::c_style | py::array::forcecast> par) {
auto x_view = make_view_1d(x);
auto par_view = make_view_1d(par);
auto y =
new NDArray<double, 1>{aare::func::scurve2(x_view, par_view)};
auto y = new NDArray<double, 1>{aare::func::scurve2(x_view, par_view)};
return return_image_data(y);
},
R"(
@ -140,6 +139,7 @@ n_threads : int, optional
py::array_t<double, py::array::c_style | py::array::forcecast> y,
py::array_t<double, py::array::c_style | py::array::forcecast> y_err,
int n_threads) {
if (y.ndim() == 3) {
// Allocate memory for the output
// Need to have pointers to allow python to manage
@ -173,6 +173,7 @@ n_threads : int, optional
auto y_view_err = make_view_1d(y_err);
auto x_view = make_view_1d(x);
double chi2 = 0;
aare::fit_gaus(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2);
@ -247,10 +248,11 @@ n_threads : int, optional
aare::fit_pol1(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2->view(), n_threads);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"par_err"_a = return_image_data(par_err),
"chi2"_a = return_image_data(chi2),
"Ndf"_a = y.shape(2) - 2);
} else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>({2});
auto par_err = new NDArray<double, 1>({2});
@ -287,7 +289,7 @@ n_threads : int, optional
)",
py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
//=========
//=========
m.def(
"fit_scurve",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
@ -331,12 +333,13 @@ n_threads : int, optional
auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
aare::fit_scurve(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2->view(), n_threads);
par_err->view(), chi2->view(), n_threads);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"par_err"_a = return_image_data(par_err),
"chi2"_a = return_image_data(chi2),
"Ndf"_a = y.shape(2) - 2);
} else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>({2});
auto par_err = new NDArray<double, 1>({2});
@ -348,7 +351,7 @@ n_threads : int, optional
double chi2 = 0;
aare::fit_scurve(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2);
par_err->view(), chi2);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"chi2"_a = chi2, "Ndf"_a = y.size() - 2);
@ -372,6 +375,7 @@ n_threads : int, optional
The number of threads to use. Default is 4.
)",
py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
m.def(
"fit_scurve2",
@ -416,12 +420,13 @@ n_threads : int, optional
auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
aare::fit_scurve2(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2->view(), n_threads);
par_err->view(), chi2->view(), n_threads);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"par_err"_a = return_image_data(par_err),
"chi2"_a = return_image_data(chi2),
"Ndf"_a = y.shape(2) - 2);
} else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>({6});
auto par_err = new NDArray<double, 1>({6});
@ -433,7 +438,7 @@ n_threads : int, optional
double chi2 = 0;
aare::fit_scurve2(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2);
par_err->view(), chi2);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"chi2"_a = chi2, "Ndf"_a = y.size() - 2);

View File

@ -21,7 +21,10 @@ using namespace ::aare;
auto read_dat_frame(JungfrauDataFile &self) {
py::array_t<JungfrauDataHeader> header(1);
py::array_t<uint16_t> image({self.rows(), self.cols()});
py::array_t<uint16_t> image({
self.rows(),
self.cols()
});
self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()),
header.mutable_data());
@ -37,7 +40,9 @@ auto read_n_dat_frames(JungfrauDataFile &self, size_t n_frames) {
}
py::array_t<JungfrauDataHeader> header(n_frames);
py::array_t<uint16_t> image({n_frames, self.rows(), self.cols()});
py::array_t<uint16_t> image({
n_frames, self.rows(),
self.cols()});
self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()),
n_frames, header.mutable_data());

View File

@ -1,26 +1,22 @@
// Files with bindings to the different classes
// New style file naming
#include "bind_Cluster.hpp"
#include "bind_ClusterCollector.hpp"
#include "bind_ClusterFile.hpp"
#include "bind_ClusterFileSink.hpp"
#include "bind_ClusterFinder.hpp"
#include "bind_ClusterFinderMT.hpp"
//New style file naming
#include "bind_ClusterVector.hpp"
// TODO! migrate the other names
//TODO! migrate the other names
#include "cluster.hpp"
#include "cluster_file.hpp"
#include "ctb_raw_file.hpp"
#include "file.hpp"
#include "fit.hpp"
#include "interpolation.hpp"
#include "jungfrau_data_file.hpp"
#include "pedestal.hpp"
#include "pixel_map.hpp"
#include "raw_file.hpp"
#include "raw_master_file.hpp"
#include "raw_sub_file.hpp"
#include "raw_master_file.hpp"
#include "raw_file.hpp"
#include "pixel_map.hpp"
#include "var_cluster.hpp"
#include "pedestal.hpp"
#include "jungfrau_data_file.hpp"
// Pybind stuff
#include <pybind11/pybind11.h>
@ -28,26 +24,6 @@
namespace py = pybind11;
/* MACRO that defines Cluster bindings for a specific size and type
T - Storage type of the cluster data (int, float, double)
N - Number of rows in the cluster
M - Number of columns in the cluster
U - Type of the pixel data (e.g., uint16_t)
TYPE_CODE - A character representing the type code (e.g., 'i' for int, 'd' for
double, 'f' for float)
*/
#define DEFINE_CLUSTER_BINDINGS(T, N, M, U, TYPE_CODE) \
define_ClusterFile<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE); \
define_ClusterVector<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE); \
define_ClusterFinder<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE); \
define_ClusterFinderMT<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE); \
define_ClusterFileSink<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE); \
define_ClusterCollector<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE); \
define_Cluster<T, N, M, U>(m, #N "x" #M #TYPE_CODE); \
register_calculate_eta<T, N, M, U>(m);
PYBIND11_MODULE(_aare, m) {
define_file_io_bindings(m);
define_raw_file_io_bindings(m);
@ -62,23 +38,59 @@ PYBIND11_MODULE(_aare, m) {
define_interpolation_bindings(m);
define_jungfrau_data_file_io_bindings(m);
DEFINE_CLUSTER_BINDINGS(int, 3, 3, uint16_t, i);
DEFINE_CLUSTER_BINDINGS(double, 3, 3, uint16_t, d);
DEFINE_CLUSTER_BINDINGS(float, 3, 3, uint16_t, f);
define_cluster_file_io_bindings<int, 3, 3, uint16_t>(m, "Cluster3x3i");
define_cluster_file_io_bindings<double, 3, 3, uint16_t>(m, "Cluster3x3d");
define_cluster_file_io_bindings<float, 3, 3, uint16_t>(m, "Cluster3x3f");
define_cluster_file_io_bindings<int, 2, 2, uint16_t>(m, "Cluster2x2i");
define_cluster_file_io_bindings<float, 2, 2, uint16_t>(m, "Cluster2x2f");
define_cluster_file_io_bindings<double, 2, 2, uint16_t>(m, "Cluster2x2d");
DEFINE_CLUSTER_BINDINGS(int, 2, 2, uint16_t, i);
DEFINE_CLUSTER_BINDINGS(double, 2, 2, uint16_t, d);
DEFINE_CLUSTER_BINDINGS(float, 2, 2, uint16_t, f);
define_ClusterVector<int, 3, 3, uint16_t>(m, "Cluster3x3i");
define_ClusterVector<double, 3, 3, uint16_t>(m, "Cluster3x3d");
define_ClusterVector<float, 3, 3, uint16_t>(m, "Cluster3x3f");
define_ClusterVector<int, 2, 2, uint16_t>(m, "Cluster2x2i");
define_ClusterVector<double, 2, 2, uint16_t>(m, "Cluster2x2d");
define_ClusterVector<float, 2, 2, uint16_t>(m, "Cluster2x2f");
DEFINE_CLUSTER_BINDINGS(int, 5, 5, uint16_t, i);
DEFINE_CLUSTER_BINDINGS(double, 5, 5, uint16_t, d);
DEFINE_CLUSTER_BINDINGS(float, 5, 5, uint16_t, f);
define_cluster_finder_bindings<int, 3, 3, uint16_t>(m, "Cluster3x3i");
define_cluster_finder_bindings<double, 3, 3, uint16_t>(m, "Cluster3x3d");
define_cluster_finder_bindings<float, 3, 3, uint16_t>(m, "Cluster3x3f");
define_cluster_finder_bindings<int, 2, 2, uint16_t>(m, "Cluster2x2i");
define_cluster_finder_bindings<double, 2, 2, uint16_t>(m, "Cluster2x2d");
define_cluster_finder_bindings<float, 2, 2, uint16_t>(m, "Cluster2x2f");
DEFINE_CLUSTER_BINDINGS(int, 7, 7, uint16_t, i);
DEFINE_CLUSTER_BINDINGS(double, 7, 7, uint16_t, d);
DEFINE_CLUSTER_BINDINGS(float, 7, 7, uint16_t, f);
define_cluster_finder_mt_bindings<int, 3, 3, uint16_t>(m, "Cluster3x3i");
define_cluster_finder_mt_bindings<double, 3, 3, uint16_t>(m, "Cluster3x3d");
define_cluster_finder_mt_bindings<float, 3, 3, uint16_t>(m, "Cluster3x3f");
define_cluster_finder_mt_bindings<int, 2, 2, uint16_t>(m, "Cluster2x2i");
define_cluster_finder_mt_bindings<double, 2, 2, uint16_t>(m, "Cluster2x2d");
define_cluster_finder_mt_bindings<float, 2, 2, uint16_t>(m, "Cluster2x2f");
DEFINE_CLUSTER_BINDINGS(int, 9, 9, uint16_t, i);
DEFINE_CLUSTER_BINDINGS(double, 9, 9, uint16_t, d);
DEFINE_CLUSTER_BINDINGS(float, 9, 9, uint16_t, f);
define_cluster_file_sink_bindings<int, 3, 3, uint16_t>(m, "Cluster3x3i");
define_cluster_file_sink_bindings<double, 3, 3, uint16_t>(m, "Cluster3x3d");
define_cluster_file_sink_bindings<float, 3, 3, uint16_t>(m, "Cluster3x3f");
define_cluster_file_sink_bindings<int, 2, 2, uint16_t>(m, "Cluster2x2i");
define_cluster_file_sink_bindings<double, 2, 2, uint16_t>(m, "Cluster2x2d");
define_cluster_file_sink_bindings<float, 2, 2, uint16_t>(m, "Cluster2x2f");
define_cluster_collector_bindings<int, 3, 3, uint16_t>(m, "Cluster3x3i");
define_cluster_collector_bindings<double, 3, 3, uint16_t>(m, "Cluster3x3f");
define_cluster_collector_bindings<float, 3, 3, uint16_t>(m, "Cluster3x3d");
define_cluster_collector_bindings<int, 2, 2, uint16_t>(m, "Cluster2x2i");
define_cluster_collector_bindings<double, 2, 2, uint16_t>(m, "Cluster2x2f");
define_cluster_collector_bindings<float, 2, 2, uint16_t>(m, "Cluster2x2d");
define_cluster<int, 3, 3, uint16_t>(m, "3x3i");
define_cluster<float, 3, 3, uint16_t>(m, "3x3f");
define_cluster<double, 3, 3, uint16_t>(m, "3x3d");
define_cluster<int, 2, 2, uint16_t>(m, "2x2i");
define_cluster<float, 2, 2, uint16_t>(m, "2x2f");
define_cluster<double, 2, 2, uint16_t>(m, "2x2d");
register_calculate_eta<int, 3, 3, uint16_t>(m);
register_calculate_eta<float, 3, 3, uint16_t>(m);
register_calculate_eta<double, 3, 3, uint16_t>(m);
register_calculate_eta<int, 2, 2, uint16_t>(m);
register_calculate_eta<float, 2, 2, uint16_t>(m);
register_calculate_eta<double, 2, 2, uint16_t>(m);
}

View File

@ -9,8 +9,7 @@
namespace py = pybind11;
template <typename SUM_TYPE>
void define_pedestal_bindings(py::module &m, const std::string &name) {
template <typename SUM_TYPE> void define_pedestal_bindings(py::module &m, const std::string &name) {
py::class_<Pedestal<SUM_TYPE>>(m, name.c_str())
.def(py::init<int, int, int>())
.def(py::init<int, int>())
@ -20,18 +19,16 @@ void define_pedestal_bindings(py::module &m, const std::string &name) {
*mea = self.mean();
return return_image_data(mea);
})
.def("variance",
[](Pedestal<SUM_TYPE> &self) {
auto var = new NDArray<SUM_TYPE, 2>{};
*var = self.variance();
return return_image_data(var);
})
.def("std",
[](Pedestal<SUM_TYPE> &self) {
auto std = new NDArray<SUM_TYPE, 2>{};
*std = self.std();
return return_image_data(std);
})
.def("variance", [](Pedestal<SUM_TYPE> &self) {
auto var = new NDArray<SUM_TYPE, 2>{};
*var = self.variance();
return return_image_data(var);
})
.def("std", [](Pedestal<SUM_TYPE> &self) {
auto std = new NDArray<SUM_TYPE, 2>{};
*std = self.std();
return return_image_data(std);
})
.def("clear", py::overload_cast<>(&Pedestal<SUM_TYPE>::clear))
.def_property_readonly("rows", &Pedestal<SUM_TYPE>::rows)
.def_property_readonly("cols", &Pedestal<SUM_TYPE>::cols)
@ -42,19 +39,14 @@ void define_pedestal_bindings(py::module &m, const std::string &name) {
[&](Pedestal<SUM_TYPE> &pedestal) {
return Pedestal<SUM_TYPE>(pedestal);
})
// TODO! add push for other data types
.def("push",
[](Pedestal<SUM_TYPE> &pedestal, py::array_t<uint16_t> &f) {
auto v = make_view_2d(f);
pedestal.push(v);
})
.def(
"push_no_update",
[](Pedestal<SUM_TYPE> &pedestal,
py::array_t<uint16_t, py::array::c_style> &f) {
auto v = make_view_2d(f);
pedestal.push_no_update(v);
},
py::arg().noconvert())
//TODO! add push for other data types
.def("push", [](Pedestal<SUM_TYPE> &pedestal, py::array_t<uint16_t> &f) {
auto v = make_view_2d(f);
pedestal.push(v);
})
.def("push_no_update", [](Pedestal<SUM_TYPE> &pedestal, py::array_t<uint16_t, py::array::c_style> &f) {
auto v = make_view_2d(f);
pedestal.push_no_update(v);
}, py::arg().noconvert())
.def("update_mean", &Pedestal<SUM_TYPE>::update_mean);
}

View File

@ -1,46 +1,41 @@
#include "aare/PixelMap.hpp"
#include "np_helper.hpp"
#include <cstdint>
#include <pybind11/numpy.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
namespace py = pybind11;
using namespace ::aare;
using namespace::aare;
void define_pixel_map_bindings(py::module &m) {
m.def("GenerateMoench03PixelMap",
[]() {
auto ptr = new NDArray<ssize_t, 2>(GenerateMoench03PixelMap());
return return_image_data(ptr);
})
.def("GenerateMoench05PixelMap",
[]() {
auto ptr = new NDArray<ssize_t, 2>(GenerateMoench05PixelMap());
return return_image_data(ptr);
})
.def("GenerateMoench05PixelMap1g",
[]() {
auto ptr =
new NDArray<ssize_t, 2>(GenerateMoench05PixelMap1g());
return return_image_data(ptr);
})
.def("GenerateMoench05PixelMapOld",
[]() {
auto ptr =
new NDArray<ssize_t, 2>(GenerateMoench05PixelMapOld());
return return_image_data(ptr);
})
.def("GenerateMH02SingleCounterPixelMap",
[]() {
auto ptr = new NDArray<ssize_t, 2>(
GenerateMH02SingleCounterPixelMap());
return return_image_data(ptr);
})
.def("GenerateMH02FourCounterPixelMap", []() {
auto ptr =
new NDArray<ssize_t, 3>(GenerateMH02FourCounterPixelMap());
return return_image_data(ptr);
});
m.def("GenerateMoench03PixelMap", []() {
auto ptr = new NDArray<ssize_t,2>(GenerateMoench03PixelMap());
return return_image_data(ptr);
})
.def("GenerateMoench05PixelMap", []() {
auto ptr = new NDArray<ssize_t,2>(GenerateMoench05PixelMap());
return return_image_data(ptr);
})
.def("GenerateMoench05PixelMap1g", []() {
auto ptr = new NDArray<ssize_t,2>(GenerateMoench05PixelMap1g());
return return_image_data(ptr);
})
.def("GenerateMoench05PixelMapOld", []() {
auto ptr = new NDArray<ssize_t,2>(GenerateMoench05PixelMapOld());
return return_image_data(ptr);
})
.def("GenerateMH02SingleCounterPixelMap", []() {
auto ptr = new NDArray<ssize_t,2>(GenerateMH02SingleCounterPixelMap());
return return_image_data(ptr);
})
.def("GenerateMH02FourCounterPixelMap", []() {
auto ptr = new NDArray<ssize_t,3>(GenerateMH02FourCounterPixelMap());
return return_image_data(ptr);
});
}

View File

@ -32,7 +32,7 @@ void define_raw_file_io_bindings(py::module &m) {
shape.push_back(self.cols());
// return headers from all subfiles
py::array_t<DetectorHeader> header(self.n_modules());
py::array_t<DetectorHeader> header(self.n_mod());
const uint8_t item_size = self.bytes_per_pixel();
if (item_size == 1) {
@ -58,14 +58,13 @@ void define_raw_file_io_bindings(py::module &m) {
throw std::runtime_error("No frames left in file");
}
std::vector<size_t> shape{n_frames, self.rows(), self.cols()};
// return headers from all subfiles
py::array_t<DetectorHeader> header;
if (self.n_modules() == 1) {
if (self.n_mod() == 1) {
header = py::array_t<DetectorHeader>(n_frames);
} else {
header = py::array_t<DetectorHeader>(
{self.n_modules(), n_frames});
header = py::array_t<DetectorHeader>({self.n_mod(), n_frames});
}
// py::array_t<DetectorHeader> header({self.n_mod(), n_frames});
@ -101,7 +100,7 @@ void define_raw_file_io_bindings(py::module &m) {
.def_property_readonly("cols", &RawFile::cols)
.def_property_readonly("bitdepth", &RawFile::bitdepth)
.def_property_readonly("geometry", &RawFile::geometry)
.def_property_readonly("n_modules", &RawFile::n_modules)
.def_property_readonly("n_mod", &RawFile::n_mod)
.def_property_readonly("detector_type", &RawFile::detector_type)
.def_property_readonly("master", &RawFile::master);
}

View File

@ -57,8 +57,7 @@ void define_raw_master_file_bindings(py::module &m) {
.def_property_readonly("total_frames_expected",
&RawMasterFile::total_frames_expected)
.def_property_readonly("geometry", &RawMasterFile::geometry)
.def_property_readonly("analog_samples", &RawMasterFile::analog_samples,
R"(
.def_property_readonly("analog_samples", &RawMasterFile::analog_samples, R"(
Number of analog samples
Returns
@ -67,7 +66,7 @@ void define_raw_master_file_bindings(py::module &m) {
The number of analog samples in the file (or None if not enabled)
)")
.def_property_readonly("digital_samples",
&RawMasterFile::digital_samples, R"(
&RawMasterFile::digital_samples, R"(
Number of digital samples
Returns

View File

@ -24,8 +24,8 @@ auto read_frame_from_RawSubFile(RawSubFile &self) {
py::array_t<DetectorHeader> header(1);
const uint8_t item_size = self.bytes_per_pixel();
std::vector<ssize_t> shape{static_cast<ssize_t>(self.rows()),
static_cast<ssize_t>(self.cols())};
static_cast<ssize_t>(self.cols())};
py::array image;
if (item_size == 1) {
image = py::array_t<uint8_t>(shape);
@ -43,10 +43,12 @@ auto read_frame_from_RawSubFile(RawSubFile &self) {
auto read_n_frames_from_RawSubFile(RawSubFile &self, size_t n_frames) {
py::array_t<DetectorHeader> header(n_frames);
const uint8_t item_size = self.bytes_per_pixel();
std::vector<ssize_t> shape{static_cast<ssize_t>(n_frames),
static_cast<ssize_t>(self.rows()),
static_cast<ssize_t>(self.cols())};
std::vector<ssize_t> shape{
static_cast<ssize_t>(n_frames),
static_cast<ssize_t>(self.rows()),
static_cast<ssize_t>(self.cols())
};
py::array image;
if (item_size == 1) {
image = py::array_t<uint8_t>(shape);
@ -55,14 +57,15 @@ auto read_n_frames_from_RawSubFile(RawSubFile &self, size_t n_frames) {
} else if (item_size == 4) {
image = py::array_t<uint32_t>(shape);
}
self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()),
n_frames, header.mutable_data());
self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()), n_frames,
header.mutable_data());
return py::make_tuple(header, image);
}
// Disable warnings for unused parameters, as we ignore some
// in the __exit__ method
//Disable warnings for unused parameters, as we ignore some
//in the __exit__ method
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
@ -73,7 +76,7 @@ void define_raw_sub_file_io_bindings(py::module &m) {
.def_property_readonly("bytes_per_frame", &RawSubFile::bytes_per_frame)
.def_property_readonly("pixels_per_frame",
&RawSubFile::pixels_per_frame)
.def_property_readonly("bytes_per_pixel", &RawSubFile::bytes_per_pixel)
.def_property_readonly("bytes_per_pixel", &RawSubFile::bytes_per_pixel)
.def("seek", &RawSubFile::seek)
.def("tell", &RawSubFile::tell)
.def_property_readonly("rows", &RawSubFile::rows)
@ -81,17 +84,18 @@ void define_raw_sub_file_io_bindings(py::module &m) {
.def_property_readonly("frames_in_file", &RawSubFile::frames_in_file)
.def("read_frame", &read_frame_from_RawSubFile)
.def("read_n", &read_n_frames_from_RawSubFile)
.def("read",
[](RawSubFile &self) {
self.seek(0);
auto n_frames = self.frames_in_file();
return read_n_frames_from_RawSubFile(self, n_frames);
})
.def("read", [](RawSubFile &self){
self.seek(0);
auto n_frames = self.frames_in_file();
return read_n_frames_from_RawSubFile(self, n_frames);
})
.def("__enter__", [](RawSubFile &self) { return &self; })
.def("__exit__",
[](RawSubFile &self, const std::optional<pybind11::type> &exc_type,
[](RawSubFile &self,
const std::optional<pybind11::type> &exc_type,
const std::optional<pybind11::object> &exc_value,
const std::optional<pybind11::object> &traceback) {})
const std::optional<pybind11::object> &traceback) {
})
.def("__iter__", [](RawSubFile &self) { return &self; })
.def("__next__", [](RawSubFile &self) {
try {
@ -100,6 +104,7 @@ void define_raw_sub_file_io_bindings(py::module &m) {
throw py::stop_iteration();
}
});
}
#pragma GCC diagnostic pop

View File

@ -12,8 +12,10 @@
// #include <pybind11/stl/filesystem.h>
// #include <string>
namespace py = pybind11;
using namespace ::aare;
using namespace::aare;
void define_var_cluster_finder_bindings(py::module &m) {
PYBIND11_NUMPY_DTYPE(VarClusterFinder<double>::Hit, size, row, col,
@ -27,12 +29,12 @@ void define_var_cluster_finder_bindings(py::module &m) {
return return_image_data(ptr);
})
.def("set_noiseMap",
[](VarClusterFinder<double> &self,
[](VarClusterFinder<double> &self,
py::array_t<double, py::array::c_style | py::array::forcecast>
noise_map) {
auto noise_map_span = make_view_2d(noise_map);
self.set_noiseMap(noise_map_span);
})
auto noise_map_span = make_view_2d(noise_map);
self.set_noiseMap(noise_map_span);
})
.def("set_peripheralThresholdFactor",
&VarClusterFinder<double>::set_peripheralThresholdFactor)
.def("find_clusters",
@ -63,7 +65,9 @@ void define_var_cluster_finder_bindings(py::module &m) {
return return_vector(ptr);
})
.def("clear_hits",
[](VarClusterFinder<double> &self) { self.clear_hits(); })
[](VarClusterFinder<double> &self) {
self.clear_hits();
})
.def("steal_hits",
[](VarClusterFinder<double> &self) {
auto ptr = new std::vector<VarClusterFinder<double>::Hit>(
@ -71,4 +75,5 @@ void define_var_cluster_finder_bindings(py::module &m) {
return return_vector(ptr);
})
.def("total_clusters", &VarClusterFinder<double>::total_clusters);
}

View File

@ -5,35 +5,32 @@ from aare import RawSubFile, DetectorType
@pytest.mark.files
def test_read_a_jungfrau_RawSubFile(test_data_path):
# Starting with f1 there is now 7 frames left in the series of files
with RawSubFile(test_data_path / "raw/jungfrau/jungfrau_single_d0_f1_0.raw", DetectorType.Jungfrau, 512, 1024, 16) as f:
assert f.frames_in_file == 7
assert f.frames_in_file == 3
headers, frames = f.read()
assert headers.size == 7
assert frames.shape == (7, 512, 1024)
assert headers.size == 3
assert frames.shape == (3, 512, 1024)
for i,h in zip(range(4,11,1), headers):
# Frame numbers in this file should be 4, 5, 6
for i,h in zip(range(4,7,1), headers):
assert h["frameNumber"] == i
# Compare to canned data using numpy
data = np.load(test_data_path / "raw/jungfrau/jungfrau_single_0.npy")
assert np.all(data[3:] == frames)
assert np.all(data[3:6] == frames)
@pytest.mark.files
def test_iterate_over_a_jungfrau_RawSubFile(test_data_path):
data = np.load(test_data_path / "raw/jungfrau/jungfrau_single_0.npy")
# Given the first subfile in a series we can read all frames from f0, f1, f2...fN
with RawSubFile(test_data_path / "raw/jungfrau/jungfrau_single_d0_f0_0.raw", DetectorType.Jungfrau, 512, 1024, 16) as f:
i = 0
for header, frame in f:
assert header["frameNumber"] == i+1
assert np.all(frame == data[i])
i += 1
assert i == 10
assert header["frameNumber"] == 10
assert i == 3
assert header["frameNumber"] == 3

377
src/AngleCalibration.cpp Normal file
View File

@ -0,0 +1,377 @@
#include "aare/AngleCalibration.hpp"
namespace aare {
AngleCalibration::AngleCalibration(
std::shared_ptr<MythenDetectorSpecifications> mythen_detector_,
std::shared_ptr<FlatField> flat_field_,
std::shared_ptr<MythenFileReader> mythen_file_reader_)
: mythen_detector(mythen_detector_), flat_field(flat_field_),
mythen_file_reader(mythen_file_reader_) {
centers.reserve(mythen_detector->max_modules());
conversions.reserve(mythen_detector->max_modules());
offsets.reserve(mythen_detector->max_modules());
num_bins = std::floor(mythen_detector->max_angle() / histogram_bin_width) -
std::floor(mythen_detector->min_angle() / histogram_bin_width) +
1; // TODO only works if negative
// and positive angle
}
void AngleCalibration::set_histogram_bin_width(double bin_width) {
histogram_bin_width = bin_width;
num_bins = std::floor(mythen_detector->max_angle() / histogram_bin_width) -
std::floor(mythen_detector->min_angle() / histogram_bin_width) +
1; // TODO only works if negative
// and positive angle
}
double AngleCalibration::get_histogram_bin_width() const {
return histogram_bin_width;
}
ssize_t AngleCalibration::get_new_num_bins() const { return num_bins; }
std::vector<double> AngleCalibration::get_centers() const { return centers; }
std::vector<double> AngleCalibration::get_conversions() const {
return conversions;
}
std::vector<double> AngleCalibration::get_offsets() const { return offsets; }
NDView<double, 1> AngleCalibration::get_new_photon_counts() const {
return new_photon_counts.view();
}
NDView<double, 1> AngleCalibration::get_new_statistical_errors() const {
return new_photon_count_errors.view();
}
void AngleCalibration::read_initial_calibration_from_file(
const std::string &filename) {
std::string line;
uint32_t module_number{};
try {
std::ifstream file(filename, std::ios_base::in);
if (!file.good()) {
throw std::logic_error("file does not exist");
}
std::stringstream file_buffer;
file_buffer << file.rdbuf();
while (file_buffer >> line) {
if (line == "module") {
file_buffer >> line;
module_number = std::stoi(line);
}
if (line == "center") {
file_buffer >> line;
centers.insert(centers.begin() + module_number,
std::stod(line));
}
if (line == "conversion") {
file_buffer >> line;
conversions.insert(conversions.begin() + module_number,
std::stod(line));
}
if (line == "offset") {
file_buffer >> line;
offsets.insert(offsets.begin() + module_number,
std::stod(line));
}
}
file.close();
} catch (const std::exception &e) {
std::cerr << "Error: " << e.what() << std::endl;
}
}
parameters AngleCalibration::convert_to_EE_parameters() const {
// normal distance between sample and detector (R)
std::vector<double> normal_distances(centers.size());
// distances between intersection point of sample normal and module origin
// (D)
std::vector<double> module_center_distances(centers.size());
// angles between undiffracted beam and orthogonal sample projection on
// detector (phi)
std::vector<double> angles(centers.size());
for (size_t i = 0; i < centers.size(); ++i) {
auto [module_center_distance, normal_distance, angle] =
convert_to_EE_parameters(i);
normal_distances[i] = normal_distance;
module_center_distances[i] = module_center_distance;
angles[i] = angle;
}
return std::make_tuple(module_center_distances, normal_distances, angles);
}
std::tuple<double, double, double>
AngleCalibration::convert_to_EE_parameters(const size_t module_index) const {
return convert_to_EE_parameters(centers[module_index],
conversions[module_index],
offsets[module_index]);
}
std::tuple<double, double, double> AngleCalibration::convert_to_EE_parameters(
const double center, const double conversion, const double offset) const {
const double module_center_distance =
center * MythenDetectorSpecifications::pitch();
const double normal_distance =
MythenDetectorSpecifications::pitch() / std::abs(conversion);
const double angle = offset + 180.0 / M_PI * center * std::abs(conversion);
return std::make_tuple(module_center_distance, normal_distance, angle);
}
size_t AngleCalibration::global_to_local_strip_index_conversion(
const size_t global_strip_index) const {
const size_t module_index =
global_strip_index / MythenDetectorSpecifications::strips_per_module();
// local strip index in module
size_t local_strip_index =
global_strip_index -
module_index * MythenDetectorSpecifications::strips_per_module();
// switch if indexing is in clock-wise direction
local_strip_index =
std::signbit(conversions[module_index])
? MythenDetectorSpecifications::strips_per_module() - 1 -
local_strip_index
: local_strip_index;
return local_strip_index;
}
/*
parameters
AngleCalibration::convert_to_BC_parameters() {}
*/
double AngleCalibration::diffraction_angle_from_DG_parameters(
const double center, const double conversion, const double offset,
const size_t strip_index, const double distance_to_strip) const {
return offset + 180.0 / M_PI *
(center * std::abs(conversion) -
atan((center - (strip_index + distance_to_strip)) *
std::abs(conversion)));
}
double AngleCalibration::diffraction_angle_from_EE_parameters(
const double module_center_distance, const double normal_distance,
const double angle, const size_t strip_index,
const double distance_to_strip) const {
return angle - 180.0 / M_PI *
atan((module_center_distance -
MythenDetectorSpecifications::pitch() *
(strip_index + distance_to_strip)) /
normal_distance); // TODO: why is it minus
// is it defined counter
// clockwise? thought
// should have a flipped
// sign
}
double AngleCalibration::angular_strip_width_from_DG_parameters(
const double center, const double conversion, const double offset,
const size_t local_strip_index) const {
return std::abs(diffraction_angle_from_DG_parameters(
center, conversion, offset, local_strip_index, -0.5) -
diffraction_angle_from_DG_parameters(
center, conversion, offset, local_strip_index, 0.5));
}
double AngleCalibration::angular_strip_width_from_EE_parameters(
const double module_center_distance, const double normal_distance,
const double angle, const size_t local_strip_index) const {
return std::abs(diffraction_angle_from_EE_parameters(
module_center_distance, normal_distance, angle,
local_strip_index, -0.5) -
diffraction_angle_from_EE_parameters(
module_center_distance, normal_distance, angle,
local_strip_index, 0.5));
// TODO: again not sure about division order - taking abs anyway
}
void AngleCalibration::calculate_fixed_bin_angle_width_histogram(
const size_t start_frame_index, const size_t end_frame_index) {
new_photon_counts = NDArray<double, 1>(std::array<ssize_t, 1>{num_bins});
new_photon_count_errors =
NDArray<double, 1>(std::array<ssize_t, 1>{num_bins});
// TODO: maybe group these into a 2d array - better cache usage
NDArray<double, 1> bin_counts(std::array<ssize_t, 1>{num_bins}, 0.0);
NDArray<double, 1> new_statistical_weights(std::array<ssize_t, 1>{num_bins},
0.0);
NDArray<double, 1> new_errors(std::array<ssize_t, 1>{num_bins}, 0.0);
NDArray<double, 1> inverse_normalized_flatfield =
flat_field->inverse_normalized_flatfield();
for (size_t frame_index = start_frame_index; frame_index < end_frame_index;
++frame_index) {
MythenFrame frame = mythen_file_reader->read_frame(frame_index);
redistribute_photon_counts_to_fixed_angle_bins(
frame, bin_counts.view(), new_statistical_weights.view(),
new_errors.view(), inverse_normalized_flatfield.view());
}
for (ssize_t i = 0; i < new_photon_counts.size(); ++i) {
new_photon_counts[i] = (new_statistical_weights[i] <=
std::numeric_limits<double>::epsilon())
? 0.
: bin_counts[i] / new_statistical_weights[i];
new_photon_count_errors[i] =
(bin_counts[i] <= std::numeric_limits<double>::epsilon())
? 0.
: 1.0 / std::sqrt(bin_counts[i]);
}
}
void AngleCalibration::redistribute_photon_counts_to_fixed_angle_bins(
const MythenFrame &frame, NDView<double, 1> bin_counts,
NDView<double, 1> new_statistical_weights, NDView<double, 1> new_errors,
NDView<double, 1> inverse_normalized_flatfield) const {
ssize_t channel = 0; // TODO handle mask - FlatField still 1d
if (frame.photon_counts.shape()[0] != mythen_detector->num_strips()) {
throw std::runtime_error("wrong number of strips read");
}
ssize_t num_bins1 = mythen_detector->min_angle() / histogram_bin_width;
ssize_t num_bins2 = mythen_detector->max_angle() / histogram_bin_width;
std::cout << "position: " << frame.detector_angle
<< std::endl; // replace with log
double exposure_rate = 1. / mythen_detector->exposure_time();
for (ssize_t strip_index = 0; strip_index < mythen_detector->num_strips();
++strip_index) {
size_t module_index =
strip_index / MythenDetectorSpecifications::strips_per_module();
if (mythen_detector->get_bad_channels()[strip_index] ||
!mythen_detector->get_connected_modules()[module_index])
continue;
double poisson_error = std::sqrt(frame.photon_counts(strip_index)) *
inverse_normalized_flatfield(strip_index) *
exposure_rate;
double corrected_photon_count =
frame.photon_counts(strip_index) *
inverse_normalized_flatfield(strip_index) * exposure_rate;
size_t local_strip_index =
global_to_local_strip_index_conversion(strip_index);
double diffraction_angle = diffraction_angle_from_DG_parameters(
centers[module_index], conversions[module_index],
offsets[module_index], local_strip_index);
diffraction_angle += (frame.detector_angle + mythen_detector->dtt0() +
mythen_detector->bloffset());
if (diffraction_angle < mythen_detector->min_angle() ||
diffraction_angle > mythen_detector->max_angle())
continue;
double angle_covered_by_strip = angular_strip_width_from_DG_parameters(
centers[module_index], conversions[module_index],
offsets[module_index], local_strip_index);
double photon_count_per_bin = histogram_bin_width *
corrected_photon_count /
angle_covered_by_strip;
double error_photon_count_per_bin =
histogram_bin_width * poisson_error / angle_covered_by_strip;
double statistical_weights =
1.0 / std::pow(error_photon_count_per_bin, 2); // 1./sigma²
double strip_boundary_left =
diffraction_angle - 0.5 * angle_covered_by_strip;
double strip_boundary_right =
diffraction_angle + 0.5 * angle_covered_by_strip;
ssize_t left_bin_index = std::max(
num_bins1,
static_cast<ssize_t>(
std::floor(strip_boundary_left / histogram_bin_width) - 1));
ssize_t right_bin_index = std::min(
num_bins2,
static_cast<ssize_t>(
std::ceil(strip_boundary_right / histogram_bin_width) + 1));
// TODO should it be < or <=
for (ssize_t bin = left_bin_index; bin <= right_bin_index; ++bin) {
double bin_coverage = std::min(strip_boundary_right,
(bin + 0.5) * histogram_bin_width) -
std::max(strip_boundary_left,
(bin - 0.5) * histogram_bin_width);
double bin_coverage_factor = bin_coverage / histogram_bin_width;
ssize_t bin_index = bin - num_bins1;
// TODO: maybe have this threshold configurable
if (bin_coverage >= 0.0001) {
new_statistical_weights(bin_index) +=
statistical_weights * bin_coverage_factor;
bin_counts(bin_index) += statistical_weights *
bin_coverage_factor *
photon_count_per_bin;
new_errors(bin_index) += statistical_weights *
bin_coverage_factor *
std::pow(photon_count_per_bin, 2);
}
}
}
}
void AngleCalibration::write_to_file(
const std::string &filename, const bool store_nonzero_bins,
const std::filesystem::path &filepath) const {
std::ofstream output_file(filepath / filename);
if (!output_file) {
std::cerr << "Error opening file!"
<< std::endl; // TODO: replace with log
}
output_file << std::fixed << std::setprecision(15);
for (ssize_t i = 0; i < num_bins; ++i) {
if (new_photon_counts[i] <= std::numeric_limits<double>::epsilon() &&
store_nonzero_bins) {
continue;
}
output_file << std::floor(mythen_detector->min_angle() /
histogram_bin_width) *
histogram_bin_width +
i * histogram_bin_width
<< " " << new_photon_counts[i] << " "
<< new_photon_count_errors[i] << std::endl;
}
output_file.close();
}
} // namespace aare

View File

@ -0,0 +1,234 @@
/************************************************
* @file AngleCalibration.test.cpp
* @short test case for angle calibration class
***********************************************/
#include "aare/AngleCalibration.hpp"
#include <filesystem>
#include "test_config.hpp"
#include <iomanip>
#include <type_traits>
#include <catch2/catch_all.hpp>
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
using namespace aare;
TEST_CASE("read initial angle calibration file",
"[.anglecalibration] [.files]") {
std::shared_ptr<MythenDetectorSpecifications> mythen_detector_ptr =
std::make_shared<MythenDetectorSpecifications>();
AngleCalibration anglecalibration(mythen_detector_ptr,
std::shared_ptr<FlatField>{},
std::shared_ptr<MythenFileReader>{});
std::string filename = test_data_path() / "AngleCalibration_Test_Data" /
"Angcal_2E_Feb2023_P29.off";
REQUIRE(std::filesystem::exists(filename));
anglecalibration.read_initial_calibration_from_file(filename);
auto centers = anglecalibration.get_centers();
auto conversions = anglecalibration.get_conversions();
auto offsets = anglecalibration.get_offsets();
std::cout.precision(17);
CHECK(centers.size() == 48);
CHECK(conversions.size() == 48);
CHECK(offsets.size() == 48);
CHECK(centers[9] == Catch::Approx(660.342326));
CHECK(offsets[47] == Catch::Approx(5.8053312));
CHECK(conversions[27] == Catch::Approx(-0.6581179125e-4));
}
TEST_CASE("read bad channels",
"[.anglecalibration][.mythenspecifications][.files]") {
MythenDetectorSpecifications mythen_detector;
std::string bad_channels_filename = test_data_path() /
"AngleCalibration_Test_Data" /
"bc2023_003_RING.chans";
REQUIRE(std::filesystem::exists(bad_channels_filename));
mythen_detector.read_bad_channels_from_file(bad_channels_filename);
CHECK(mythen_detector.get_bad_channels().size() == 61440);
CHECK(mythen_detector.get_bad_channels()[61437] == true);
CHECK(std::all_of(mythen_detector.get_bad_channels().begin() + 30720,
mythen_detector.get_bad_channels().begin() + 61439,
[](const bool element) { return element; }));
}
TEST_CASE("read unconnected modules",
"[.anglecalibration][.mythenspecifications][.files]") {
MythenDetectorSpecifications mythen_detector;
std::string unconnected_modules_filename =
test_data_path() / "AngleCalibration_Test_Data" / "ModOut.txt";
REQUIRE(std::filesystem::exists(unconnected_modules_filename));
mythen_detector.read_unconnected_modules_from_file(
unconnected_modules_filename);
CHECK(mythen_detector.get_connected_modules().size() == 48);
CHECK(std::all_of(mythen_detector.get_connected_modules().begin(),
mythen_detector.get_connected_modules().end(),
[](const bool element) { return element; }));
}
TEST_CASE("read flatfield", "[.anglecalibration][.flatfield][.files]") {
std::shared_ptr<MythenDetectorSpecifications> mythen_detector_ptr =
std::make_shared<MythenDetectorSpecifications>();
FlatField flatfield(mythen_detector_ptr);
std::string flatfield_filename =
test_data_path() / "AngleCalibration_Test_Data" /
"Flatfield_E22p0keV_T11000eV_up_48M_a_LONG_Feb2023_open_WS_SUMC.raw";
REQUIRE(std::filesystem::exists(flatfield_filename));
flatfield.read_flatfield_from_file(flatfield_filename);
auto flatfield_data = flatfield.get_flatfield();
CHECK(flatfield_data.size() == 61440);
CHECK(flatfield_data[0] == 0);
CHECK(flatfield_data[21] == 4234186);
}
TEST_CASE("compare result with python code", "[.anglecalibration] [.files]") {
auto fpath = test_data_path() / "AngleCalibration_Test_Data";
REQUIRE(std::filesystem::exists(fpath));
std::shared_ptr<MythenDetectorSpecifications> mythen_detector_ptr =
std::make_shared<MythenDetectorSpecifications>();
std::string bad_channels_filename = fpath / "bc2023_003_RING.chans";
REQUIRE(std::filesystem::exists(bad_channels_filename));
mythen_detector_ptr->read_bad_channels_from_file(bad_channels_filename);
std::string unconnected_modules_filename = fpath / "ModOut.txt";
REQUIRE(std::filesystem::exists(unconnected_modules_filename));
mythen_detector_ptr->read_unconnected_modules_from_file(
unconnected_modules_filename);
std::shared_ptr<FlatField> flat_field_ptr =
std::make_shared<FlatField>(mythen_detector_ptr);
std::string flatfield_filename =
fpath /
"Flatfield_E22p0keV_T11000eV_up_48M_a_LONG_Feb2023_open_WS_SUMC.raw";
REQUIRE(std::filesystem::exists(flatfield_filename));
flat_field_ptr->read_flatfield_from_file(flatfield_filename);
std::shared_ptr<MythenFileReader> mythen_file_reader_ptr =
std::make_shared<MythenFileReader>(fpath,
"ang1up_22keV_LaB60p3mm_48M_a_0");
AngleCalibration anglecalibration(mythen_detector_ptr, flat_field_ptr,
mythen_file_reader_ptr);
std::string initial_angles_filename = fpath / "Angcal_2E_Feb2023_P29.off";
REQUIRE(std::filesystem::exists(initial_angles_filename));
anglecalibration.read_initial_calibration_from_file(
initial_angles_filename);
anglecalibration.calculate_fixed_bin_angle_width_histogram(320, 340);
// anglecalibration.write_to_file("cpp_new_photon_counts.xye");
auto expected_filename_photons =
test_data_path() / "AngleCalibration_Test_Data" / "new_photons.bin";
REQUIRE(std::filesystem::exists(expected_filename_photons));
auto expected_filename_errors =
test_data_path() / "AngleCalibration_Test_Data" / "new_errors.bin";
REQUIRE(std::filesystem::exists(expected_filename_errors));
ssize_t new_num_bins = anglecalibration.get_new_num_bins();
auto python_output_errors = load<double, 1>(
expected_filename_errors, std::array<ssize_t, 1>{new_num_bins});
auto python_output_photons = load<double, 1>(
expected_filename_photons, std::array<ssize_t, 1>{new_num_bins});
CHECK(anglecalibration.get_new_photon_counts().equals(
python_output_photons.view(),
1e-8)); // not sure about precision does not exactly match to all
// decimal digits
CHECK(anglecalibration.get_new_statistical_errors().equals(
python_output_errors.view(),
1e-8)); //
}
TEST_CASE("check conversion from DG to EE parameters", "[.anglecalibration]") {
std::shared_ptr<MythenDetectorSpecifications> mythen_detector_ptr =
std::make_shared<MythenDetectorSpecifications>();
AngleCalibration anglecalibration(mythen_detector_ptr,
std::shared_ptr<FlatField>{},
std::shared_ptr<MythenFileReader>{});
// DG test parameters
const double center = 642.197591224993;
const double conversion = 0.657694036246975e-4;
const double offset = 5.004892881251670;
const ssize_t local_strip_index = 1;
double diffraction_angle_DG_param =
anglecalibration.diffraction_angle_from_DG_parameters(
center, conversion, offset, local_strip_index);
auto [distance_center, normal_distance, angle] =
anglecalibration.convert_to_EE_parameters(center, conversion, offset);
double diffraction_angle_EE_param =
anglecalibration.diffraction_angle_from_EE_parameters(
distance_center, normal_distance, angle, local_strip_index);
CHECK(diffraction_angle_EE_param ==
Catch::Approx(diffraction_angle_DG_param));
double strip_width_DG_param =
anglecalibration.angular_strip_width_from_DG_parameters(
center, conversion, offset, local_strip_index);
double strip_width_EE_param =
anglecalibration.angular_strip_width_from_EE_parameters(
distance_center, normal_distance, angle, local_strip_index);
CHECK(strip_width_DG_param == Catch::Approx(strip_width_EE_param));
}

View File

@ -1,395 +0,0 @@
#include "aare/ClusterFile.hpp"
#include <algorithm>
namespace aare {
ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size,
const std::string &mode)
: m_chunk_size(chunk_size), m_mode(mode) {
if (mode == "r") {
fp = fopen(fname.c_str(), "rb");
if (!fp) {
throw std::runtime_error("Could not open file for reading: " +
fname.string());
}
} else if (mode == "w") {
fp = fopen(fname.c_str(), "wb");
if (!fp) {
throw std::runtime_error("Could not open file for writing: " +
fname.string());
}
} else if (mode == "a") {
fp = fopen(fname.c_str(), "ab");
if (!fp) {
throw std::runtime_error("Could not open file for appending: " +
fname.string());
}
} else {
throw std::runtime_error("Unsupported mode: " + mode);
}
}
void ClusterFile::set_roi(ROI roi) { m_roi = roi; }
void ClusterFile::set_noise_map(const NDView<int32_t, 2> noise_map) {
m_noise_map = NDArray<int32_t, 2>(noise_map);
}
void ClusterFile::set_gain_map(const NDView<double, 2> gain_map) {
m_gain_map = NDArray<double, 2>(gain_map);
// Gain map is passed as ADU/keV to avoid dividing in when applying the gain
// map we invert it here
for (auto &item : m_gain_map->view()) {
item = 1.0 / item;
}
}
ClusterFile::~ClusterFile() { close(); }
void ClusterFile::close() {
if (fp) {
fclose(fp);
fp = nullptr;
}
}
void ClusterFile::write_frame(const ClusterVector<int32_t> &clusters) {
if (m_mode != "w" && m_mode != "a") {
throw std::runtime_error("File not opened for writing");
}
if (!(clusters.cluster_size_x() == 3) &&
!(clusters.cluster_size_y() == 3)) {
throw std::runtime_error("Only 3x3 clusters are supported");
}
// First write the frame number - 4 bytes
int32_t frame_number = clusters.frame_number();
if (fwrite(&frame_number, sizeof(frame_number), 1, fp) != 1) {
throw std::runtime_error(LOCATION + "Could not write frame number");
}
// Then write the number of clusters - 4 bytes
uint32_t n_clusters = clusters.size();
if (fwrite(&n_clusters, sizeof(n_clusters), 1, fp) != 1) {
throw std::runtime_error(LOCATION +
"Could not write number of clusters");
}
// Now write the clusters in the frame
if (fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp) !=
clusters.size()) {
throw std::runtime_error(LOCATION + "Could not write clusters");
}
}
ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters) {
if (m_mode != "r") {
throw std::runtime_error("File not opened for reading");
}
if (m_noise_map || m_roi) {
return read_clusters_with_cut(n_clusters);
} else {
return read_clusters_without_cut(n_clusters);
}
}
ClusterVector<int32_t>
ClusterFile::read_clusters_without_cut(size_t n_clusters) {
if (m_mode != "r") {
throw std::runtime_error("File not opened for reading");
}
ClusterVector<int32_t> clusters(3, 3, n_clusters);
int32_t iframe = 0; // frame number needs to be 4 bytes!
size_t nph_read = 0;
uint32_t nn = m_num_left;
uint32_t nph = m_num_left; // number of clusters in frame needs to be 4
// auto buf = reinterpret_cast<Cluster3x3 *>(clusters.data());
auto buf = clusters.data();
// if there are photons left from previous frame read them first
if (nph) {
if (nph > n_clusters) {
// if we have more photons left in the frame then photons to read we
// read directly the requested number
nn = n_clusters;
} else {
nn = nph;
}
nph_read += fread((buf + nph_read * clusters.item_size()),
clusters.item_size(), nn, fp);
m_num_left = nph - nn; // write back the number of photons left
}
if (nph_read < n_clusters) {
// keep on reading frames and photons until reaching n_clusters
while (fread(&iframe, sizeof(iframe), 1, fp)) {
clusters.set_frame_number(iframe);
// read number of clusters in frame
if (fread(&nph, sizeof(nph), 1, fp)) {
if (nph > (n_clusters - nph_read))
nn = n_clusters - nph_read;
else
nn = nph;
nph_read += fread((buf + nph_read * clusters.item_size()),
clusters.item_size(), nn, fp);
m_num_left = nph - nn;
}
if (nph_read >= n_clusters)
break;
}
}
// Resize the vector to the number of clusters.
// No new allocation, only change bounds.
clusters.resize(nph_read);
if (m_gain_map)
clusters.apply_gain_map(m_gain_map->view());
return clusters;
}
ClusterVector<int32_t> ClusterFile::read_clusters_with_cut(size_t n_clusters) {
ClusterVector<int32_t> clusters(3, 3);
clusters.reserve(n_clusters);
// if there are photons left from previous frame read them first
if (m_num_left) {
while (m_num_left && clusters.size() < n_clusters) {
Cluster3x3 c = read_one_cluster();
if (is_selected(c)) {
clusters.push_back(c.x, c.y,
reinterpret_cast<std::byte *>(c.data));
}
}
}
// we did not have enough clusters left in the previous frame
// keep on reading frames until reaching n_clusters
if (clusters.size() < n_clusters) {
// sanity check
if (m_num_left) {
throw std::runtime_error(
LOCATION + "Entered second loop with clusters left\n");
}
int32_t frame_number = 0; // frame number needs to be 4 bytes!
while (fread(&frame_number, sizeof(frame_number), 1, fp)) {
if (fread(&m_num_left, sizeof(m_num_left), 1, fp)) {
clusters.set_frame_number(
frame_number); // cluster vector will hold the last frame
// number
while (m_num_left && clusters.size() < n_clusters) {
Cluster3x3 c = read_one_cluster();
if (is_selected(c)) {
clusters.push_back(
c.x, c.y, reinterpret_cast<std::byte *>(c.data));
}
}
}
// we have enough clusters, break out of the outer while loop
if (clusters.size() >= n_clusters)
break;
}
}
if (m_gain_map)
clusters.apply_gain_map(m_gain_map->view());
return clusters;
}
Cluster3x3 ClusterFile::read_one_cluster() {
Cluster3x3 c;
auto rc = fread(&c, sizeof(c), 1, fp);
if (rc != 1) {
throw std::runtime_error(LOCATION + "Could not read cluster");
}
--m_num_left;
return c;
}
ClusterVector<int32_t> ClusterFile::read_frame() {
if (m_mode != "r") {
throw std::runtime_error(LOCATION + "File not opened for reading");
}
if (m_noise_map || m_roi) {
return read_frame_with_cut();
} else {
return read_frame_without_cut();
}
}
ClusterVector<int32_t> ClusterFile::read_frame_without_cut() {
if (m_mode != "r") {
throw std::runtime_error("File not opened for reading");
}
if (m_num_left) {
throw std::runtime_error(
"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");
}
int32_t n_clusters; // Saved as 32bit integer in the cluster file
if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) {
throw std::runtime_error(LOCATION +
"Could not read number of clusters");
}
ClusterVector<int32_t> clusters(3, 3, n_clusters);
clusters.set_frame_number(frame_number);
if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) !=
static_cast<size_t>(n_clusters)) {
throw std::runtime_error(LOCATION + "Could not read clusters");
}
clusters.resize(n_clusters);
if (m_gain_map)
clusters.apply_gain_map(m_gain_map->view());
return clusters;
}
ClusterVector<int32_t> ClusterFile::read_frame_with_cut() {
if (m_mode != "r") {
throw std::runtime_error("File not opened for reading");
}
if (m_num_left) {
throw std::runtime_error(
"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("Could not read frame number");
}
if (fread(&m_num_left, sizeof(m_num_left), 1, fp) != 1) {
throw std::runtime_error("Could not read number of clusters");
}
ClusterVector<int32_t> clusters(3, 3);
clusters.reserve(m_num_left);
clusters.set_frame_number(frame_number);
while (m_num_left) {
Cluster3x3 c = read_one_cluster();
if (is_selected(c)) {
clusters.push_back(c.x, c.y, reinterpret_cast<std::byte *>(c.data));
}
}
if (m_gain_map)
clusters.apply_gain_map(m_gain_map->view());
return clusters;
}
bool ClusterFile::is_selected(Cluster3x3 &cl) {
// Should fail fast
if (m_roi) {
if (!(m_roi->contains(cl.x, cl.y))) {
return false;
}
}
if (m_noise_map) {
int32_t sum_1x1 = cl.data[4]; // central pixel
int32_t sum_2x2 = cl.sum_2x2(); // highest sum of 2x2 subclusters
int32_t sum_3x3 = cl.sum(); // sum of all pixels
auto noise =
(*m_noise_map)(cl.y, cl.x); // TODO! check if this is correct
if (sum_1x1 <= noise || sum_2x2 <= 2 * noise || sum_3x3 <= 3 * noise) {
return false;
}
}
// we passed all checks
return true;
}
NDArray<double, 2> calculate_eta2(ClusterVector<int> &clusters) {
// TOTO! make work with 2x2 clusters
NDArray<double, 2> eta2({static_cast<int64_t>(clusters.size()), 2});
if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) {
for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_eta2(clusters.at<Cluster3x3>(i));
eta2(i, 0) = e.x;
eta2(i, 1) = e.y;
}
} else if (clusters.cluster_size_x() == 2 ||
clusters.cluster_size_y() == 2) {
for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_eta2(clusters.at<Cluster2x2>(i));
eta2(i, 0) = e.x;
eta2(i, 1) = e.y;
}
} else {
throw std::runtime_error("Only 3x3 and 2x2 clusters are supported");
}
return eta2;
}
/**
* @brief Calculate the eta2 values for a 3x3 cluster and return them in a Eta2
* struct containing etay, etax and the corner of the cluster.
*/
Eta2 calculate_eta2(Cluster3x3 &cl) {
Eta2 eta{};
std::array<int32_t, 4> tot2;
tot2[0] = cl.data[0] + cl.data[1] + cl.data[3] + cl.data[4];
tot2[1] = cl.data[1] + cl.data[2] + cl.data[4] + cl.data[5];
tot2[2] = cl.data[3] + cl.data[4] + cl.data[6] + cl.data[7];
tot2[3] = cl.data[4] + cl.data[5] + cl.data[7] + cl.data[8];
auto c = std::max_element(tot2.begin(), tot2.end()) - tot2.begin();
eta.sum = tot2[c];
switch (c) {
case cBottomLeft:
if ((cl.data[3] + cl.data[4]) != 0)
eta.x = static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
if ((cl.data[1] + cl.data[4]) != 0)
eta.y = static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
eta.c = cBottomLeft;
break;
case cBottomRight:
if ((cl.data[2] + cl.data[5]) != 0)
eta.x = static_cast<double>(cl.data[5]) / (cl.data[4] + cl.data[5]);
if ((cl.data[1] + cl.data[4]) != 0)
eta.y = static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
eta.c = cBottomRight;
break;
case cTopLeft:
if ((cl.data[7] + cl.data[4]) != 0)
eta.x = static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
if ((cl.data[7] + cl.data[4]) != 0)
eta.y = static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
eta.c = cTopLeft;
break;
case cTopRight:
if ((cl.data[5] + cl.data[4]) != 0)
eta.x = static_cast<double>(cl.data[5]) / (cl.data[5] + cl.data[4]);
if ((cl.data[7] + cl.data[4]) != 0)
eta.y = static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
eta.c = cTopRight;
break;
// no default to allow compiler to warn about missing cases
}
return eta;
}
Eta2 calculate_eta2(Cluster2x2 &cl) {
Eta2 eta{};
if ((cl.data[0] + cl.data[1]) != 0)
eta.x = static_cast<double>(cl.data[1]) / (cl.data[0] + cl.data[1]);
if ((cl.data[0] + cl.data[2]) != 0)
eta.y = static_cast<double>(cl.data[2]) / (cl.data[0] + cl.data[2]);
eta.sum = cl.data[0] + cl.data[1] + cl.data[2] + cl.data[3];
eta.c = cBottomLeft; // TODO! This is not correct, but need to put something
return eta;
}
} // namespace aare

View File

@ -10,8 +10,9 @@ using aare::Cluster;
using aare::ClusterFile;
using aare::ClusterVector;
TEST_CASE("Read one frame from a cluster file", "[.files]") {
// We know that the frame has 97 clusters
//We know that the frame has 97 clusters
auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust";
REQUIRE(std::filesystem::exists(fpath));
@ -26,6 +27,7 @@ TEST_CASE("Read one frame from a cluster file", "[.files]") {
std::begin(expected_cluster_data)));
}
TEST_CASE("Read one frame using ROI", "[.files]") {
// We know that the frame has 97 clusters
auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust";
@ -58,6 +60,8 @@ TEST_CASE("Read one frame using ROI", "[.files]") {
std::begin(expected_cluster_data)));
}
TEST_CASE("Read clusters from single frame file", "[.files]") {
// frame_number, num_clusters [135] 97

View File

@ -14,24 +14,22 @@ CtbRawFile::CtbRawFile(const std::filesystem::path &fname) : m_master(fname) {
m_file.open(m_master.data_fname(0, 0), std::ios::binary);
}
void CtbRawFile::read_into(std::byte *image_buf, DetectorHeader *header) {
if (m_current_frame >= m_master.frames_in_file()) {
void CtbRawFile::read_into(std::byte *image_buf, DetectorHeader* header) {
if(m_current_frame >= m_master.frames_in_file()){
throw std::runtime_error(LOCATION + " End of file reached");
}
if (m_current_frame != 0 &&
m_current_frame % m_master.max_frames_per_file() == 0) {
open_data_file(m_current_subfile + 1);
if(m_current_frame != 0 && m_current_frame % m_master.max_frames_per_file() == 0){
open_data_file(m_current_subfile+1);
}
if (header) {
if(header){
m_file.read(reinterpret_cast<char *>(header), sizeof(DetectorHeader));
} else {
}else{
m_file.seekg(sizeof(DetectorHeader), std::ios::cur);
}
m_file.read(reinterpret_cast<char *>(image_buf),
m_master.image_size_in_bytes());
m_file.read(reinterpret_cast<char *>(image_buf), m_master.image_size_in_bytes());
m_current_frame++;
}
@ -40,16 +38,13 @@ void CtbRawFile::seek(size_t frame_number) {
open_data_file(index);
}
size_t frame_number_in_file = frame_number % m_master.max_frames_per_file();
m_file.seekg((sizeof(DetectorHeader) + m_master.image_size_in_bytes()) *
frame_number_in_file);
m_file.seekg((sizeof(DetectorHeader)+m_master.image_size_in_bytes()) * frame_number_in_file);
m_current_frame = frame_number;
}
size_t CtbRawFile::tell() const { return m_current_frame; }
size_t CtbRawFile::image_size_in_bytes() const {
return m_master.image_size_in_bytes();
}
size_t CtbRawFile::image_size_in_bytes() const { return m_master.image_size_in_bytes(); }
size_t CtbRawFile::frames_in_file() const { return m_master.frames_in_file(); }
@ -68,11 +63,12 @@ void CtbRawFile::open_data_file(size_t subfile_index) {
throw std::runtime_error(LOCATION + "Subfile index out of range");
}
m_current_subfile = subfile_index;
m_file = std::ifstream(m_master.data_fname(0, subfile_index),
std::ios::binary); // only one module for CTB
m_file = std::ifstream(m_master.data_fname(0, subfile_index), std::ios::binary); // only one module for CTB
if (!m_file.is_open()) {
throw std::runtime_error(LOCATION + "Could not open data file");
}
}
} // namespace aare

View File

@ -10,8 +10,7 @@ namespace aare {
* @brief Construct a DType object from a type_info object
* @param t type_info object
* @throw runtime_error if the type is not supported
* @note supported types are: int8_t, uint8_t, int16_t, uint16_t, int32_t,
* uint32_t, int64_t, uint64_t, float, double
* @note supported types are: int8_t, uint8_t, int16_t, uint16_t, int32_t, uint32_t, int64_t, uint64_t, float, double
* @note the type_info object is obtained using typeid (e.g. typeid(int))
*/
Dtype::Dtype(const std::type_info &t) {
@ -36,8 +35,7 @@ Dtype::Dtype(const std::type_info &t) {
else if (t == typeid(double))
m_type = TypeIndex::DOUBLE;
else
throw std::runtime_error(
"Could not construct data type. Type not supported.");
throw std::runtime_error("Could not construct data type. Type not supported.");
}
/**
@ -65,8 +63,7 @@ uint8_t Dtype::bitdepth() const {
case TypeIndex::NONE:
return 0;
default:
throw std::runtime_error(LOCATION +
"Could not get bitdepth. Type not supported.");
throw std::runtime_error(LOCATION + "Could not get bitdepth. Type not supported.");
}
}
@ -141,8 +138,7 @@ Dtype Dtype::from_bitdepth(uint8_t bitdepth) {
case 64:
return Dtype(TypeIndex::UINT64);
default:
throw std::runtime_error(
"Could not construct data type from bitdepth.");
throw std::runtime_error("Could not construct data type from bitdepth.");
}
}
/**
@ -179,27 +175,17 @@ std::string Dtype::to_string() const {
case TypeIndex::DOUBLE:
return "f8";
case TypeIndex::ERROR:
throw std::runtime_error(
"Could not get string representation. Type not supported.");
throw std::runtime_error("Could not get string representation. Type not supported.");
case TypeIndex::NONE:
throw std::runtime_error(
"Could not get string representation. Type not supported.");
throw std::runtime_error("Could not get string representation. Type not supported.");
}
return {};
}
bool Dtype::operator==(const Dtype &other) const noexcept {
return m_type == other.m_type;
}
bool Dtype::operator!=(const Dtype &other) const noexcept {
return !(*this == other);
}
bool Dtype::operator==(const Dtype &other) const noexcept { return m_type == other.m_type; }
bool Dtype::operator!=(const Dtype &other) const noexcept { return !(*this == other); }
bool Dtype::operator==(const std::type_info &t) const {
return Dtype(t) == *this;
}
bool Dtype::operator!=(const std::type_info &t) const {
return Dtype(t) != *this;
}
bool Dtype::operator==(const std::type_info &t) const { return Dtype(t) == *this; }
bool Dtype::operator!=(const std::type_info &t) const { return Dtype(t) != *this; }
} // namespace aare

View File

@ -51,6 +51,4 @@ TEST_CASE("Construct from string with endianess") {
REQUIRE_THROWS(Dtype(">i4") == typeid(int32_t));
}
TEST_CASE("Convert to string") {
REQUIRE(Dtype(typeid(int)).to_string() == "<i4");
}
TEST_CASE("Convert to string") { REQUIRE(Dtype(typeid(int)).to_string() == "<i4"); }

View File

@ -19,24 +19,28 @@ File::File(const std::filesystem::path &fname, const std::string &mode,
fmt::format("File does not exist: {}", fname.string()));
}
// Assuming we are pointing at a master file?
// Assuming we are pointing at a master file?
// TODO! How do we read raw files directly?
if (fname.extension() == ".raw" || fname.extension() == ".json") {
// file_impl = new RawFile(fname, mode, cfg);
file_impl = std::make_unique<RawFile>(fname, mode);
} else if (fname.extension() == ".npy") {
}
else if (fname.extension() == ".npy") {
// file_impl = new NumpyFile(fname, mode, cfg);
file_impl = std::make_unique<NumpyFile>(fname, mode, cfg);
} else if (fname.extension() == ".dat") {
}else if(fname.extension() == ".dat"){
file_impl = std::make_unique<JungfrauDataFile>(fname);
} else {
throw std::runtime_error("Unsupported file type");
}
}
File::File(File &&other) noexcept { std::swap(file_impl, other.file_impl); }
File &File::operator=(File &&other) noexcept {
File::File(File &&other) noexcept{
std::swap(file_impl, other.file_impl);
}
File& File::operator=(File &&other) noexcept {
if (this != &other) {
File tmp(std::move(other));
std::swap(file_impl, tmp.file_impl);
@ -66,16 +70,15 @@ size_t File::frame_number(size_t frame_index) {
}
size_t File::bytes_per_frame() const { return file_impl->bytes_per_frame(); }
size_t File::pixels_per_frame() const { return file_impl->pixels_per_frame(); }
size_t File::pixels_per_frame() const{ return file_impl->pixels_per_frame(); }
void File::seek(size_t frame_index) { file_impl->seek(frame_index); }
size_t File::tell() const { return file_impl->tell(); }
size_t File::rows() const { return file_impl->rows(); }
size_t File::cols() const { return file_impl->cols(); }
size_t File::bitdepth() const { return file_impl->bitdepth(); }
size_t File::bytes_per_pixel() const {
return file_impl->bitdepth() / bits_per_byte;
}
size_t File::bytes_per_pixel() const { return file_impl->bitdepth() / bits_per_byte; }
DetectorType File::detector_type() const { return file_impl->detector_type(); }
} // namespace aare

View File

@ -6,12 +6,10 @@
namespace aare {
FilePtr::FilePtr(const std::filesystem::path &fname,
const std::string &mode = "rb") {
FilePtr::FilePtr(const std::filesystem::path& fname, const std::string& mode = "rb") {
fp_ = fopen(fname.c_str(), mode.c_str());
if (!fp_)
throw std::runtime_error(
fmt::format("Could not open: {}", fname.c_str()));
throw std::runtime_error(fmt::format("Could not open: {}", fname.c_str()));
}
FilePtr::FilePtr(FilePtr &&other) { std::swap(fp_, other.fp_); }
@ -26,16 +24,15 @@ FILE *FilePtr::get() { return fp_; }
ssize_t FilePtr::tell() {
auto pos = ftell(fp_);
if (pos == -1)
throw std::runtime_error(
fmt::format("Error getting file position: {}", error_msg()));
throw std::runtime_error(fmt::format("Error getting file position: {}", error_msg()));
return pos;
}
}
FilePtr::~FilePtr() {
if (fp_)
fclose(fp_); // check?
}
std::string FilePtr::error_msg() {
std::string FilePtr::error_msg(){
if (feof(fp_)) {
return "End of file reached";
}

View File

@ -1,12 +1,13 @@
#include "aare/Fit.hpp"
#include "aare/utils/par.hpp"
#include "aare/utils/task.hpp"
#include "aare/utils/par.hpp"
#include <lmcurve2.h>
#include <lmfit.hpp>
#include <thread>
#include <array>
namespace aare {
namespace func {
@ -33,10 +34,8 @@ NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par) {
return y;
}
double scurve(const double x, const double *par) {
return (par[0] + par[1] * x) +
0.5 * (1 + erf((x - par[2]) / (sqrt(2) * par[3]))) *
(par[4] + par[5] * (x - par[2]));
double scurve(const double x, const double * par) {
return (par[0] + par[1] * x) + 0.5 * (1 + erf((x - par[2]) / (sqrt(2) * par[3]))) * (par[4] + par[5] * (x - par[2]));
}
NDArray<double, 1> scurve(NDView<double, 1> x, NDView<double, 1> par) {
@ -47,10 +46,8 @@ NDArray<double, 1> scurve(NDView<double, 1> x, NDView<double, 1> par) {
return y;
}
double scurve2(const double x, const double *par) {
return (par[0] + par[1] * x) +
0.5 * (1 - erf((x - par[2]) / (sqrt(2) * par[3]))) *
(par[4] + par[5] * (x - par[2]));
double scurve2(const double x, const double * par) {
return (par[0] + par[1] * x) + 0.5 * (1 - erf((x - par[2]) / (sqrt(2) * par[3]))) * (par[4] + par[5] * (x - par[2]));
}
NDArray<double, 1> scurve2(NDView<double, 1> x, NDView<double, 1> par) {
@ -94,8 +91,7 @@ NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
return result;
}
std::array<double, 3> gaus_init_par(const NDView<double, 1> x,
const NDView<double, 1> y) {
std::array<double, 3> gaus_init_par(const NDView<double, 1> x, const NDView<double, 1> y) {
std::array<double, 3> start_par{0, 0, 0};
auto e = std::max_element(y.begin(), y.end());
auto idx = std::distance(y.begin(), e);
@ -107,29 +103,31 @@ std::array<double, 3> gaus_init_par(const NDView<double, 1> x,
// For sigma we estimate the fwhm and divide by 2.35
// assuming equally spaced x values
auto delta = x[1] - x[0];
start_par[2] = std::count_if(y.begin(), y.end(),
[e](double val) { return val > *e / 2; }) *
delta / 2.35;
start_par[2] =
std::count_if(y.begin(), y.end(),
[e, delta](double val) { return val > *e / 2; }) *
delta / 2.35;
return start_par;
}
std::array<double, 2> pol1_init_par(const NDView<double, 1> x,
const NDView<double, 1> y) {
// Estimate the initial parameters for the fit
std::array<double, 2> start_par{0, 0};
auto y2 = std::max_element(y.begin(), y.end());
auto x2 = x[std::distance(y.begin(), y2)];
auto y1 = std::min_element(y.begin(), y.end());
auto x1 = x[std::distance(y.begin(), y1)];
std::array<double, 2> pol1_init_par(const NDView<double, 1> x, const NDView<double, 1> y){
// Estimate the initial parameters for the fit
std::array<double, 2> start_par{0, 0};
start_par[0] =
(*y2 - *y1) / (x2 - x1); // For amplitude we use the maximum value
start_par[1] =
*y1 - ((*y2 - *y1) / (x2 - x1)) *
x1; // For the mean we use the x value of the maximum value
return start_par;
auto y2 = std::max_element(y.begin(), y.end());
auto x2 = x[std::distance(y.begin(), y2)];
auto y1 = std::min_element(y.begin(), y.end());
auto x1 = x[std::distance(y.begin(), y1)];
start_par[0] =
(*y2 - *y1) / (x2 - x1); // For amplitude we use the maximum value
start_par[1] =
*y1 - ((*y2 - *y1) / (x2 - x1)) *
x1; // For the mean we use the x value of the maximum value
return start_par;
}
void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
@ -143,6 +141,7 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
"and par_out, par_err_out must have size 3");
}
// /* Collection of output parameters for status info. */
// typedef struct {
// double fnorm; /* norm of the residue vector fvec. */
@ -154,32 +153,23 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
// */
// } lm_status_struct;
lm_status_struct status;
par_out = gaus_init_par(x, y);
std::array<double, 9> cov{0, 0, 0, 0, 0, 0, 0, 0, 0};
std::array<double, 9> cov{0, 0, 0, 0, 0, 0, 0 , 0 , 0};
// void lmcurve2( const int n_par, double *par, double *parerr, double
// *covar, const int m_dat, const double *t, const double *y, const double
// *dy, double (*f)( const double ti, const double *par ), const
// lm_control_struct *control, lm_status_struct *status); n_par - Number of
// free variables. Length of parameter vector par. par - Parameter vector.
// On input, it must contain a reasonable guess. On output, it contains the
// solution found to minimize ||r||. parerr - Parameter uncertainties
// vector. Array of length n_par or NULL. On output, unless it or covar is
// NULL, it contains the weighted parameter uncertainties for the found
// parameters. covar - Covariance matrix. Array of length n_par * n_par or
// NULL. On output, unless it is NULL, it contains the covariance matrix.
// m_dat - Number of data points. Length of vectors t, y, dy. Must statisfy
// n_par <= m_dat. t - Array of length m_dat. Contains the abcissae (time,
// or "x") for which function f will be evaluated. y - Array of length
// m_dat. Contains the ordinate values that shall be fitted. dy - Array of
// length m_dat. Contains the standard deviations of the values y. f - A
// user-supplied parametric function f(ti;par). control - Parameter
// collection for tuning the fit procedure. In most cases, the default
// &lm_control_double is adequate. If f is only computed with
// single-precision accuracy, &lm_control_float should be used. Parameters
// are explained in lmmin2(3). status - A record used to return information
// about the minimization process: For details, see lmmin2(3).
// void lmcurve2( const int n_par, double *par, double *parerr, double *covar, const int m_dat, const double *t, const double *y, const double *dy, double (*f)( const double ti, const double *par ), const lm_control_struct *control, lm_status_struct *status);
// n_par - Number of free variables. Length of parameter vector par.
// par - Parameter vector. On input, it must contain a reasonable guess. On output, it contains the solution found to minimize ||r||.
// parerr - Parameter uncertainties vector. Array of length n_par or NULL. On output, unless it or covar is NULL, it contains the weighted parameter uncertainties for the found parameters.
// covar - Covariance matrix. Array of length n_par * n_par or NULL. On output, unless it is NULL, it contains the covariance matrix.
// m_dat - Number of data points. Length of vectors t, y, dy. Must statisfy n_par <= m_dat.
// t - Array of length m_dat. Contains the abcissae (time, or "x") for which function f will be evaluated.
// y - Array of length m_dat. Contains the ordinate values that shall be fitted.
// dy - Array of length m_dat. Contains the standard deviations of the values y.
// f - A user-supplied parametric function f(ti;par).
// control - Parameter collection for tuning the fit procedure. In most cases, the default &lm_control_double is adequate. If f is only computed with single-precision accuracy, &lm_control_float should be used. Parameters are explained in lmmin2(3).
// status - A record used to return information about the minimization process: For details, see lmmin2(3).
lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
x.size(), x.data(), y.data(), y_err.data(), aare::func::gaus,
@ -188,14 +178,12 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
// Calculate chi2
chi2 = 0;
for (ssize_t i = 0; i < y.size(); i++) {
chi2 +=
std::pow((y(i) - func::gaus(x(i), par_out.data())) / y_err(i), 2);
chi2 += std::pow((y(i) - func::gaus(x(i), par_out.data())) / y_err(i), 2);
}
}
void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out,
NDView<double, 2> chi2_out,
NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads) {
@ -209,9 +197,10 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
{par_out.shape(2)});
NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
{par_err_out.shape(2)});
fit_gaus(x, y_view, y_err_view, par_out_view, par_err_out_view,
chi2_out(row, col));
}
}
};
@ -221,8 +210,7 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
}
void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> par_out, NDView<double, 1> par_err_out,
double &chi2) {
NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2) {
// Check that we have the correct sizes
if (y.size() != x.size() || y.size() != y_err.size() ||
@ -242,14 +230,13 @@ void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
// Calculate chi2
chi2 = 0;
for (ssize_t i = 0; i < y.size(); i++) {
chi2 +=
std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
chi2 += std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
}
}
void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out,
NDView<double, 2> chi2_out, int n_threads) {
NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads) {
auto process = [&](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) {
@ -262,14 +249,15 @@ void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
{par_err_out.shape(2)});
fit_pol1(x, y_view, y_err_view, par_out_view, par_err_out_view,
chi2_out(row, col));
fit_pol1(x, y_view, y_err_view, par_out_view, par_err_out_view, chi2_out(row, col));
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks);
}
NDArray<double, 1> fit_pol1(NDView<double, 1> x, NDView<double, 1> y) {
@ -312,29 +300,27 @@ NDArray<double, 3> fit_pol1(NDView<double, 1> x, NDView<double, 3> y,
// ~~ S-CURVES ~~
// SCURVE --
std::array<double, 6> scurve_init_par(const NDView<double, 1> x,
const NDView<double, 1> y) {
// Estimate the initial parameters for the fit
std::array<double, 6> start_par{0, 0, 0, 0, 0, 0};
std::array<double, 6> scurve_init_par(const NDView<double, 1> x, const NDView<double, 1> y){
// Estimate the initial parameters for the fit
std::array<double, 6> start_par{0, 0, 0, 0, 0, 0};
auto ymax = std::max_element(y.begin(), y.end());
auto ymin = std::min_element(y.begin(), y.end());
start_par[4] = *ymin + (*ymax - *ymin) / 2;
// Find the first x where the corresponding y value is above the threshold
// (start_par[4])
for (ssize_t i = 0; i < y.size(); ++i) {
if (y[i] >= start_par[4]) {
start_par[2] = x[i];
break; // Exit the loop after finding the first valid x
auto ymax = std::max_element(y.begin(), y.end());
auto ymin = std::min_element(y.begin(), y.end());
start_par[4] = *ymin + (*ymax - *ymin) / 2;
// Find the first x where the corresponding y value is above the threshold (start_par[4])
for (ssize_t i = 0; i < y.size(); ++i) {
if (y[i] >= start_par[4]) {
start_par[2] = x[i];
break; // Exit the loop after finding the first valid x
}
}
}
start_par[3] = 2 * sqrt(start_par[2]);
start_par[0] = 100;
start_par[1] = 0.25;
start_par[5] = 1;
return start_par;
start_par[3] = 2 * sqrt(start_par[2]);
start_par[0] = 100;
start_par[1] = 0.25;
start_par[5] = 1;
return start_par;
}
// - No error
@ -348,8 +334,7 @@ NDArray<double, 1> fit_scurve(NDView<double, 1> x, NDView<double, 1> y) {
return result;
}
NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
int n_threads) {
NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y, int n_threads) {
NDArray<double, 3> result({y.shape(0), y.shape(1), 6}, 0);
auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) {
@ -373,9 +358,8 @@ NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
}
// - Error
void fit_scurve(NDView<double, 1> x, NDView<double, 1> y,
NDView<double, 1> y_err, NDView<double, 1> par_out,
NDView<double, 1> par_err_out, double &chi2) {
void fit_scurve(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2) {
// Check that we have the correct sizes
if (y.size() != x.size() || y.size() != y_err.size() ||
@ -396,15 +380,13 @@ void fit_scurve(NDView<double, 1> x, NDView<double, 1> y,
// Calculate chi2
chi2 = 0;
for (ssize_t i = 0; i < y.size(); i++) {
chi2 +=
std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
chi2 += std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
}
}
void fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
NDView<double, 3> y_err, NDView<double, 3> par_out,
NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads) {
void fit_scurve(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads) {
auto process = [&](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) {
@ -417,41 +399,40 @@ void fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
{par_err_out.shape(2)});
fit_scurve(x, y_view, y_err_view, par_out_view,
par_err_out_view, chi2_out(row, col));
fit_scurve(x, y_view, y_err_view, par_out_view, par_err_out_view, chi2_out(row, col));
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks);
}
// SCURVE2 ---
std::array<double, 6> scurve2_init_par(const NDView<double, 1> x,
const NDView<double, 1> y) {
// Estimate the initial parameters for the fit
std::array<double, 6> start_par{0, 0, 0, 0, 0, 0};
std::array<double, 6> scurve2_init_par(const NDView<double, 1> x, const NDView<double, 1> y){
// Estimate the initial parameters for the fit
std::array<double, 6> start_par{0, 0, 0, 0, 0, 0};
auto ymax = std::max_element(y.begin(), y.end());
auto ymin = std::min_element(y.begin(), y.end());
start_par[4] = *ymin + (*ymax - *ymin) / 2;
// Find the first x where the corresponding y value is above the threshold
// (start_par[4])
for (ssize_t i = 0; i < y.size(); ++i) {
if (y[i] <= start_par[4]) {
start_par[2] = x[i];
break; // Exit the loop after finding the first valid x
auto ymax = std::max_element(y.begin(), y.end());
auto ymin = std::min_element(y.begin(), y.end());
start_par[4] = *ymin + (*ymax - *ymin) / 2;
// Find the first x where the corresponding y value is above the threshold (start_par[4])
for (ssize_t i = 0; i < y.size(); ++i) {
if (y[i] <= start_par[4]) {
start_par[2] = x[i];
break; // Exit the loop after finding the first valid x
}
}
}
start_par[3] = 2 * sqrt(start_par[2]);
start_par[0] = 100;
start_par[1] = 0.25;
start_par[5] = -1;
return start_par;
start_par[3] = 2 * sqrt(start_par[2]);
start_par[0] = 100;
start_par[1] = 0.25;
start_par[5] = -1;
return start_par;
}
// - No error
@ -465,8 +446,7 @@ NDArray<double, 1> fit_scurve2(NDView<double, 1> x, NDView<double, 1> y) {
return result;
}
NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
int n_threads) {
NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, int n_threads) {
NDArray<double, 3> result({y.shape(0), y.shape(1), 6}, 0);
auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) {
@ -490,9 +470,8 @@ NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
}
// - Error
void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y,
NDView<double, 1> y_err, NDView<double, 1> par_out,
NDView<double, 1> par_err_out, double &chi2) {
void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2) {
// Check that we have the correct sizes
if (y.size() != x.size() || y.size() != y_err.size() ||
@ -513,15 +492,13 @@ void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y,
// Calculate chi2
chi2 = 0;
for (ssize_t i = 0; i < y.size(); i++) {
chi2 +=
std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
chi2 += std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
}
}
void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
NDView<double, 3> y_err, NDView<double, 3> par_out,
NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads) {
void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads) {
auto process = [&](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) {
@ -534,14 +511,15 @@ void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
{par_err_out.shape(2)});
fit_scurve2(x, y_view, y_err_view, par_out_view,
par_err_out_view, chi2_out(row, col));
fit_scurve2(x, y_view, y_err_view, par_out_view, par_err_out_view, chi2_out(row, col));
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks);
}
} // namespace aare

View File

@ -29,7 +29,8 @@ uint64_t Frame::size() const { return m_rows * m_cols; }
size_t Frame::bytes() const { return m_rows * m_cols * m_dtype.bytes(); }
std::byte *Frame::data() const { return m_data; }
std::byte *Frame::pixel_ptr(uint32_t row, uint32_t col) const {
std::byte *Frame::pixel_ptr(uint32_t row, uint32_t col) const{
if ((row >= m_rows) || (col >= m_cols)) {
std::cerr << "Invalid row or column index" << '\n';
return nullptr;
@ -37,6 +38,7 @@ std::byte *Frame::pixel_ptr(uint32_t row, uint32_t col) const {
return m_data + (row * m_cols + col) * (m_dtype.bytes());
}
Frame &Frame::operator=(Frame &&other) noexcept {
if (this == &other) {
return *this;
@ -68,4 +70,5 @@ Frame Frame::clone() const {
return frame;
}
} // namespace aare

View File

@ -65,8 +65,7 @@ TEST_CASE("Set a value in a 64 bit frame") {
// only the value we did set should be non-zero
for (size_t i = 0; i < rows; i++) {
for (size_t j = 0; j < cols; j++) {
uint64_t *data =
reinterpret_cast<uint64_t *>(frame.pixel_ptr(i, j));
uint64_t *data = reinterpret_cast<uint64_t *>(frame.pixel_ptr(i, j));
REQUIRE(data != nullptr);
if (i == 5 && j == 7) {
REQUIRE(*data == value);
@ -151,3 +150,4 @@ TEST_CASE("test explicit copy constructor") {
REQUIRE(frame2.bytes() == rows * cols * bitdepth / 8);
REQUIRE(frame2.data() != data);
}

110
src/Hdf5FileReader.test.cpp Normal file
View File

@ -0,0 +1,110 @@
/************************************************
* @file Hdf5FileReader.test.cpp
* @short test case for reading hdf5 files
***********************************************/
#include <filesystem>
#include "test_config.hpp"
#include <H5Cpp.h>
#include "aare/Hdf5FileReader.hpp"
#include "aare/NDArray.hpp"
#include <catch2/catch_all.hpp>
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
using namespace aare;
TEST_CASE("read hdf5 file", "[.hdf5file][.files]") {
// TODO generalize datasetpath
std::string filename = test_data_path() / "AngleCalibration_Test_Data" /
"ang1up_22keV_LaB60p3mm_48M_a_0320.h5";
REQUIRE(std::filesystem::exists(filename));
HDF5FileReader file_reader;
file_reader.open_file(filename);
auto dataset = file_reader.get_dataset("/entry/data/data");
auto shape = dataset.get_shape();
CHECK(shape[0] == 61440);
auto type = dataset.get_datatype();
const std::type_info *type_info = dataset.get_cpp_type();
CHECK(*type_info == typeid(uint32_t));
SECTION("read dataset into NDArray") {
NDArray<uint32_t, 1> dataset_array =
dataset.store_as_ndarray<uint32_t, 1>();
CHECK(dataset_array(0) == 866);
CHECK(dataset_array(61439) == 1436);
}
SECTION("read dataset into Frame") {
Frame frame = dataset.store_as_frame();
CHECK(*(reinterpret_cast<uint32_t *>(frame.pixel_ptr(0, 0))) == 866);
CHECK(*(reinterpret_cast<uint32_t *>(frame.pixel_ptr(0, 61439))) ==
1436);
}
SECTION("read subset of dataset") {
Frame frame(1, 10, Dtype(typeid(uint32_t)));
Subset subset{std::vector<hsize_t>{10}, std::vector<hsize_t>{10}};
dataset.read_into_buffer(frame.data(), subset);
CHECK(*(reinterpret_cast<uint32_t *>(frame.pixel_ptr(0, 0))) == 664);
CHECK(*(reinterpret_cast<uint32_t *>(frame.pixel_ptr(0, 9))) == 654);
}
/*
SECTION("read scalar") {
}
*/
}
TEST_CASE("test datatypes", "[.hdf5file]") {
auto [dtype, expected_type_info] = GENERATE(
std::make_tuple(H5::DataType(H5::PredType::NATIVE_INT), &typeid(int)),
std::make_tuple(H5::DataType(H5::PredType::NATIVE_INT8),
&typeid(int8_t)),
std::make_tuple(H5::DataType(H5::PredType::NATIVE_UINT16),
&typeid(uint16_t)),
std::make_tuple(H5::DataType(H5::PredType::NATIVE_INT16),
&typeid(int16_t)),
std::make_tuple(H5::DataType(H5::PredType::STD_U32LE),
&typeid(uint32_t)),
std::make_tuple(H5::DataType(H5::PredType::STD_I32LE),
&typeid(int32_t)),
std::make_tuple(H5::DataType(H5::PredType::NATIVE_INT32),
&typeid(int32_t)),
std::make_tuple(H5::DataType(H5::PredType::IEEE_F64LE),
&typeid(double)),
std::make_tuple(H5::DataType(H5::PredType::IEEE_F32LE), &typeid(float)),
std::make_tuple(H5::DataType(H5::PredType::NATIVE_FLOAT),
&typeid(float)),
std::make_tuple(H5::DataType(H5::PredType::NATIVE_DOUBLE),
&typeid(double)),
std::make_tuple(H5::DataType(H5::PredType::NATIVE_CHAR),
&typeid(int8_t)));
const std::type_info &type_info = deduce_cpp_type(dtype);
CHECK(type_info == *expected_type_info);
// TODO: handle bit swapping
REQUIRE_THROWS(deduce_cpp_type(
H5::DataType(H5::PredType::IEEE_F32BE))); // does not convert from big
// to little endian
}

View File

@ -19,15 +19,16 @@ JungfrauDataFile::JungfrauDataFile(const std::filesystem::path &fname) {
open_file(m_current_file_index);
}
// FileInterface
Frame JungfrauDataFile::read_frame() {
Frame JungfrauDataFile::read_frame(){
Frame f(rows(), cols(), Dtype::UINT16);
read_into(reinterpret_cast<std::byte *>(f.data()), nullptr);
return f;
}
Frame JungfrauDataFile::read_frame(size_t frame_number) {
Frame JungfrauDataFile::read_frame(size_t frame_number){
seek(frame_number);
Frame f(rows(), cols(), Dtype::UINT16);
read_into(reinterpret_cast<std::byte *>(f.data()), nullptr);
@ -36,7 +37,7 @@ Frame JungfrauDataFile::read_frame(size_t frame_number) {
std::vector<Frame> JungfrauDataFile::read_n(size_t n_frames) {
std::vector<Frame> frames;
for (size_t i = 0; i < n_frames; ++i) {
for(size_t i = 0; i < n_frames; ++i){
frames.push_back(read_frame());
}
return frames;
@ -47,7 +48,7 @@ void JungfrauDataFile::read_into(std::byte *image_buf) {
}
void JungfrauDataFile::read_into(std::byte *image_buf, size_t n_frames) {
read_into(image_buf, n_frames, nullptr);
}
}
size_t JungfrauDataFile::frame_number(size_t frame_index) {
seek(frame_index);
@ -58,9 +59,7 @@ std::array<ssize_t, 2> JungfrauDataFile::shape() const {
return {static_cast<ssize_t>(rows()), static_cast<ssize_t>(cols())};
}
DetectorType JungfrauDataFile::detector_type() const {
return DetectorType::Jungfrau;
}
DetectorType JungfrauDataFile::detector_type() const { return DetectorType::Jungfrau; }
std::string JungfrauDataFile::base_name() const { return m_base_name; }
@ -196,23 +195,22 @@ void JungfrauDataFile::read_into(std::byte *image_buf, size_t n_frames,
JungfrauDataHeader *header) {
if (header) {
for (size_t i = 0; i < n_frames; ++i)
read_into(image_buf + i * m_bytes_per_frame, header + i);
} else {
read_into(image_buf + i * m_bytes_per_frame, header + i);
}else{
for (size_t i = 0; i < n_frames; ++i)
read_into(image_buf + i * m_bytes_per_frame, nullptr);
}
}
void JungfrauDataFile::read_into(NDArray<uint16_t> *image,
JungfrauDataHeader *header) {
if (image->shape() != shape()) {
throw std::runtime_error(
LOCATION + "Image shape does not match file size: " +
std::to_string(rows()) + "x" + std::to_string(cols()));
void JungfrauDataFile::read_into(NDArray<uint16_t>* image, JungfrauDataHeader* header) {
if(image->shape()!=shape()){
throw std::runtime_error(LOCATION +
"Image shape does not match file size: " + std::to_string(rows()) + "x" + std::to_string(cols()));
}
read_into(reinterpret_cast<std::byte *>(image->data()), header);
}
JungfrauDataHeader JungfrauDataFile::read_header() {
JungfrauDataHeader header;
if (auto rc = fread(&header, 1, sizeof(header), m_fp.get());

View File

@ -1,21 +1,21 @@
#include "aare/JungfrauDataFile.hpp"
#include "test_config.hpp"
#include <catch2/catch_test_macros.hpp>
#include "test_config.hpp"
using aare::JungfrauDataFile;
using aare::JungfrauDataHeader;
TEST_CASE("Open a Jungfrau data file", "[.files]") {
// we know we have 4 files with 7, 7, 7, and 3 frames
// firs frame number if 1 and the bunch id is frame_number**2
// so we can check the header
//we know we have 4 files with 7, 7, 7, and 3 frames
//firs frame number if 1 and the bunch id is frame_number**2
//so we can check the header
auto fpath = test_data_path() / "dat" / "AldoJF500k_000000.dat";
REQUIRE(std::filesystem::exists(fpath));
JungfrauDataFile f(fpath);
REQUIRE(f.rows() == 512);
REQUIRE(f.cols() == 1024);
REQUIRE(f.bytes_per_frame() == 1048576);
REQUIRE(f.bytes_per_frame() == 1048576);
REQUIRE(f.pixels_per_frame() == 524288);
REQUIRE(f.bytes_per_pixel() == 2);
REQUIRE(f.bitdepth() == 16);
@ -25,7 +25,7 @@ TEST_CASE("Open a Jungfrau data file", "[.files]") {
REQUIRE(f.total_frames() == 24);
REQUIRE(f.current_file() == fpath);
// Check that the frame number and buch id is read correctly
//Check that the frame number and buch id is read correctly
for (size_t i = 0; i < 24; ++i) {
JungfrauDataHeader header;
aare::NDArray<uint16_t> image(f.shape());
@ -37,64 +37,65 @@ TEST_CASE("Open a Jungfrau data file", "[.files]") {
}
}
TEST_CASE("Seek in a JungfrauDataFile", "[.files]") {
TEST_CASE("Seek in a JungfrauDataFile", "[.files]"){
auto fpath = test_data_path() / "dat" / "AldoJF65k_000000.dat";
REQUIRE(std::filesystem::exists(fpath));
JungfrauDataFile f(fpath);
// The file should have 113 frames
//The file should have 113 frames
f.seek(19);
REQUIRE(f.tell() == 19);
auto h = f.read_header();
REQUIRE(h.framenum == 19 + 1);
REQUIRE(h.framenum == 19+1);
// Reading again does not change the file pointer
//Reading again does not change the file pointer
auto h2 = f.read_header();
REQUIRE(h2.framenum == 19 + 1);
REQUIRE(h2.framenum == 19+1);
f.seek(59);
REQUIRE(f.tell() == 59);
auto h3 = f.read_header();
REQUIRE(h3.framenum == 59 + 1);
REQUIRE(h3.framenum == 59+1);
JungfrauDataHeader h4;
aare::NDArray<uint16_t> image(f.shape());
f.read_into(&image, &h4);
REQUIRE(h4.framenum == 59 + 1);
REQUIRE(h4.framenum == 59+1);
// now we should be on the next frame
//now we should be on the next frame
REQUIRE(f.tell() == 60);
REQUIRE(f.read_header().framenum == 60 + 1);
REQUIRE(f.read_header().framenum == 60+1);
REQUIRE_THROWS(f.seek(86356)); // out of range
REQUIRE_THROWS(f.seek(86356)); //out of range
}
TEST_CASE("Open a Jungfrau data file with non zero file index", "[.files]") {
TEST_CASE("Open a Jungfrau data file with non zero file index", "[.files]"){
auto fpath = test_data_path() / "dat" / "AldoJF65k_000003.dat";
REQUIRE(std::filesystem::exists(fpath));
JungfrauDataFile f(fpath);
// 18 files per data file, opening the 3rd file we ignore the first 3
REQUIRE(f.total_frames() == 113 - 18 * 3);
//18 files per data file, opening the 3rd file we ignore the first 3
REQUIRE(f.total_frames() == 113-18*3);
REQUIRE(f.tell() == 0);
// Frame numbers start at 1 in the first file
REQUIRE(f.read_header().framenum == 18 * 3 + 1);
//Frame numbers start at 1 in the first file
REQUIRE(f.read_header().framenum == 18*3+1);
// moving relative to the third file
f.seek(5);
REQUIRE(f.read_header().framenum == 18 * 3 + 1 + 5);
REQUIRE(f.read_header().framenum == 18*3+1+5);
// ignoring the first 3 files
REQUIRE(f.n_files() == 4);
REQUIRE(f.current_file().stem() == "AldoJF65k_000003");
}
TEST_CASE("Read into throws if size doesn't match", "[.files]") {
TEST_CASE("Read into throws if size doesn't match", "[.files]"){
auto fpath = test_data_path() / "dat" / "AldoJF65k_000000.dat";
REQUIRE(std::filesystem::exists(fpath));
@ -108,4 +109,6 @@ TEST_CASE("Read into throws if size doesn't match", "[.files]") {
REQUIRE_THROWS(f.read_into(&image));
REQUIRE(f.tell() == 0);
}

View File

@ -0,0 +1,33 @@
/************************************************
* @file MythenFileReader.test.cpp
* @short test case for angle calibration class
***********************************************/
#include "aare/MythenFileReader.hpp"
#include <filesystem>
#include "test_config.hpp"
#include <catch2/catch_all.hpp>
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
using namespace aare;
TEST_CASE("test mythenfile_reader", "[.mythenfilereader][.files]") {
auto fpath = test_data_path() / "AngleCalibration_Test_Data";
REQUIRE(std::filesystem::exists(fpath));
MythenFileReader file_reader(fpath, "ang1up_22keV_LaB60p3mm_48M_a_0");
auto frame = file_reader.read_frame(320);
CHECK(frame.detector_angle == 0.99955);
CHECK(frame.channel_mask == std::array<uint8_t, 3>{0, 0, 1});
CHECK(frame.photon_counts.size() == 61440);
}

View File

@ -35,7 +35,7 @@ TEST_CASE("Construct from an NDView") {
}
}
TEST_CASE("3D NDArray from NDView") {
TEST_CASE("3D NDArray from NDView"){
std::vector<int> data(27);
std::iota(data.begin(), data.end(), 0);
NDView<int, 3> view(data.data(), Shape<3>{3, 3, 3});
@ -44,9 +44,9 @@ TEST_CASE("3D NDArray from NDView") {
REQUIRE(image.size() == view.size());
REQUIRE(image.data() != view.data());
for (ssize_t i = 0; i < image.shape(0); i++) {
for (ssize_t j = 0; j < image.shape(1); j++) {
for (ssize_t k = 0; k < image.shape(2); k++) {
for(ssize_t i=0; i<image.shape(0); i++){
for(ssize_t j=0; j<image.shape(1); j++){
for(ssize_t k=0; k<image.shape(2); k++){
REQUIRE(image(i, j, k) == view(i, j, k));
}
}
@ -132,7 +132,7 @@ TEST_CASE("Elementwise multiplication of 3D image") {
NDArray<int> MultiplyNDArrayUsingOperator(NDArray<int> &a, NDArray<int> &b) {
// return a * a * b * b;
NDArray<int> c = a * b;
NDArray<int>c = a*b;
return c;
}
@ -162,6 +162,7 @@ NDArray<int> AddNDArrayUsingIndex(NDArray<int> &a, NDArray<int> &b) {
return res;
}
TEST_CASE("Compare two images") {
NDArray<int> a;
NDArray<int> b;
@ -221,6 +222,7 @@ TEST_CASE("Bitwise and on data") {
REQUIRE(a(2) == 384);
}
TEST_CASE("Elementwise operations on images") {
std::array<ssize_t, 2> shape{5, 5};
double a_val = 3.0;
@ -256,8 +258,7 @@ TEST_CASE("Elementwise operations on images") {
NDArray<double> A(shape, a_val);
NDArray<double> B(shape, b_val);
NDArray<double> C = A - B;
// auto C = A - B; // This works but the result is a lazy ArraySub
// object
// auto C = A - B; // This works but the result is a lazy ArraySub object
// Value of C matches
for (uint32_t i = 0; i < C.size(); ++i) {
@ -281,8 +282,7 @@ TEST_CASE("Elementwise operations on images") {
SECTION("Multiply two images") {
NDArray<double> A(shape, a_val);
NDArray<double> B(shape, b_val);
// auto C = A * B; // This works but the result is a lazy ArrayMul
// object
// auto C = A * B; // This works but the result is a lazy ArrayMul object
NDArray<double> C = A * B;
// Value of C matches
@ -307,8 +307,7 @@ TEST_CASE("Elementwise operations on images") {
SECTION("Divide two images") {
NDArray<double> A(shape, a_val);
NDArray<double> B(shape, b_val);
// auto C = A / B; // This works but the result is a lazy ArrayDiv
// object
// auto C = A / B; // This works but the result is a lazy ArrayDiv object
NDArray<double> C = A / B;
// Value of C matches

View File

@ -2,8 +2,8 @@
#include <catch2/catch_test_macros.hpp>
#include <iostream>
#include <numeric>
#include <vector>
#include <numeric>
using aare::NDView;
using aare::Shape;
@ -151,10 +151,8 @@ TEST_CASE("divide with another span") {
std::vector<int> vec1{3, 2, 1};
std::vector<int> result{3, 6, 3};
NDView<int, 1> data0(vec0.data(),
Shape<1>{static_cast<ssize_t>(vec0.size())});
NDView<int, 1> data1(vec1.data(),
Shape<1>{static_cast<ssize_t>(vec1.size())});
NDView<int, 1> data0(vec0.data(), Shape<1>{static_cast<ssize_t>(vec0.size())});
NDView<int, 1> data1(vec1.data(), Shape<1>{static_cast<ssize_t>(vec1.size())});
data0 /= data1;
@ -183,7 +181,8 @@ TEST_CASE("compare two views") {
REQUIRE((view1 == view2));
}
TEST_CASE("Create a view over a vector") {
TEST_CASE("Create a view over a vector"){
std::vector<int> vec(12);
std::iota(vec.begin(), vec.end(), 0);
auto v = aare::make_view(vec);

View File

@ -4,16 +4,16 @@
namespace aare {
NumpyFile::NumpyFile(const std::filesystem::path &fname,
const std::string &mode, FileConfig cfg) {
NumpyFile::NumpyFile(const std::filesystem::path &fname, const std::string &mode, FileConfig cfg) {
// TODO! add opts to constructor
m_mode = mode;
if (mode == "r") {
fp = fopen(fname.string().c_str(), "rb");
if (!fp) {
throw std::runtime_error(
fmt::format("Could not open: {} for reading", fname.string()));
throw std::runtime_error(fmt::format("Could not open: {} for reading", fname.string()));
}
load_metadata();
} else if (mode == "w") {
@ -24,15 +24,11 @@ NumpyFile::NumpyFile(const std::filesystem::path &fname,
m_header.shape = {0, cfg.rows, cfg.cols};
fp = fopen(fname.string().c_str(), "wb");
if (!fp) {
throw std::runtime_error(
fmt::format("Could not open: {} for reading", fname.string()));
throw std::runtime_error(fmt::format("Could not open: {} for reading", fname.string()));
}
initial_header_len = aare::NumpyHelpers::write_header(
std::filesystem::path(fname.c_str()), m_header);
initial_header_len = aare::NumpyHelpers::write_header(std::filesystem::path(fname.c_str()), m_header);
}
m_pixels_per_frame =
std::accumulate(m_header.shape.begin() + 1, m_header.shape.end(), 1,
std::multiplies<>());
m_pixels_per_frame = std::accumulate(m_header.shape.begin() + 1, m_header.shape.end(), 1, std::multiplies<>());
m_bytes_per_frame = m_header.dtype.bitdepth() / 8 * m_pixels_per_frame;
}
@ -67,8 +63,7 @@ void NumpyFile::get_frame_into(size_t frame_number, std::byte *image_buf) {
if (frame_number > m_header.shape[0]) {
throw std::invalid_argument("Frame number out of range");
}
if (fseek(fp, header_size + frame_number * m_bytes_per_frame,
SEEK_SET)) // NOLINT
if (fseek(fp, header_size + frame_number * m_bytes_per_frame, SEEK_SET)) // NOLINT
throw std::runtime_error("Could not seek to frame");
size_t const rc = fread(image_buf, m_bytes_per_frame, 1, fp);
@ -118,8 +113,7 @@ NumpyFile::~NumpyFile() noexcept {
// write header
size_t const rc = fwrite(header_str.c_str(), header_str.size(), 1, fp);
if (rc != 1) {
std::cout << "Error writing header to numpy file in destructor"
<< std::endl;
std::cout << "Error writing header to numpy file in destructor" << std::endl;
}
}
@ -146,10 +140,8 @@ void NumpyFile::load_metadata() {
}
// read version
rc =
fread(reinterpret_cast<char *>(&major_ver_), sizeof(major_ver_), 1, fp);
rc +=
fread(reinterpret_cast<char *>(&minor_ver_), sizeof(minor_ver_), 1, fp);
rc = fread(reinterpret_cast<char *>(&major_ver_), sizeof(major_ver_), 1, fp);
rc += fread(reinterpret_cast<char *>(&minor_ver_), sizeof(minor_ver_), 1, fp);
if (rc != 2) {
throw std::runtime_error("Error reading numpy version");
}
@ -167,8 +159,7 @@ void NumpyFile::load_metadata() {
if (rc != 1) {
throw std::runtime_error("Error reading header length");
}
header_size = aare::NumpyHelpers::magic_string_length + 2 +
header_len_size + header_len;
header_size = aare::NumpyHelpers::magic_string_length + 2 + header_len_size + header_len;
if (header_size % 16 != 0) {
fmt::print("Warning: header length is not a multiple of 16\n");
}

View File

@ -1,8 +1,8 @@
#include "aare/NumpyFile.hpp"
#include "aare/NDArray.hpp"
#include "test_config.hpp"
#include <catch2/catch_test_macros.hpp>
#include "test_config.hpp"
using aare::Dtype;
using aare::NumpyFile;
@ -23,7 +23,7 @@ TEST_CASE("Read a 1D numpy file with int32 data type", "[.integration]") {
REQUIRE(data(i) == i);
}
}
TEST_CASE("Read a 3D numpy file with np.double data type", "[.integration]") {
auto fpath = test_data_path() / "numpy" / "test_3d_double.npy";

View File

@ -29,8 +29,7 @@ namespace aare {
std::string NumpyHeader::to_string() const {
std::stringstream sstm;
sstm << "dtype: " << dtype.to_string()
<< ", fortran_order: " << fortran_order << ' ';
sstm << "dtype: " << dtype.to_string() << ", fortran_order: " << fortran_order << ' ';
sstm << "shape: (";
for (auto item : shape)
sstm << item << ',';
@ -38,10 +37,10 @@ std::string NumpyHeader::to_string() const {
return sstm.str();
}
namespace NumpyHelpers {
std::unordered_map<std::string, std::string>
parse_dict(std::string in, const std::vector<std::string> &keys) {
std::unordered_map<std::string, std::string> parse_dict(std::string in, const std::vector<std::string> &keys) {
std::unordered_map<std::string, std::string> map;
if (keys.empty())
return map;
@ -101,8 +100,7 @@ aare::Dtype parse_descr(std::string typestring) {
constexpr char little_endian_char = '<';
constexpr char big_endian_char = '>';
constexpr char no_endian_char = '|';
constexpr std::array<char, 3> endian_chars = {
little_endian_char, big_endian_char, no_endian_char};
constexpr std::array<char, 3> endian_chars = {little_endian_char, big_endian_char, no_endian_char};
constexpr std::array<char, 4> numtype_chars = {'f', 'i', 'u', 'c'};
const char byteorder_c = typestring[0];
@ -141,9 +139,7 @@ std::string get_value_from_map(const std::string &mapstr) {
return trim(tmp);
}
bool is_digits(const std::string &str) {
return std::all_of(str.begin(), str.end(), ::isdigit);
}
bool is_digits(const std::string &str) { return std::all_of(str.begin(), str.end(), ::isdigit); }
std::vector<std::string> parse_tuple(std::string in) {
std::vector<std::string> v;
@ -219,25 +215,20 @@ inline std::string write_boolean(bool b) {
return "False";
}
inline std::string write_header_dict(const std::string &descr,
bool fortran_order,
const std::vector<size_t> &shape) {
inline std::string write_header_dict(const std::string &descr, bool fortran_order, const std::vector<size_t> &shape) {
std::string const s_fortran_order = write_boolean(fortran_order);
std::string const shape_s = write_tuple(shape);
return "{'descr': '" + descr + "', 'fortran_order': " + s_fortran_order +
", 'shape': " + shape_s + ", }";
return "{'descr': '" + descr + "', 'fortran_order': " + s_fortran_order + ", 'shape': " + shape_s + ", }";
}
size_t write_header(const std::filesystem::path &fname,
const NumpyHeader &header) {
size_t write_header(const std::filesystem::path &fname, const NumpyHeader &header) {
std::ofstream out(fname, std::ios::binary | std::ios::out);
return write_header(out, header);
}
size_t write_header(std::ostream &out, const NumpyHeader &header) {
std::string const header_dict = write_header_dict(
header.dtype.to_string(), header.fortran_order, header.shape);
std::string const header_dict = write_header_dict(header.dtype.to_string(), header.fortran_order, header.shape);
size_t length = magic_string_length + 2 + 2 + header_dict.length() + 1;
@ -256,22 +247,17 @@ size_t write_header(std::ostream &out, const NumpyHeader &header) {
// write header length
if (version_major == 1 && version_minor == 0) {
auto header_len =
static_cast<uint16_t>(header_dict.length() + padding.length() + 1);
auto header_len = static_cast<uint16_t>(header_dict.length() + padding.length() + 1);
std::array<uint8_t, 2> header_len_le16{
static_cast<uint8_t>((header_len >> 0) & 0xff),
static_cast<uint8_t>((header_len >> 8) & 0xff)};
std::array<uint8_t, 2> header_len_le16{static_cast<uint8_t>((header_len >> 0) & 0xff),
static_cast<uint8_t>((header_len >> 8) & 0xff)};
out.write(reinterpret_cast<char *>(header_len_le16.data()), 2);
} else {
auto header_len =
static_cast<uint32_t>(header_dict.length() + padding.length() + 1);
auto header_len = static_cast<uint32_t>(header_dict.length() + padding.length() + 1);
std::array<uint8_t, 4> header_len_le32{
static_cast<uint8_t>((header_len >> 0) & 0xff),
static_cast<uint8_t>((header_len >> 8) & 0xff),
static_cast<uint8_t>((header_len >> 16) & 0xff),
static_cast<uint8_t>((header_len >> 24) & 0xff)};
static_cast<uint8_t>((header_len >> 0) & 0xff), static_cast<uint8_t>((header_len >> 8) & 0xff),
static_cast<uint8_t>((header_len >> 16) & 0xff), static_cast<uint8_t>((header_len >> 24) & 0xff)};
out.write(reinterpret_cast<char *>(header_len_le32.data()), 4);
}

View File

@ -19,9 +19,7 @@ TEST_CASE("Check for quotes and return stripped string") {
REQUIRE(parse_str("''") == "");
}
TEST_CASE("parsing a string without quotes throws") {
REQUIRE_THROWS(parse_str("hej"));
}
TEST_CASE("parsing a string without quotes throws") { REQUIRE_THROWS(parse_str("hej")); }
TEST_CASE("trim whitespace") {
REQUIRE(trim(" hej ") == "hej");
@ -55,8 +53,7 @@ TEST_CASE("is element in array") {
}
TEST_CASE("Parse numpy dict") {
std::string in =
"{'descr': '<f4', 'fortran_order': False, 'shape': (3, 4)}";
std::string in = "{'descr': '<f4', 'fortran_order': False, 'shape': (3, 4)}";
std::vector<std::string> keys{"descr", "fortran_order", "shape"};
auto map = parse_dict(in, keys);
REQUIRE(map["descr"] == "'<f4'");

View File

@ -1,7 +1,8 @@
#include "aare/Pedestal.hpp"
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <catch2/catch_test_macros.hpp>
#include <chrono>
#include <random>
@ -57,8 +58,7 @@ TEST_CASE("test pedestal push") {
if (k < 5) {
REQUIRE(pedestal.cur_samples()(i, j) == k + 1);
REQUIRE(pedestal.get_sum()(i, j) == (k + 1) * (i + j));
REQUIRE(pedestal.get_sum2()(i, j) ==
(k + 1) * (i + j) * (i + j));
REQUIRE(pedestal.get_sum2()(i, j) == (k + 1) * (i + j) * (i + j));
} else {
REQUIRE(pedestal.cur_samples()(i, j) == 5);
REQUIRE(pedestal.get_sum()(i, j) == 5 * (i + j));
@ -95,12 +95,9 @@ TEST_CASE("test pedestal with normal distribution") {
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 5; j++) {
REQUIRE_THAT(mean(i, j),
Catch::Matchers::WithinAbs(MEAN, MEAN * TOLERANCE));
REQUIRE_THAT(variance(i, j),
Catch::Matchers::WithinAbs(VAR, VAR * TOLERANCE));
REQUIRE_THAT(standard_deviation(i, j),
Catch::Matchers::WithinAbs(STD, STD * TOLERANCE));
REQUIRE_THAT(mean(i, j), Catch::Matchers::WithinAbs(MEAN, MEAN * TOLERANCE));
REQUIRE_THAT(variance(i, j), Catch::Matchers::WithinAbs(VAR, VAR * TOLERANCE));
REQUIRE_THAT(standard_deviation(i, j), Catch::Matchers::WithinAbs(STD, STD * TOLERANCE));
}
}
}

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