solved merge conflict

src/ClusterFile.cpp

@@ -0,0 +1,396 @@
+#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);
+}
+
+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

src/ClusterFile.test.cpp

@@ -8,10 +8,9 @@
 using aare::Cluster;
 using aare::ClusterFile;
 
-TEST_CASE("Read one frame from a a cluster file", "[.integration]") {
+TEST_CASE("Read one frame from a a cluster file", "[.files]") {
     // We know that the frame has 97 clusters
-    auto fpath =
-        test_data_path() / "clusters" / "single_frame_97_clustrers.clust";
+    auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust";
     REQUIRE(std::filesystem::exists(fpath));
 
     ClusterFile<Cluster<int32_t, 3, 3>> f(fpath);
@@ -20,10 +19,9 @@ TEST_CASE("Read one frame from a a cluster file", "[.integration]") {
     REQUIRE(clusters.frame_number() == 135);
 }
 
-TEST_CASE("Read one frame using ROI", "[.integration]") {
+TEST_CASE("Read one frame using ROI", "[.files]") {
     // We know that the frame has 97 clusters
-    auto fpath =
-        test_data_path() / "clusters" / "single_frame_97_clustrers.clust";
+    auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust";
     REQUIRE(std::filesystem::exists(fpath));
 
     ClusterFile<Cluster<int32_t, 3, 3>> f(fpath);
@@ -47,10 +45,9 @@ TEST_CASE("Read one frame using ROI", "[.integration]") {
     }
 }
 
-TEST_CASE("Read clusters from single frame file", "[.integration]") {
+TEST_CASE("Read clusters from single frame file", "[.files]") {
 
-    auto fpath =
-        test_data_path() / "clusters" / "single_frame_97_clustrers.clust";
+    auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust";
     REQUIRE(std::filesystem::exists(fpath));
 
     SECTION("Read fewer clusters than available") {

src/File.cpp

@@ -1,4 +1,5 @@
 #include "aare/File.hpp"
+#include "aare/JungfrauDataFile.hpp"
 #include "aare/NumpyFile.hpp"
 #include "aare/RawFile.hpp"
 
@@ -27,6 +28,8 @@ File::File(const std::filesystem::path &fname, const std::string &mode,
     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"){
+        file_impl = std::make_unique<JungfrauDataFile>(fname);
     } else {
         throw std::runtime_error("Unsupported file type");
     }

src/FilePtr.cpp

@@ -0,0 +1,44 @@
+
+#include "aare/FilePtr.hpp"
+#include <fmt/format.h>
+#include <stdexcept>
+#include <utility>
+
+namespace aare {
+
+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()));
+}
+
+FilePtr::FilePtr(FilePtr &&other) { std::swap(fp_, other.fp_); }
+
+FilePtr &FilePtr::operator=(FilePtr &&other) {
+    std::swap(fp_, other.fp_);
+    return *this;
+}
+
+FILE *FilePtr::get() { return fp_; }
+
+int64_t FilePtr::tell() {
+    auto pos = ftell(fp_);
+    if (pos == -1)
+        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(){
+    if (feof(fp_)) {
+        return "End of file reached";
+    }
+    if (ferror(fp_)) {
+        return fmt::format("Error reading file: {}", std::strerror(errno));
+    }
+    return "";
+}
+} // namespace aare

src/Fit.cpp

@@ -18,7 +18,7 @@ double gaus(const double x, const double *par) {
 
 NDArray<double, 1> gaus(NDView<double, 1> x, NDView<double, 1> par) {
     NDArray<double, 1> y({x.shape(0)}, 0);
-    for (size_t i = 0; i < x.size(); i++) {
+    for (ssize_t i = 0; i < x.size(); i++) {
         y(i) = gaus(x(i), par.data());
     }
     return y;
@@ -28,7 +28,7 @@ double pol1(const double x, const double *par) { return par[0] * x + par[1]; }
 
 NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par) {
     NDArray<double, 1> y({x.shape()}, 0);
-    for (size_t i = 0; i < x.size(); i++) {
+    for (ssize_t i = 0; i < x.size(); i++) {
         y(i) = pol1(x(i), par.data());
     }
     return y;
@@ -153,7 +153,7 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
 
     // Calculate chi2
     chi2 = 0;
-    for (size_t i = 0; i < y.size(); i++) {
+    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);
     }
 }
@@ -205,7 +205,7 @@ void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
 
     // Calculate chi2
     chi2 = 0;
-    for (size_t i = 0; i < y.size(); i++) {
+    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);
     }
 }

src/Interpolator.cpp

@@ -53,4 +53,90 @@ Interpolator::Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
     }
 }
 
+<<<<<<< HEAD
+=======
+std::vector<Photon> Interpolator::interpolate(const ClusterVector<int32_t>& clusters) {
+    std::vector<Photon> photons;
+    photons.reserve(clusters.size());
+
+    if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) {
+        for (size_t i = 0; i<clusters.size(); i++){
+
+            auto cluster = clusters.at<Cluster3x3>(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<clusters.size(); i++){
+            auto cluster = clusters.at<Cluster2x2>(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

src/JungfrauDataFile.cpp

@@ -0,0 +1,238 @@
+#include "aare/JungfrauDataFile.hpp"
+#include "aare/algorithm.hpp"
+#include "aare/defs.hpp"
+
+#include <cerrno>
+#include <fmt/format.h>
+
+namespace aare {
+
+JungfrauDataFile::JungfrauDataFile(const std::filesystem::path &fname) {
+
+    if (!std::filesystem::exists(fname)) {
+        throw std::runtime_error(LOCATION +
+                                 "File does not exist: " + fname.string());
+    }
+    find_frame_size(fname);
+    parse_fname(fname);
+    scan_files();
+    open_file(m_current_file_index);
+}
+
+
+// FileInterface
+
+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){
+    seek(frame_number);
+    Frame f(rows(), cols(), Dtype::UINT16);
+    read_into(reinterpret_cast<std::byte *>(f.data()), nullptr);
+    return f;
+}
+
+std::vector<Frame> JungfrauDataFile::read_n(size_t n_frames) {
+    std::vector<Frame> frames;
+    for(size_t i = 0; i < n_frames; ++i){
+        frames.push_back(read_frame());
+    }
+    return frames;
+}
+
+void JungfrauDataFile::read_into(std::byte *image_buf) {
+    read_into(image_buf, nullptr);
+}
+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);
+    return read_header().framenum;
+}
+
+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; }
+
+std::string JungfrauDataFile::base_name() const { return m_base_name; }
+
+size_t JungfrauDataFile::bytes_per_frame() { return m_bytes_per_frame; }
+
+size_t JungfrauDataFile::pixels_per_frame() { return m_rows * m_cols; }
+
+size_t JungfrauDataFile::bytes_per_pixel() const { return sizeof(pixel_type); }
+
+size_t JungfrauDataFile::bitdepth() const {
+    return bytes_per_pixel() * bits_per_byte;
+}
+
+void JungfrauDataFile::seek(size_t frame_index) {
+    if (frame_index >= m_total_frames) {
+        throw std::runtime_error(LOCATION + "Frame index out of range: " +
+                                 std::to_string(frame_index));
+    }
+    m_current_frame_index = frame_index;
+    auto file_index = first_larger(m_last_frame_in_file, frame_index);
+
+    if (file_index != m_current_file_index)
+        open_file(file_index);
+
+    auto frame_offset = (file_index)
+                            ? frame_index - m_last_frame_in_file[file_index - 1]
+                            : frame_index;
+    auto byte_offset = frame_offset * (m_bytes_per_frame + header_size);
+    m_fp.seek(byte_offset);
+};
+
+size_t JungfrauDataFile::tell() { return m_current_frame_index; }
+size_t JungfrauDataFile::total_frames() const { return m_total_frames; }
+size_t JungfrauDataFile::rows() const { return m_rows; }
+size_t JungfrauDataFile::cols() const { return m_cols; }
+
+size_t JungfrauDataFile::n_files() const { return m_last_frame_in_file.size(); }
+
+void JungfrauDataFile::find_frame_size(const std::filesystem::path &fname) {
+
+    static constexpr size_t module_data_size =
+        header_size + sizeof(pixel_type) * 512 * 1024;
+    static constexpr size_t half_data_size =
+        header_size + sizeof(pixel_type) * 256 * 1024;
+    static constexpr size_t chip_data_size =
+        header_size + sizeof(pixel_type) * 256 * 256;
+
+    auto file_size = std::filesystem::file_size(fname);
+    if (file_size == 0) {
+        throw std::runtime_error(LOCATION +
+                                 "Cannot guess frame size: file is empty");
+    }
+
+    if (file_size % module_data_size == 0) {
+        m_rows = 512;
+        m_cols = 1024;
+        m_bytes_per_frame = module_data_size - header_size;
+    } else if (file_size % half_data_size == 0) {
+        m_rows = 256;
+        m_cols = 1024;
+        m_bytes_per_frame = half_data_size - header_size;
+    } else if (file_size % chip_data_size == 0) {
+        m_rows = 256;
+        m_cols = 256;
+        m_bytes_per_frame = chip_data_size - header_size;
+    } else {
+        throw std::runtime_error(LOCATION +
+                                 "Cannot find frame size: file size is not a "
+                                 "multiple of any known frame size");
+    }
+}
+
+void JungfrauDataFile::parse_fname(const std::filesystem::path &fname) {
+    m_path = fname.parent_path();
+    m_base_name = fname.stem();
+
+    // find file index, then remove if from the base name
+    if (auto pos = m_base_name.find_last_of('_'); pos != std::string::npos) {
+        m_offset = std::stoul(m_base_name.substr(pos + 1));
+        m_base_name.erase(pos);
+    }
+}
+
+void JungfrauDataFile::scan_files() {
+    // find how many files we have and the number of frames in each file
+    m_last_frame_in_file.clear();
+    size_t file_index = m_offset;
+    while (std::filesystem::exists(fpath(file_index))) {
+        auto n_frames = std::filesystem::file_size(fpath(file_index)) /
+                        (m_bytes_per_frame + header_size);
+        m_last_frame_in_file.push_back(n_frames);
+        ++file_index;
+    }
+
+    // find where we need to open the next file and total number of frames
+    m_last_frame_in_file = cumsum(m_last_frame_in_file);
+    m_total_frames = m_last_frame_in_file.back();
+}
+
+void JungfrauDataFile::read_into(std::byte *image_buf,
+                                 JungfrauDataHeader *header) {
+
+    // read header if not passed nullptr
+    if (header) {
+        if (auto rc = fread(header, sizeof(JungfrauDataHeader), 1, m_fp.get());
+            rc != 1) {
+            throw std::runtime_error(
+                LOCATION +
+                "Could not read header from file:" + m_fp.error_msg());
+        }
+    } else {
+        m_fp.seek(header_size, SEEK_CUR);
+    }
+
+    // read data
+    if (auto rc = fread(image_buf, 1, m_bytes_per_frame, m_fp.get());
+        rc != m_bytes_per_frame) {
+        throw std::runtime_error(LOCATION + "Could not read image from file" +
+                                 m_fp.error_msg());
+    }
+
+    // prepare for next read
+    // if we are at the end of the file, open the next file
+    ++m_current_frame_index;
+    if (m_current_frame_index >= m_last_frame_in_file[m_current_file_index] &&
+        (m_current_frame_index < m_total_frames)) {
+        ++m_current_file_index;
+        open_file(m_current_file_index);
+    }
+}
+
+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{
+        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()));
+    }
+    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());
+        rc != sizeof(header)) {
+        throw std::runtime_error(LOCATION + "Could not read header from file" +
+                                 m_fp.error_msg());
+    }
+    m_fp.seek(-header_size, SEEK_CUR);
+    return header;
+}
+
+void JungfrauDataFile::open_file(size_t file_index) {
+    // fmt::print(stderr, "Opening file: {}\n",
+    // fpath(file_index+m_offset).string());
+    m_fp = FilePtr(fpath(file_index + m_offset), "rb");
+    m_current_file_index = file_index;
+}
+
+std::filesystem::path JungfrauDataFile::fpath(size_t file_index) const {
+    auto fname = fmt::format("{}_{:0{}}.dat", m_base_name, file_index,
+                             n_digits_in_file_index);
+    return m_path / fname;
+}
+
+} // namespace aare

