diff --git a/include/aare/ClusterFile.hpp b/include/aare/ClusterFile.hpp index 1995a16..ff5d338 100644 --- a/include/aare/ClusterFile.hpp +++ b/include/aare/ClusterFile.hpp @@ -11,7 +11,6 @@ namespace aare { -<<<<<<< HEAD /* Binary cluster file. Expects data to be layed out as: int32_t frame_number @@ -21,44 +20,6 @@ int32_t frame_number uint32_t number_of_clusters .... */ -======= - -// TODO! Legacy enums, migrate to enum class -typedef enum { - cBottomLeft = 0, - cBottomRight = 1, - cTopLeft = 2, - cTopRight = 3 -} corner; - -typedef enum { - pBottomLeft = 0, - pBottom = 1, - pBottomRight = 2, - pLeft = 3, - pCenter = 4, - pRight = 5, - pTopLeft = 6, - pTop = 7, - pTopRight = 8 -} pixel; - -struct Eta2 { - double x; - double y; - corner c; - int32_t sum; -}; - -struct ClusterAnalysis { - uint32_t c; - int32_t tot; - double etax; - double etay; -}; - - ->>>>>>> developer // TODO: change to support any type of clusters, e.g. header line with // clsuter_size_x, cluster_size_y, @@ -109,8 +70,6 @@ class ClusterFile { */ ClusterVector read_clusters(size_t n_clusters); - ClusterVector read_clusters(size_t n_clusters, ROI roi); - /** * @brief Read a single frame from the file and return the clusters. The * cluster vector will have the frame number set. diff --git a/python/aare/__init__.py b/python/aare/__init__.py index 36aac14..8c51d73 100644 --- a/python/aare/__init__.py +++ b/python/aare/__init__.py @@ -3,7 +3,7 @@ from . import _aare from ._aare import File, RawMasterFile, RawSubFile, JungfrauDataFile -from ._aare import Pedestal_d, Pedestal_f, ClusterFinder, VarClusterFinder +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 diff --git a/python/src/cluster_file.hpp b/python/src/cluster_file.hpp index 7ece8e6..3e7aa48 100644 --- a/python/src/cluster_file.hpp +++ b/python/src/cluster_file.hpp @@ -19,11 +19,12 @@ namespace py = pybind11; using namespace ::aare; -template +template void define_cluster_file_io_bindings(py::module &m, const std::string &typestr) { - // PYBIND11_NUMPY_DTYPE(Cluster, x, y, - // data); // is this used - maybe use as cluster type + + using ClusterType = Cluster; auto class_name = fmt::format("ClusterFile_{}", typestr); diff --git a/python/src/module.cpp b/python/src/module.cpp index 78fd283..8d5b5ab 100644 --- a/python/src/module.cpp +++ b/python/src/module.cpp @@ -32,12 +32,12 @@ PYBIND11_MODULE(_aare, m) { define_interpolation_bindings(m); define_jungfrau_data_file_io_bindings(m); - define_cluster_file_io_bindings>(m, "Cluster3x3i"); - define_cluster_file_io_bindings>(m, "Cluster3x3d"); - define_cluster_file_io_bindings>(m, "Cluster3x3f"); - define_cluster_file_io_bindings>(m, "Cluster2x2i"); - define_cluster_file_io_bindings>(m, "Cluster2x2f"); - define_cluster_file_io_bindings>(m, "Cluster2x2d"); + define_cluster_file_io_bindings(m, "Cluster3x3i"); + define_cluster_file_io_bindings(m, "Cluster3x3d"); + define_cluster_file_io_bindings(m, "Cluster3x3f"); + define_cluster_file_io_bindings(m, "Cluster2x2i"); + define_cluster_file_io_bindings(m, "Cluster2x2f"); + define_cluster_file_io_bindings(m, "Cluster2x2d"); define_cluster_vector(m, "Cluster3x3i"); define_cluster_vector(m, "Cluster3x3d"); diff --git a/python/tests/conftest.py b/python/tests/conftest.py index 5badf13..fbcfeb3 100644 --- a/python/tests/conftest.py +++ b/python/tests/conftest.py @@ -25,5 +25,10 @@ def pytest_collection_modifyitems(config, items): @pytest.fixture def test_data_path(): - return Path(os.environ["AARE_TEST_DATA"]) + env_value = os.environ.get("AARE_TEST_DATA") + if not env_value: + raise RuntimeError("Environment variable AARE_TEST_DATA is not set or is empty") + + return Path(env_value) + diff --git a/python/tests/test_Cluster.py b/python/tests/test_Cluster.py index 3bb4828..29d5ad9 100644 --- a/python/tests/test_Cluster.py +++ b/python/tests/test_Cluster.py @@ -1,7 +1,9 @@ import pytest import numpy as np -import aare._aare as aare #import ClusterVector_Cluster3x3i, ClusterVector_Cluster2x2i, Interpolator, Cluster3x3i, ClusterFinder_Cluster3x3i, Cluster2x2i, ClusterFile_Cluster3x3i, Cluster3x3f, calculate_eta2 +import aare._aare as aare +from conftest import test_data_path + def test_ClusterVector(): """Test ClusterVector""" @@ -64,15 +66,19 @@ def test_Interpolator(): assert interpolated_photons[0]["energy"] == 4 @pytest.mark.files -def test_cluster_file(): +def test_cluster_file(test_data_path): """Test ClusterFile""" - cluster_file = aare.ClusterFile_Cluster3x3i(test_data_path() / "clust/single_frame_97_clustrers.clust") - clustervector = cluster_file.read_clusters() #conversion does not work + cluster_file = aare.ClusterFile_Cluster3x3i(test_data_path / "clust/single_frame_97_clustrers.clust") + clustervector = cluster_file.read_clusters(10) #conversion does not work cluster_file.close() + assert clustervector.size == 10 + ###reading with wrong file - cluster_file = ClusterFile_Cluster2x2i(test_data_path() / "clust/single_frame_97_clustrers.clust") #TODO check behavior! + with pytest.raises(TypeError): + cluster_file = aare.ClusterFile_Cluster2x2i(test_data_path / "clust/single_frame_97_clustrers.clust") + cluster_file.close() def test_calculate_eta(): """Calculate Eta""" @@ -103,6 +109,7 @@ def test_cluster_finder(): assert clusters.size == 0 +#TODO dont understand behavior def test_cluster_collector(): """Test ClusterCollector""" diff --git a/src/ClusterFile.cpp b/src/ClusterFile.cpp deleted file mode 100644 index f77ac92..0000000 --- a/src/ClusterFile.cpp +++ /dev/null @@ -1,396 +0,0 @@ -#include "aare/ClusterFile.hpp" - -#include - -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 noise_map){ - m_noise_map = NDArray(noise_map); -} - -void ClusterFile::set_gain_map(const NDView gain_map){ - m_gain_map = NDArray(gain_map); -} - -ClusterFile::~ClusterFile() { close(); } - -void ClusterFile::close() { - if (fp) { - fclose(fp); - fp = nullptr; - } -} - -void ClusterFile::write_frame(const ClusterVector &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 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 ClusterFile::read_clusters_without_cut(size_t n_clusters) { - if (m_mode != "r") { - throw std::runtime_error("File not opened for reading"); - } - - ClusterVector 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(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 ClusterFile::read_clusters_with_cut(size_t n_clusters) { - ClusterVector 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(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(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 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 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 clusters(3, 3, n_clusters); - clusters.set_frame_number(frame_number); - - if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) != - static_cast(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 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 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(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 calculate_eta2(ClusterVector &clusters) { - //TOTO! make work with 2x2 clusters - NDArray eta2({static_cast(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(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(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 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(cl.data[4]) / (cl.data[3] + cl.data[4]); - if ((cl.data[1] + cl.data[4]) != 0) - eta.y = - static_cast(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(cl.data[5]) / (cl.data[4] + cl.data[5]); - if ((cl.data[1] + cl.data[4]) != 0) - eta.y = - static_cast(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(cl.data[4]) / (cl.data[3] + cl.data[4]); - if ((cl.data[7] + cl.data[4]) != 0) - eta.y = - static_cast(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(cl.data[5]) / (cl.data[5] + cl.data[4]); - if ((cl.data[7] + cl.data[4]) != 0) - eta.y = - static_cast(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(cl.data[1]) / (cl.data[0] + cl.data[1]); - if ((cl.data[0] + cl.data[2]) != 0) - eta.y = static_cast(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 \ No newline at end of file diff --git a/src/Interpolator.cpp b/src/Interpolator.cpp index ecc3e61..4bc2b34 100644 --- a/src/Interpolator.cpp +++ b/src/Interpolator.cpp @@ -53,90 +53,4 @@ Interpolator::Interpolator(NDView etacube, NDView xbins, } } -<<<<<<< HEAD -======= -std::vector Interpolator::interpolate(const ClusterVector& clusters) { - std::vector photons; - photons.reserve(clusters.size()); - - if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) { - for (size_t i = 0; i(i); - Eta2 eta= calculate_eta2(cluster); - - Photon photon; - photon.x = cluster.x; - photon.y = cluster.y; - photon.energy = eta.sum; - - - //Finding the index of the last element that is smaller - //should work fine as long as we have many bins - auto ie = last_smaller(m_energy_bins, photon.energy); - auto ix = last_smaller(m_etabinsx, eta.x); - auto iy = last_smaller(m_etabinsy, eta.y); - - double dX{}, dY{}; - // cBottomLeft = 0, - // cBottomRight = 1, - // cTopLeft = 2, - // cTopRight = 3 - switch (eta.c) { - case cTopLeft: - dX = -1.; - dY = 0.; - break; - case cTopRight:; - dX = 0.; - dY = 0.; - break; - case cBottomLeft: - dX = -1.; - dY = -1.; - break; - case cBottomRight: - dX = 0.; - dY = -1.; - break; - } - photon.x += m_ietax(ix, iy, ie)*2 + dX; - photon.y += m_ietay(ix, iy, ie)*2 + dY; - photons.push_back(photon); - } - }else if(clusters.cluster_size_x() == 2 || clusters.cluster_size_y() == 2){ - for (size_t i = 0; i(i); - Eta2 eta= calculate_eta2(cluster); - - Photon photon; - photon.x = cluster.x; - photon.y = cluster.y; - photon.energy = eta.sum; - - //Now do some actual interpolation. - //Find which energy bin the cluster is in - // auto ie = nearest_index(m_energy_bins, photon.energy)-1; - // auto ix = nearest_index(m_etabinsx, eta.x)-1; - // auto iy = nearest_index(m_etabinsy, eta.y)-1; - //Finding the index of the last element that is smaller - //should work fine as long as we have many bins - auto ie = last_smaller(m_energy_bins, photon.energy); - auto ix = last_smaller(m_etabinsx, eta.x); - auto iy = last_smaller(m_etabinsy, eta.y); - - photon.x += m_ietax(ix, iy, ie)*2; //eta goes between 0 and 1 but we could move the hit anywhere in the 2x2 - photon.y += m_ietay(ix, iy, ie)*2; - photons.push_back(photon); - } - - }else{ - throw std::runtime_error("Only 3x3 and 2x2 clusters are supported for interpolation"); - } - - - return photons; -} - ->>>>>>> developer } // namespace aare \ No newline at end of file