src/JungfrauDataFile.test.cpp

@@ -0,0 +1,114 @@
+#include "aare/JungfrauDataFile.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
+    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.pixels_per_frame() == 524288);
+    REQUIRE(f.bytes_per_pixel() == 2);
+    REQUIRE(f.bitdepth() == 16);
+    REQUIRE(f.base_name() == "AldoJF500k");
+    REQUIRE(f.n_files() == 4);
+    REQUIRE(f.tell() == 0);
+    REQUIRE(f.total_frames() == 24);
+    REQUIRE(f.current_file() == fpath);
+
+    //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());
+        f.read_into(&image, &header);
+        REQUIRE(header.framenum == i + 1);
+        REQUIRE(header.bunchid == (i + 1) * (i + 1));
+        REQUIRE(image.shape(0) == 512);
+        REQUIRE(image.shape(1) == 1024);
+    }
+}
+
+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
+    f.seek(19);
+    REQUIRE(f.tell() == 19);
+    auto h = f.read_header();
+    REQUIRE(h.framenum == 19+1);
+
+    //Reading again does not change the file pointer
+    auto h2 = f.read_header();
+    REQUIRE(h2.framenum == 19+1);
+
+    f.seek(59);
+    REQUIRE(f.tell() == 59);
+    auto h3 = f.read_header();
+    REQUIRE(h3.framenum == 59+1);
+
+    JungfrauDataHeader h4;
+    aare::NDArray<uint16_t> image(f.shape());
+    f.read_into(&image, &h4);
+    REQUIRE(h4.framenum == 59+1);
+
+    //now we should be on the next frame
+    REQUIRE(f.tell() == 60);
+    REQUIRE(f.read_header().framenum == 60+1);
+
+    REQUIRE_THROWS(f.seek(86356)); //out of range
+}
+
+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);
+    REQUIRE(f.tell() == 0);
+
+    //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);
+
+    // 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]"){
+    auto fpath = test_data_path() / "dat" / "AldoJF65k_000000.dat";
+    REQUIRE(std::filesystem::exists(fpath));
+
+    JungfrauDataFile f(fpath);
+
+    aare::NDArray<uint16_t> image({39, 85});
+    JungfrauDataHeader header;
+
+    REQUIRE_THROWS(f.read_into(&image, &header));
+    REQUIRE_THROWS(f.read_into(&image, nullptr));
+    REQUIRE_THROWS(f.read_into(&image));
+
+    REQUIRE(f.tell() == 0);
+
+
+}

src/NDArray.test.cpp

@@ -183,14 +183,14 @@ TEST_CASE("Size and shape matches") {
     int64_t h = 75;
     std::array<int64_t, 2> shape{w, h};
     NDArray<double> a{shape};
-    REQUIRE(a.size() == static_cast<uint64_t>(w * h));
+    REQUIRE(a.size() == w * h);
     REQUIRE(a.shape() == shape);
 }
 
 TEST_CASE("Initial value matches for all elements") {
     double v = 4.35;
     NDArray<double> a{{5, 5}, v};
-    for (uint32_t i = 0; i < a.size(); ++i) {
+    for (int i = 0; i < a.size(); ++i) {
         REQUIRE(a(i) == v);
     }
 }

src/algorithm.test.cpp

@@ -5,7 +5,7 @@
 
 TEST_CASE("Find the closed index in a 1D array", "[algorithm]") {
     aare::NDArray<double, 1> arr({5});
-    for (size_t i = 0; i < arr.size(); i++) {
+    for (ssize_t i = 0; i < arr.size(); i++) {
         arr[i] = i;
     }
     // arr 0, 1, 2, 3, 4
@@ -18,7 +18,7 @@ TEST_CASE("Find the closed index in a 1D array", "[algorithm]") {
 
 TEST_CASE("Passing integers to nearest_index works", "[algorithm]") {
     aare::NDArray<int, 1> arr({5});
-    for (size_t i = 0; i < arr.size(); i++) {
+    for (ssize_t i = 0; i < arr.size(); i++) {
         arr[i] = i;
     }
     // arr 0, 1, 2, 3, 4
@@ -47,9 +47,22 @@ TEST_CASE("nearest index works with std::array", "[algorithm]") {
     REQUIRE(aare::nearest_index(arr, -10.0) == 0);
 }
 
+TEST_CASE("nearest index when there is no different uses the first element",
+          "[algorithm]") {
+    std::vector<int> vec = {5, 5, 5, 5, 5};
+    REQUIRE(aare::nearest_index(vec, 5) == 0);
+}
+
+TEST_CASE("nearest index when there is no different uses the first element "
+          "also when all smaller",
+          "[algorithm]") {
+    std::vector<int> vec = {5, 5, 5, 5, 5};
+    REQUIRE(aare::nearest_index(vec, 10) == 0);
+}
+
 TEST_CASE("last smaller", "[algorithm]") {
     aare::NDArray<double, 1> arr({5});
-    for (size_t i = 0; i < arr.size(); i++) {
+    for (ssize_t i = 0; i < arr.size(); i++) {
         arr[i] = i;
     }
     // arr 0, 1, 2, 3, 4
@@ -61,9 +74,89 @@ TEST_CASE("last smaller", "[algorithm]") {
 
 TEST_CASE("returns last bin strictly smaller", "[algorithm]") {
     aare::NDArray<double, 1> arr({5});
-    for (size_t i = 0; i < arr.size(); i++) {
+    for (ssize_t i = 0; i < arr.size(); i++) {
         arr[i] = i;
     }
     // arr 0, 1, 2, 3, 4
-    REQUIRE(aare::last_smaller(arr, 2.0) == 2);
-}
+    REQUIRE(aare::last_smaller(arr, 2.0) == 1);
+}
+
+TEST_CASE("last_smaller with all elements smaller returns last element",
+          "[algorithm]") {
+    aare::NDArray<double, 1> arr({5});
+    for (ssize_t i = 0; i < arr.size(); i++) {
+        arr[i] = i;
+    }
+    // arr 0, 1, 2, 3, 4
+    REQUIRE(aare::last_smaller(arr, 50.) == 4);
+}
+
+TEST_CASE("last_smaller with all elements bigger returns first element",
+          "[algorithm]") {
+    aare::NDArray<double, 1> arr({5});
+    for (ssize_t i = 0; i < arr.size(); i++) {
+        arr[i] = i;
+    }
+    // arr 0, 1, 2, 3, 4
+    REQUIRE(aare::last_smaller(arr, -50.) == 0);
+}
+
+TEST_CASE("last smaller with all elements equal returns the first element",
+          "[algorithm]") {
+    std::vector<int> vec = {5, 5, 5, 5, 5, 5, 5};
+    REQUIRE(aare::last_smaller(vec, 5) == 0);
+}
+
+TEST_CASE("first_lager with vector", "[algorithm]") {
+    std::vector<double> vec = {0, 1, 2, 3, 4};
+    REQUIRE(aare::first_larger(vec, 2.5) == 3);
+}
+
+TEST_CASE("first_lager with all elements smaller returns last element",
+          "[algorithm]") {
+    std::vector<double> vec = {0, 1, 2, 3, 4};
+    REQUIRE(aare::first_larger(vec, 50.) == 4);
+}
+
+TEST_CASE("first_lager with all elements bigger returns first element",
+          "[algorithm]") {
+    std::vector<double> vec = {0, 1, 2, 3, 4};
+    REQUIRE(aare::first_larger(vec, -50.) == 0);
+}
+
+TEST_CASE("first_lager with all elements the same as the check returns last",
+          "[algorithm]") {
+    std::vector<int> vec = {14, 14, 14, 14, 14};
+    REQUIRE(aare::first_larger(vec, 14) == 4);
+}
+
+TEST_CASE("first larger with the same element", "[algorithm]") {
+    std::vector<int> vec = {7, 8, 9, 10, 11};
+    REQUIRE(aare::first_larger(vec, 9) == 3);
+}
+
+TEST_CASE("cumsum works", "[algorithm]") {
+    std::vector<double> vec = {0, 1, 2, 3, 4};
+    auto result = aare::cumsum(vec);
+    REQUIRE(result.size() == vec.size());
+    REQUIRE(result[0] == 0);
+    REQUIRE(result[1] == 1);
+    REQUIRE(result[2] == 3);
+    REQUIRE(result[3] == 6);
+    REQUIRE(result[4] == 10);
+}
+TEST_CASE("cumsum works with empty vector", "[algorithm]") {
+    std::vector<double> vec = {};
+    auto result = aare::cumsum(vec);
+    REQUIRE(result.size() == 0);
+}
+TEST_CASE("cumsum works with negative numbers", "[algorithm]") {
+    std::vector<double> vec = {0, -1, -2, -3, -4};
+    auto result = aare::cumsum(vec);
+    REQUIRE(result.size() == vec.size());
+    REQUIRE(result[0] == 0);
+    REQUIRE(result[1] == -1);
+    REQUIRE(result[2] == -3);
+    REQUIRE(result[3] == -6);
+    REQUIRE(result[4] == -10);
+}