Merge branch 'developer' into jf_h5reader

- avoid default shallow copy of pointer members triggered by copy of derived class
- as a result avoid double-free bugs

slsDetectorCalibration/analogDetector.h

@@ -3,14 +3,14 @@
 #ifndef ANALOGDETECTOR_H
 #define ANALOGDETECTOR_H
 
-//#include <mutex>
+#include <mutex>
 
 #include "commonModeSubtractionNew.h"
 #include "ghostSummation.h"
 #include "pedestalSubtraction.h"
+#include "sls/tiffIO.h"
 #include "slsDetectorData.h"
 #include "slsInterpolation.h"
-#include "sls/tiffIO.h"
 #include <pthread.h>
 
 #ifdef ROOTSPECTRUM
@@ -90,6 +90,8 @@ template <class dataType> class analogDetector {
         ymin = 0;
         ymax = ny;
         fMode = ePedestal;
+        dMode = eInterpolating;
+        // std::cout << "dMode " << dMode << std::endl;
         thr = 0;
         myFile = NULL;
 #ifdef ROOTSPECTRUM
@@ -111,15 +113,20 @@ template <class dataType> class analogDetector {
        destructor. Deletes the pdestalSubtraction array and the image
     */
     virtual ~analogDetector() {
+        // std::cout << "#### Debug: Destructing analogDetector! ####"
+        //           << std::endl;
         for (int i = 0; i < ny; i++) {
-            delete[] stat[i];
+            if (stat[i]) { delete[] stat[i]; stat[i] = nullptr; }
+            // delete[] stat[i];
             /* delete [] pedMean[i];  */
             /* delete [] pedVariance[i]; */
-        };
+        }
         /* delete [] pedMean;  */
         /* delete [] pedVariance; */
-        delete[] stat;
-        delete[] image;
+        // delete[] stat;
+        // delete[] image;
+        if (stat) { delete[] stat; stat = nullptr; }
+        if (image) { delete[] image; image = nullptr; }
 #ifdef ROOTSPECTRUM
         delete hs;
 #ifdef ROOTCLUST
@@ -137,6 +144,8 @@ template <class dataType> class analogDetector {
      */
     analogDetector(analogDetector *orig) {
         /* copy construction from orig*/
+        // std::cout << "#### Debug: Calling analogDetector cloning method! ####"
+        //           << std::endl;
         det = orig->det;
         nx = orig->nx;
         ny = orig->ny;
@@ -152,6 +161,8 @@ template <class dataType> class analogDetector {
         thr = orig->thr;
         // nSigma=orig->nSigma;
         fMode = orig->fMode;
+        dMode = orig->dMode;
+        // std::cout << "dMode " << dMode << std::endl;
         myFile = orig->myFile;
 
         stat = new pedestalSubtraction *[ny];
@@ -216,12 +227,75 @@ template <class dataType> class analogDetector {
             ghSum = NULL;
     }
 
+    /**
+       constructor creating a deep copy of another analog detector
+       \param other analog Detector structure to be copied
+     */
+    analogDetector(const analogDetector &other)
+        : det(other.det), nx(other.nx), ny(other.ny), dataSign(other.dataSign),
+          iframe(other.iframe), gmap(other.gmap), id(other.id),
+          xmin(other.xmin), xmax(other.xmax), ymin(other.ymin),
+          ymax(other.ymax), thr(other.thr), fMode(other.fMode),
+          dMode(other.dMode), myFile(NULL) {
+
+        // std::cout << "#### Debug: Calling analogDetector copy constructor! ####"
+        //           << std::endl;
+
+        // Deep copy the stat array
+        stat = new pedestalSubtraction *[ny];
+        for (int i = 0; i < ny; i++) {
+            stat[i] = new pedestalSubtraction[nx];
+            std::copy(other.stat[i], other.stat[i] + nx, stat[i]);
+        }
+
+        // Deep copy image array
+        image = new int[nx * ny];
+        std::copy(other.image, other.image + (nx * ny), image);
+
+        // Copy common-mode subtraction object (if it exists)
+        if (other.cmSub) {
+            cmSub = other.cmSub->Clone();
+            std::cout << "Copying cmSub" << std::endl;
+        } else {
+            cmSub = nullptr;
+        }
+
+        // Copy ghost summation object (if it exists)
+        if (other.ghSum) {
+            ghSum = other.ghSum->Clone();
+            std::cout << "Copying ghSum" << std::endl;
+        } else {
+            ghSum = nullptr;
+        }
+
+        // Ensure pedestal values are copied properly
+        int nped = other.GetNPedestals(0, 0);
+        for (int iy = 0; iy < ny; ++iy) {
+            for (int ix = 0; ix < nx; ++ix) {
+                stat[iy][ix].SetNPedestals(nped);
+                setPedestal(ix, iy, other.getPedestal(ix, iy),
+                            other.getPedestalRMS(ix, iy),
+                            other.GetNPedestals(ix, iy));
+            }
+        }
+    }
+
     /**
        clone. Must be virtual!
        \returns a clone of the original analog detector
      */
     virtual analogDetector *Clone() = 0;
 
+    /**
+       Deep copy. Must be virtual!
+       This is a new addition because of multithreaded storage cell data (where
+       each sc has its own mutex). If the pure virtual function exists here,
+       EVERY derived class has to overwrite it! That means a Copy() function
+       must also be implemented in any derived class. \returns a deep copy of
+       the original analog detector
+     */
+    virtual analogDetector *Copy() = 0;
+
     /**
        Gives an id to the structure. For debugging purposes in case of
        multithreading. \param i is to be set \returns current id
@@ -548,6 +622,11 @@ template <class dataType> class analogDetector {
         return ped;
     };
 
+    // Const version for use in the copy constructor
+    virtual double getPedestal(int ix, int iy, int cm = 0) const {
+        return stat[iy][ix].getPedestal();
+    };
+
     /**
        gets  pedestal rms (i.e. noise)
        \param ix pixel x coordinate
@@ -566,6 +645,11 @@ template <class dataType> class analogDetector {
         return ped;
     };
 
+    // Const version for use in the copy constructor
+    virtual double getPedestalRMS(int ix, int iy) const {
+        return stat[iy][ix].getPedestalRMS();
+    };
+
     /**
        sets  pedestal
        \param ix pixel x coordinate
@@ -1226,7 +1310,7 @@ template <class dataType> class analogDetector {
     /** gets number of samples for moving average pedestal calculation
            \returns actual number of samples
        */
-    int GetNPedestals(int ix, int iy) {
+    int GetNPedestals(int ix, int iy) const {
         if (ix >= 0 && ix < nx && iy >= 0 && iy < ny)
             return stat[iy][ix].GetNPedestals();
         else

slsDetectorCalibration/circularFifo.h

@@ -50,6 +50,8 @@ template <typename Element> class CircularFifo {
     mutable sem_t data_mutex;
     mutable sem_t free_mutex;
     unsigned int increment(unsigned int idx_) const;
+    int id_;
+    int thread_id_;
 };
 
 template <typename Element> int CircularFifo<Element>::getDataValue() const {
@@ -74,14 +76,18 @@ template <typename Element> int CircularFifo<Element>::getFreeValue() const {
 template <typename Element>
 bool CircularFifo<Element>::push(Element *&item_, bool no_block) {
     // check for fifo full
-    if (no_block && isFull())
-        return false;
+    if (no_block && isFull()) {
+        //std::cout << "Full Fifo at push. Returning." << std::endl;
+        return false; // No space, return immediately
+    }
 
-    sem_wait(&free_mutex);
-    array[tail] = item_;
-    tail = increment(tail);
-    sem_post(&data_mutex);
-    return true;
+    //std::cout << "Thread " << thread_id_ <<" Push Fifo " << id_ << " item " << static_cast<void*>(item_) << std::endl;
+
+    sem_wait(&free_mutex); // Wait for space
+    array[tail] = item_; // Add item to the buffer
+    tail = increment(tail); // Move the tail pointer
+    sem_post(&data_mutex); // Signal that there is new data
+    return true; // Success
 }
 
 /** Consumer only: Removes and returns item from the queue
@@ -94,14 +100,18 @@ bool CircularFifo<Element>::push(Element *&item_, bool no_block) {
 template <typename Element>
 bool CircularFifo<Element>::pop(Element *&item_, bool no_block) {
     // check for fifo empty
-    if (no_block && isEmpty())
-        return false;
+    if (no_block && isEmpty()) {
+        //std::cout << "Empty Fifo at pop. Returning." << std::endl;
+        return false; // No data in fifo, return immediately
+    }
 
-    sem_wait(&data_mutex);
-    item_ = array[head];
-    head = increment(head);
-    sem_post(&free_mutex);
-    return true;
+    //std::cout << "Thread " << thread_id_ << " Pop Fifo " << id_ << " item " << static_cast<void*>(item_) << std::endl;
+
+    sem_wait(&data_mutex); // Wait for data
+    item_ = array[head]; // Retreive item from the current head of the buffer
+    head = increment(head); // Move the head pointer (to point to the next item)
+    sem_post(&free_mutex); // Signal that there is new free space available 
+    return true; //Success
 }
 
 /** Useful for testinng and Consumer check of status

slsDetectorCalibration/dataStructures/HDF5File.cpp

@@ -0,0 +1,482 @@
+#include "HDF5File.h"
+
+#include "ansi.h"
+
+#include <algorithm>
+#include <fmt/ranges.h>
+
+/*
+ * No class member helper functions
+ */
+std::string vectorToString(std::vector<hsize_t> const& v) {
+    return fmt::format("({})", fmt::join(v, ", "));
+}
+
+/*
+ * increment frame offset (if s dimension exists, loop through all before incrementing z) 
+ * should also work if file_dims[1] is not s but x (in that case we ignore it)
+ */
+void conditionalIncrement(std::vector<hsize_t>& vec, hsize_t max_value) {
+
+    if (vec.size() < 3) {
+        throw std::invalid_argument("Vector must have at least 3 elements.");
+    }
+
+    // If vector has 4 elements, increment vec[1] first
+    if (vec.size() == 4) {
+        if (++vec[1] >= max_value) { //max_value is never reached!
+            vec[1] = 0;  // Reset and increment vec[0]
+            ++vec[0];
+        }
+    }
+    // If vector has 3 elements, increment vec[0] directly
+    else if (vec.size() == 3) {
+        ++vec[0];
+    }
+}
+
+void printDatatypeSize(hid_t dataset) {
+
+	hid_t datatype = H5Dget_type(dataset);
+	H5T_class_t class_id = H5Tget_class(datatype);
+	size_t type_size = H5Tget_size(datatype);
+
+	H5Tclose(datatype);
+
+	std::cout << " dataset type class: " << class_id 
+          	  << ", size: " << type_size << " bytes\n";
+
+	//Ensure the read datatype matches a system native type correctly
+	//hid_t read_type = (type_size == 8) ? H5T_NATIVE_LLONG : H5T_NATIVE_UINT;
+
+	//return read_type;
+
+}
+
+
+/* **********************
+ * Class member functions
+ * **********************
+ */
+
+// Default constructor
+/*
+HDF5File::HDF5File () {
+  //InitializeParameters(); //old
+}
+*/
+
+/*
+HDF5File::~HDF5File () {
+	
+	if(current_image)
+		delete [] current_image;
+	
+}
+*/
+
+void HDF5File::SetImageDataPath (std::string const& name) { 
+	std::cout << "Image dataset path set to " << name << std::endl;
+	data_datasetname = name; 
+}
+
+void HDF5File::SetFrameIndexPath (std::string const& name) { 
+	std::cout << "Frame index dataset path set to " << name << std::endl;
+	index_datasetname = name; 
+}
+
+void HDF5File::InitializeDimensions () {
+
+	rank = H5Sget_simple_extent_ndims(dataspace);
+    file_dims.resize(rank);
+    H5Sget_simple_extent_dims(dataspace, file_dims.data(), nullptr);
+
+    std::cout << "Dataset dimensions: " << vectorToString(file_dims) << "\n";
+
+}
+
+std::vector<hsize_t> HDF5File::GetDatasetDimensions() {
+	return file_dims;
+}
+
+std::vector<hsize_t> HDF5File::GetChunkDimensions() {
+	return chunk_dims;
+}
+
+hsize_t HDF5File::GetRank() {
+	return rank;
+}
+
+bool HDF5File::ValidateDimensions () {
+
+	// validate rank
+	if(rank != RANK) {
+		cprintf(RED,"rank found %llu. Expected %d\n", rank, RANK);
+		std::cerr << "Error: Rank could not be validated\n";
+		return false;
+	}
+
+	// validate file dimensions of x and y (assuming those are the last two dimensions of the dataset)
+	if ( (file_dims[file_dims.size()-2] != DEFAULT_X_DIMS) || (file_dims[file_dims.size()-1] != DEFAULT_Y_DIMS) ) {
+		cprintf(RED,"file dimensions of x found %llu. Expected %d\n", file_dims[file_dims.size()-2], DEFAULT_X_DIMS);
+		cprintf(RED,"file dimensions of y found %llu. Expected %d\n", file_dims[file_dims.size()-1], DEFAULT_Y_DIMS);
+		std::cerr << "Error: Dataset dimensions could not be validated\n";
+		return false;
+	}
+
+	cprintf(GREEN, "File rank & dimensions validated.");
+	return true;	
+}
+
+bool HDF5File::ReadChunkDimensions () {
+
+	// Get layout
+	hid_t plist_id = H5Dget_create_plist(dataset);
+
+	if (H5Pget_layout (plist_id) != H5D_CHUNKED)  {
+		cprintf(RED,"NOTE: Dataset is not chunked!\n");
+		std::cerr << "Error: Dataset is not chunked\n";
+		return false;
+	}
+
+	// Get Chunk Dimensions
+	int rank_chunk = H5Pget_chunk (plist_id, 0, nullptr);
+	chunk_dims.resize(rank_chunk);
+	H5Pget_chunk (plist_id, rank_chunk, chunk_dims.data());
+
+	std::cout << "Chunk dimensions: " << vectorToString(chunk_dims) << "\n";
+
+	H5Pclose (plist_id);
+
+	return true;
+
+}
+
+bool HDF5File::ValidateChunkDimensions () {
+
+	// validate rank
+	if(chunk_dims.size() != rank) {
+		cprintf(RED,"Chunk rank does not match dataset rank! Found %lu. Expected %llu\n", chunk_dims.size(), rank);
+		std::cerr << "Error: Chunk rank does not match dataset rank\n";
+		return false;
+	}
+
+	// validate chunk dimensions of x and y (assuming those are the last two dimensions of the dataset)
+	if ( (chunk_dims[chunk_dims.size()-2] != DEFAULT_CHUNK_X_DIMS) || (chunk_dims[chunk_dims.size()-1] != DEFAULT_CHUNK_Y_DIMS) ) {
+		cprintf(RED,"file dimensions of x found %llu. Expected %d\n", chunk_dims[chunk_dims.size()-2], DEFAULT_CHUNK_X_DIMS);
+		cprintf(RED,"file dimensions of y found %llu. Expected %d\n", chunk_dims[chunk_dims.size()-1], DEFAULT_CHUNK_Y_DIMS);
+		std::cerr << "Error: Chunk dimensions could not be validated\n";
+		return false;
+	}
+
+	cprintf(GREEN, "Chunk rank & dimensions validated.");
+	return true;
+
+}
+
+bool HDF5File::OpenFrameIndexDataset() {
+
+	// Get all the frame numbers
+	// Open frame index dataset
+	hid_t fi_dataset = H5Dopen2 (file, index_datasetname.c_str(), H5P_DEFAULT);
+	if (fi_dataset < 0){
+		cprintf (RED,"Could not open frame index dataset %s\n", index_datasetname.c_str());
+		std::cerr << "Error: Could not open frame index dataset\n";
+		return false;
+	}
+
+	hid_t fi_dataspace = H5Dget_space (fi_dataset);
+	int fi_rank = H5Sget_simple_extent_ndims(fi_dataspace);
+	std::vector<hsize_t> fi_dims(fi_rank);
+	H5Sget_simple_extent_dims (fi_dataspace, fi_dims.data(), nullptr);
+
+	std::cout << "Frame index dataset dimensions: " << vectorToString(fi_dims) << "\n";
+
+	// validate size
+	if (fi_dims[0] != file_dims[0]) {
+		cprintf (RED,"Frame index dimensions of z found %llu. Expected %llu\n", fi_dims[0], file_dims[0]);
+		std::cerr << "Error: Z dimension of frame index dataset does not align with z dimension of image dataset\n";
+		H5Sclose (fi_dataspace);
+		H5Dclose (fi_dataset);
+		return false;
+	}
+
+	// allocate frame index memory
+	frame_index_list.resize(fi_dims[0]); //file_dims
+
+	// read and print datatype size of dataset
+	std::cout << "Frame index";
+	printDatatypeSize(fi_dataset);
+
+	// make sure we only read the first column of the frame index dataset (not all storage cells)
+	//NOTE: For XFEL datasets, this may mean that some frame numbers are skipped
+	//(because they assign a unique frame number to every storage cell)
+	//Possibly, there is a cleaner fix for this...
+	std::vector<hsize_t> start(fi_rank,0);
+	std::vector<hsize_t> count(fi_rank,1);
+	count[0] = fi_dims[0];
+	hid_t memspace = H5Screate_simple (fi_rank, count.data(), nullptr);
+	if (memspace < 0) {
+    	std::cerr << "Error: Failed to create memory space for HDF5 read operation\n";
+		H5Sclose(memspace);
+		H5Sclose(fi_dataspace);
+		H5Dclose(fi_dataset);
+    	return false;
+	}
+
+	// Create hyperslab selection
+	if (H5Sselect_hyperslab(fi_dataspace, H5S_SELECT_SET, start.data(), nullptr, count.data(), nullptr) < 0 ) {
+		cprintf (RED,"Could not create hyperslab for %s\n", vectorToString(start).c_str());
+		std::cerr << "Error: Hyperslab creation failed for " << vectorToString(start) << "\n";
+		H5Sclose(memspace);
+		H5Sclose (fi_dataspace);
+		H5Dclose(fi_dataset);
+		return false;
+	}
+
+	//read frame index values
+	//Is u32 correct? I would think not. But I get a segmentation fault if I use u64.
+	if (H5Dread (fi_dataset, H5T_STD_U64LE, memspace, fi_dataspace, H5P_DEFAULT, frame_index_list.data()) < 0) {
+		cprintf (RED,"Could not read frame index dataset %s\n", index_datasetname.c_str());
+		std::cerr << "Error: Could not read frame index dataset\n";
+		H5Sclose(memspace);
+		H5Sclose (fi_dataspace);
+		H5Dclose (fi_dataset);
+		return false;
+	}
+
+	H5Sclose(memspace);
+	H5Sclose (fi_dataspace);
+	H5Dclose(fi_dataset);
+	return true;
+}
+
+int HDF5File::OpenResources (char const*const fname, bool validate) {
+
+	std::cout << "Debug HDF5File.cpp: Attempting to open file " << fname << std::endl;
+	// Open File
+	file = H5Fopen (fname, H5F_ACC_RDONLY, H5P_DEFAULT);
+	if (file < 0) {
+	  cprintf(RED,"Could not open hdf5 file\n");
+	  std::cerr << "Error: H5Fopen failed\n";
+	  return 0;
+	}
+	cprintf(BLUE, "Opened File: %s\n", fname);
+
+	// Open Dataset
+	dataset = H5Dopen2 (file, data_datasetname.c_str(), H5P_DEFAULT);
+	if (dataset < 0){
+	  cprintf(RED,"Could not open dataset\n");
+	  std::cerr << "Error: H5Dopen2 failed\n";
+	  CloseResources ();
+	  return 0;
+	}
+	cprintf(BLUE, "Opened Dataset: %s\n", data_datasetname.c_str());
+
+	// print datatype size of dataset
+	std::cout << "Image";
+	printDatatypeSize(dataset);
+
+	// Create Dataspace
+	dataspace = H5Dget_space (dataset);
+	if (dataspace < 0){
+	  cprintf(RED,"Could not open dataspace\n");
+	  std::cerr << "Error: H5Dget_space failed\n";
+	  CloseResources ();
+	  return 0;
+	}
+
+	// Get Dimensions
+	InitializeDimensions();
+	// Get chunk dimensions
+	if (!ReadChunkDimensions()) {
+		CloseResources();
+		return 0;
+	}
+
+	// validate file dimensions
+	if (validate) {
+		if ( !ValidateDimensions() || !ValidateChunkDimensions() ) {
+			CloseResources();
+			return 0;		
+		}
+	}
+
+	//Read frame indices
+	if (!OpenFrameIndexDataset()) {
+		CloseResources();
+		return 0;
+	}
+	
+	return 1;
+}
+
+
+void HDF5File::CloseResources () {
+	if (dataspace >=0 ) 	{
+		H5Sclose(dataspace);
+		dataspace = -1;
+	}
+	if (dataset >=0 ) 	{
+		H5Dclose(dataset);
+		dataset = -1;
+	}
+	if (file >=0 ) 		{
+		H5Fclose(file);
+		file = -1;
+	}
+}
+
+/*
+ * Function takes uint16_t* argument to make explicit that the caller has to handle memory allocation and deallocation.
+ * This is legacy caused by the structure with which the slsDetectorCalibration cluster finder is written.
+ * (Best practice for modern C++ would be using smart pointers.)
+ *
+ * Originially, this function took uint16_t** which may lead to memory management issues since image gets redirected
+ * to point to current_image, which is owned by HDF5File.
+ * (Good practice in classic C-style. HDF5File needs to clean up the resource at destruction.)
+ * 
+ * \param image pointer to uint16_t, buffer which the image is read into. (Memory handled by caller!)
+ * \param offset contains iFrame at [0] and storage cell number at [1],
+ * depending on dimensionality of the dataset, the storage cell number may not be included.
+ * Note that frame number (as read from file) may (likely) differ from frame index (in the dataset)!
+ */
+int HDF5File::ReadImage (uint16_t* image, std::vector<hsize_t>& offset ) {
+
+	// Validate input arguments
+    if (!image) {
+        std::cerr << "Error: image buffer is null.\n";
+        return -99;
+    }
+  	if ( offset.size() != rank-2 ) {
+		cprintf ( RED,"Offset vector must have size %llu. Found %lu\n", rank-2, offset.size() );
+        std::cerr << "Error: Wrong offset vector size\n";
+		CloseResources ();
+        return -99;
+    }
+
+	// Initialize frame_offset
+	if (frame_offset.empty())
+		frame_offset.resize(rank,0);
+
+	// Check if we reached the end of file
+	// Compares that the offsets of frame and storage cell (Z and S) have reached the end of file
+	// Excludes X and Y indices (of the image dataset) from the comparison
+	// As it is now, this never triggers, because frame_offset[1] is never equals file_dims[1]=16
+	/*
+	if( std::equal( frame_offset.cbegin(), frame_offset.cend()-2, file_dims.cbegin() ) ) {
+		printf("End of file reached\n");
+		return -1;
+	}
+	*/
+	if (frame_offset[0] == file_dims[0]) {
+		printf("End of file reached\n");
+		return -1;
+	}
+	/* //old
+	if (frame_offset[0] == file_dims[0]-1) {
+		printf("end of file\n");
+		return -1;
+	}
+	*/
+
+	// Validate frame_offset index
+    if (frame_offset[0] >= frame_index_list.size()) {
+        std::cerr << "Error: frame_offset[0] = " << frame_offset[0] << " of bounds.\n";
+        return -99;
+    }
+
+	// Check if images exist at the current frame offset
+	if (frame_index_list[frame_offset[0]] == 0) {
+		cprintf (RED,"No images at this frame offset %llu\n", frame_offset[0]);
+		std::cerr << "Error: Framenumber 0 at this frame offset\n";
+		CloseResources ();
+		return -99;
+	}
+
+	// Optional: Ensure dataset and dataspace are valid
+    if (dataset < 0) {
+        std::cerr << "Error: Invalid dataset ID.\n";
+        return -99;
+    }
+    if (dataspace < 0) {
+        std::cerr << "Error: Invalid dataspace.\n";
+        return -99;
+    }
+
+	// Define the size of the hyperslab to read
+	std::vector<hsize_t> frame_size(rank, 1);
+	std::copy(file_dims.begin() + rank-2, file_dims.end(), frame_size.begin() + rank-2);
+	/*
+	for ( int d=0; d < rank; ++d ) {
+		if (d < rank-2)
+			frame_size[d] = 1;
+		if ( d >= rank-2 )
+			frame_size[d] = file_dims[d];
+	}
+	*/
+
+	// Define memory space
+	hid_t memspace = H5Screate_simple (rank, frame_size.data(), nullptr);
+	if (memspace < 0) {
+    	std::cerr << "Error: Failed to create memory space for HDF5 read operation\n";
+    	CloseResources();
+    	return -99;
+	}
+
+	// Create hyperslab selection
+	// This aligns dataspace such that we read the correct frame
+	if (H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, frame_offset.data(), nullptr, frame_size.data(), nullptr) < 0 ) {
+		cprintf (RED,"Could not create hyperslab for frame offset %s\n", vectorToString(frame_offset).c_str());
+		std::cerr << "Error: Hyperslab creation failed for frame offset " << vectorToString(frame_offset) << "\n";
+		CloseResources();
+		H5Sclose(memspace);
+		return -99;
+	}
+
+	// Read dataset into image buffer (previously read to current_image owned by HDF5File)
+	if (H5Dread(dataset, H5T_STD_U16LE, memspace, dataspace, H5P_DEFAULT, image) < 0 ) {
+		cprintf (RED,"Could not read dataset for frame offset %s\n", vectorToString(frame_offset).c_str());
+		std::cerr << "Error: Reading of dataset failed for given start frame offset " << vectorToString(frame_offset) << "\n";
+		CloseResources ();
+		H5Sclose(memspace);
+		return -99;
+	}
+
+	// Clean up memory space
+	H5Sclose(memspace);
+
+	//*image = current_image; //if uint16_t** is passed, HDF5File owns the resource image points to, which is potentially dangerous
+
+	// Return frame number
+	unsigned int retval = frame_index_list[frame_offset[0]];
+	
+	// Pass updated frame offset value(s) via offset parameter vector
+	std::copy_n(frame_offset.begin(), offset.size(), offset.begin());
+	/*
+	std::transform( offset.begin(), offset.end(), offset.begin(), 
+    	[&, i = 0](size_t) mutable { return frame_offset[i++]; } );
+	*/
+
+	// Increment frame offset correctly
+	conditionalIncrement(frame_offset, file_dims[1]);
+	//++frame_offset[0]; //old
+
+	return retval;
+}
+
+void HDF5File::PrintCurrentImage (uint16_t* image) {
+	printf("\n");
+	printf("Frame %llu, Image: %llu\n", frame_offset[0]-1, frame_index_list[frame_offset[0]-1]);
+
+	hsize_t size = file_dims[rank-1] * file_dims[rank-2];
+	for (hsize_t i = 0; i < size; ++i){
+		printf("%u ", image[i]);
+		if (!((i+1) % file_dims[rank-2] ))
+			printf("\n\n");
+	}
+	printf("\n\n\n\n");
+}
+
+
+

slsDetectorCalibration/dataStructures/HDF5File.h

@@ -0,0 +1,159 @@
+#pragma once
+/************************************************
+ * @file HDF5Fle.h
+ * @short functions to open/close/read HDF5 File
+ * Adapted for generalization, accepts rank 3 and 4
+ * Supports control over storage cells
+ ***********************************************/
+/**
+ *@short functions to open/close/read HDF5 File
+ */
+
+
+#include <iostream>
+#include <vector>
+#include <string>
+
+#include "hdf5.h"
+#include "hdf5_hl.h"
+
+
+//#define MAX_STR_LENGTH	1000
+
+#define RANK	4                       // Dimension of the image dataset, only for validation
+#define DEFAULT_Z_DIMS	10000			// only for validation
+#define DEFAULT_Y_DIMS	1024			// only for validation
+#define DEFAULT_X_DIMS	512				// only for validation
+//#define DEFAULT_S_DIMS	1       	// Storage cells
+
+#define DEFAULT_CHUNK_Z_DIMS	1		// only for validation
+#define DEFAULT_CHUNK_Y_DIMS	1024	// only for validation
+#define DEFAULT_CHUNK_X_DIMS	512		// only for validation
+//#define DEFAULT_CHUNK_S_DIMS	1
+
+
+#define DATA_DATASETNAME	"/data/JF18T01V01/data" //Furka JF
+#define INDEX_DATASETNAME	"/data/JF18T01V01/frame_index"
+
+//enum{Z,S,X,Y}; //S is the storage cell //enum is not used
+
+class HDF5File {
+
+public:
+
+	/**
+	 * Constructor
+	 */
+	//HDF5File () = default; //No need to declare if it is default
+
+	/**
+	 * Destructor
+	 */
+	//~HDF5File () = default; //Since the destructor is default (and copy and move are default too)
+
+
+	std::vector<hsize_t> GetDatasetDimensions ();
+
+	std::vector<hsize_t> GetChunkDimensions ();
+
+	hsize_t GetRank ();
+
+	void SetImageDataPath (std::string const& name);
+
+	void SetFrameIndexPath (std::string const& name);
+
+
+
+	/**
+	 * Open HDF5 file and dataset,
+	 * reads frame index dataset to array
+	 * @param fname file name
+	 * @param validate true if one must validate if file is
+	 * chunked with dims [? x 128 x 512] and chunk dims [1 x 128 x 512]
+	 * @returns 1 if successful, else 0 if fail
+	 */
+	int OpenResources (const char* const fname, bool validate); 
+
+	/**
+	 * Close Open resources
+	 */
+	void CloseResources ();
+
+	/**
+	 * Read an image into current_image,
+	 * increment Z-offset (frame) and (if rank==4) storage cell
+	 * @returns frame number read from file,
+	 */
+	int ReadImage (uint16_t* image, std::vector<hsize_t>& offset);
+
+	/**
+	 * Print current image in memory
+	 */
+	void PrintCurrentImage (uint16_t* image);
+
+private:
+
+	/**
+	 * Initialize dimensions of image dataset for each new file
+	 */
+	void InitializeDimensions ();
+
+	bool ReadChunkDimensions ();
+
+	bool ValidateDimensions ();
+
+	bool ValidateChunkDimensions ();
+
+	/**
+	 * Open dataset containing the frame numbers
+	 */
+	bool OpenFrameIndexDataset ();
+
+
+	/** file name */
+	std::string file_name{};
+	/** dataset name for image data */
+	std::string data_datasetname = DATA_DATASETNAME;
+	/** dataset name for frame index data */
+	std::string index_datasetname = INDEX_DATASETNAME;
+
+	/** file handle */
+	hid_t file{};
+	/** dataspace handle */
+	hid_t dataspace{};
+	/** memory space handle */
+	//hid_t memspace; //old
+	/** dataset handle */
+	hid_t dataset{};
+
+	/** file dimensions */
+	std::vector<hsize_t> file_dims{};
+	//hsize_t file_dims[RANK]{}; //static array (dimensions are known)
+
+	/** chunk dimensions 
+	 ** not necessarily required
+	 ** useful for optimization or validation */
+	std::vector<hsize_t> chunk_dims{};
+	//hsize_t chunk_dims[RANK]{};
+
+	/** Rank of the image dataset */
+	hsize_t rank{};
+
+	/** number of frames */
+	unsigned int number_of_frames{};
+
+	/** frame index list */
+	std::vector<hsize_t> frame_index_list{};
+
+	/** Current image
+	 ** dynamic array 
+	 ** uint16_t pointer format is chosen to support use with slsDetectorCalibration cluster finder */
+	//uint16_t* current_image{nullptr};
+	//uint16_t current_chunk[DEFAULT_CHUNK_Z_DIMS][DEFAULT_CHUNK_Y_DIMS][DEFAULT_CHUNK_X_DIMS];
+
+	/** Current frame offset
+	 ** (Z-offset, S-offset, 0, 0) or (Z-offset, 0, 0), increments automatically with ReadImage */
+	std::vector<hsize_t> frame_offset{};
+	//hsize_t frame_offset[RANK]{};
+
+};

slsDetectorCalibration/dataStructures/jungfrauLGADStrixelsDataQuad.h

@@ -0,0 +1,422 @@
+// SPDX-License-Identifier: LGPL-3.0-or-other
+// Copyright (C) 2021 Contributors to the SLS Detector Package
+#ifndef JUNGFRAULGADSTRIXELSDATAQUAD_H
+#define JUNGFRAULGADSTRIXELSDATAQUAD_H
+#ifdef CINT
+#include "sls/sls_detector_defs_CINT.h"
+#else
+#include "sls/sls_detector_defs.h"
+#endif
+#include "slsDetectorData.h"
+
+// #define VERSION_V2
+/**
+    @short  structure for a Detector Packet or Image Header
+    @li frameNumber is the frame number
+    @li expLength is the subframe number (32 bit eiger) or real time exposure
+   time in 100ns (others)
+    @li packetNumber is the packet number
+    @li bunchId is the bunch id from beamline
+    @li timestamp is the time stamp with 10 MHz clock
+    @li modId is the unique module id (unique even for left, right, top, bottom)
+    @li xCoord is the x coordinate in the complete detector system
+    @li yCoord is the y coordinate in the complete detector system
+    @li zCoord is the z coordinate in the complete detector system
+    @li debug is for debugging purposes
+    @li roundRNumber is the round robin set number
+    @li detType is the detector type see :: detectorType
+    @li version is the version number of this structure format
+*/
+
+#include <algorithm>
+#include <numeric>
+#include <tuple>
+
+namespace strixelQuad {
+  constexpr int nc_rawimg = 1024; // for full images  //256;
+  constexpr int nc_quad = 512;
+  constexpr int nr_rawimg = 512;
+  constexpr int nr_chip = 256;
+  constexpr int gr = 9;
+
+  //shift due to extra pixels
+  constexpr int shift_x = 2; //left
+  
+  constexpr int nc_strixel = ( nc_quad - shift_x - 2*gr ) / 3; //164
+  constexpr int nr_strixel = ( nr_chip - 1 - gr ) * 3; //one half (-1 because double sided pixel) //738
+  constexpr int nr_center = 12; //double sided pixels to be skipped
+
+  // boundaries in ASIC coordinates (pixels at both bounds are included)
+  constexpr int xstart = 256 + gr;             // 265
+  constexpr int xend = 255 + nc_quad - gr;     // 758
+  constexpr int bottom_ystart = gr;            //   9
+  constexpr int bottom_yend = nr_chip - 2;     // 254
+  constexpr int top_ystart = nr_chip + 1;      // 257
+  constexpr int top_yend = nr_chip*2 - gr - 1; // 502
+
+  // x shift because of 2-pixel strixels on one side
+  constexpr int shift = 2;
+
+} // namespace strixelQuad
+
+//to account for module rotation
+enum rotation {
+  NORMAL = 0,
+  INVERSE = 1
+};
+
+const int rota = NORMAL;
+
+typedef struct {
+    uint64_t bunchNumber; /**< is the frame number */
+    uint64_t pre;         /**< something */
+
+} jf_header; // Aldo's header
+
+using namespace strixelQuad;
+
+class jungfrauLGADStrixelsDataQuad : public slsDetectorData<uint16_t> {
+
+  private:
+    int iframe;
+    int x0, y0, x1, y1, shifty;
+    struct {
+      uint16_t xmin;
+      uint16_t xmax;
+      uint16_t ymin;
+      uint16_t ymax;
+      int nc;
+    } globalROI;
+
+    //to account for the inverted routing of the two different quad halfs
+    enum location {
+      BOTTOM = 0,
+      TOP = 1
+    };
+
+    int multiplicator = 3;
+    std::vector<int> mods{ 0, 1, 2 };
+
+    void reverseVector( std::vector<int>& v ) {
+      std::reverse( v.begin(), v.end() );
+      std::cout << "mods reversed ";
+      for ( auto i : v )
+	std::cout << i << " ";
+      std::cout << '\n';
+    }
+
+    void setMappingShifts( const int rot, const int half ) {
+
+      x0 = xstart;
+      x1 = xend;
+
+      if (rot==NORMAL) {
+	x0 += shift;
+      } else {
+	x1-=shift;
+        reverseVector(mods);
+      }
+
+      if (half==BOTTOM) {
+	y0 = bottom_ystart;
+	y1 = bottom_yend;
+	shifty = 0;
+      } else {
+	y0 = top_ystart;
+	y1 = top_yend;
+	reverseVector(mods);
+	shifty = nr_strixel + nr_center; //double-sided pixels in the center have to be jumped
+      }
+
+    }
+
+    void remap( int xmin=0, int xmax=0, int ymin=0, int ymax=0 ) {
+
+      int ix, iy = 0;
+      // remapping loop
+        for (int ipy = y0; ipy <= y1; ipy++) {
+            for (int ipx = x0; ipx <= x1; ipx++) {
+
+                ix = int((ipx - x0) / multiplicator);
+                for (int m = 0; m < multiplicator; ++m) {
+                    if ((ipx - x0) % multiplicator == m)
+                        iy = (ipy - y0) * multiplicator + mods[m] + shifty;
+                }
+
+                //	  if (iy< 40)	  cout << iy << "  " << ix <<endl;
+		if (xmin < xmax && ymin < ymax) {
+		  if ( ipx>=xmin && ipx<=xmax && ipy>=ymin && ipy <=ymax )
+		    dataMap[iy][ix] =
+		      sizeof(header) + (globalROI.nc * (ipy - globalROI.ymin) + (ipx - globalROI.xmin)) * 2;
+		} else {
+		  dataMap[iy][ix] = sizeof(header) + (nc_rawimg * ipy + ipx) * 2;
+		}
+            }
+        }
+
+    }
+
+    void remapQuad(const int rot) {
+
+        setMappingShifts( rot, BOTTOM );
+	remap();
+	setMappingShifts( rot, TOP );
+	remap();
+
+    }
+
+
+    std::tuple< uint16_t, uint16_t, uint16_t, uint16_t > adjustROItoLimits(uint16_t xmin,
+									   uint16_t xmax,
+									   uint16_t ymin,
+									   uint16_t ymax,
+									   uint16_t lim_roi_xmin,
+									   uint16_t lim_roi_xmax,
+									   uint16_t lim_roi_ymin,
+									   uint16_t lim_roi_ymax) {
+      uint16_t xmin_roi, xmax_roi, ymin_roi, ymax_roi;
+      if ( xmin < lim_roi_xmin)
+        xmin_roi = lim_roi_xmin;
+      else
+        xmin_roi = xmin;
+      if ( xmax > lim_roi_xmax )
+        xmax_roi = lim_roi_xmax;
+      else
+        xmax_roi = xmax;
+      if ( ymin < lim_roi_ymin )
+        ymin_roi = lim_roi_ymin;
+      else
+        ymin_roi = ymin;
+      if ( ymax > lim_roi_ymax )
+        ymax_roi = lim_roi_ymax;
+      else
+        ymax_roi = ymax;
+      return std::make_tuple(xmin_roi, xmax_roi, ymin_roi, ymax_roi);
+    }
+
+    std::vector < std::tuple< int, uint16_t, uint16_t, uint16_t, uint16_t > > mapSubROIs(uint16_t xmin,
+											 uint16_t xmax,
+											 uint16_t ymin,
+											 uint16_t ymax) {
+      bool bottom = false;
+      bool top = false;
+
+      for ( int x=xmin; x!=xmax+1; ++x ) {
+	for ( int y=ymin; y!=ymax; ++y ) {
+	  if ( xstart<=x && x<=xend && bottom_ystart<=y && y<=bottom_yend )
+	    bottom = true;
+	  if ( xstart<=x && x<=xend && top_ystart<=y && y<=top_yend )
+	    top = true;
+	}
+      }
+
+      uint16_t xmin_roi{}, xmax_roi{}, ymin_roi{}, ymax_roi{};
+      std::vector < std::tuple< int, uint16_t, uint16_t, uint16_t, uint16_t > > rois{};
+
+      if (bottom) {
+	std::tie( xmin_roi, xmax_roi, ymin_roi, ymax_roi ) =
+	  adjustROItoLimits( xmin, xmax, ymin, ymax,
+			     xstart, xend, bottom_ystart, bottom_yend );
+	rois.push_back( std::make_tuple( BOTTOM, xmin_roi, xmax_roi, ymin_roi, ymax_roi ) );
+      }
+      if (top) {
+	std::tie( xmin_roi, xmax_roi, ymin_roi, ymax_roi ) =
+	  adjustROItoLimits( xmin, xmax, ymin, ymax,
+			     xstart, xend, top_ystart, top_yend );
+	rois.push_back( std::make_tuple( TOP, xmin_roi, xmax_roi, ymin_roi, ymax_roi ) );
+      }
+
+      return rois;
+    }
+
+    void remapROI(std::tuple< int, uint16_t, uint16_t, uint16_t, uint16_t > roi, const int rot ) {
+
+      int half, xmin, xmax, ymin, ymax;
+      std::tie( half, xmin, xmax, ymin, ymax ) = roi;
+	
+      setMappingShifts(rot, half);
+
+      std::cout << "remapping roi: "
+		<< ", x0: " << x0 << ", x1: " << x1 << ", y0: " << y0
+		<< ", y1: " << y1 << std::endl;
+      std::cout << "Adjusted roi: [" << xmin << ", " << xmax << ", " << ymin << ", " << ymax << "]" << std::endl;
+
+      remap( xmin, xmax, ymin, ymax );
+    }
+
+  public:
+
+    using header = sls::defs::sls_receiver_header;
+
+    jungfrauLGADStrixelsDataQuad(uint16_t xmin = 0, uint16_t xmax = 0,
+				 uint16_t ymin = 0, uint16_t ymax = 0)
+        : slsDetectorData<uint16_t>(
+              nc_strixel,
+              nr_strixel * 2 + nr_center,
+              nc_strixel * ( nr_strixel * 2 + nr_center ) * 2 + sizeof(header)) {
+        std::cout << "Jungfrau strixels quad with full module data "
+                  << std::endl;
+
+        // Fill all strixels with dummy values
+        for (int ix = 0; ix != nc_strixel; ++ix) {
+            for (int iy = 0; iy != nr_strixel * 2 + nr_center; ++iy) {
+                dataMap[iy][ix] = sizeof(header);
+            }
+        }
+
+	globalROI.xmin = xmin;
+	globalROI.xmax = xmax;
+	globalROI.ymin = ymin;
+	globalROI.ymax = ymax;
+
+        std::cout << "sizeofheader = " << sizeof(header) << std::endl;
+        std::cout << "Jungfrau strixels quad with full module data "
+                  << std::endl;
+
+        if (xmin < xmax && ymin < ymax) {
+
+	  // get ROI raw image number of columns
+	  globalROI.nc = xmax - xmin + 1;
+	  std::cout << "nc_roi = " << globalROI.nc << std::endl;
+	  
+	  dataSize =
+	    (xmax - xmin + 1) * (ymax - ymin + 1) * 2 + sizeof(header);
+	  std::cout << "datasize " << dataSize << std::endl;
+
+	  auto rois = mapSubROIs(xmin, xmax, ymin, ymax);
+	  //function to fill vector of rois from globalROI
+
+	  for ( auto roi : rois )
+	    remapROI(roi, rota);
+
+        } else {
+
+	  remapQuad( rota );
+
+        }
+
+        iframe = 0;
+        std::cout << "data struct created" << std::endl;
+    };
+
+    /**
+       Returns the value of the selected channel for the given dataset as
+       double. \param data pointer to the dataset (including headers etc) \param
+       ix pixel number in the x direction \param iy pixel number in the y
+       direction \returns data for the selected channel, with inversion if
+       required as double
+
+    */
+    virtual double getValue(char *data, int ix, int iy = 0) {
+
+        uint16_t val = getChannel(data, ix, iy) & 0x3fff;
+        return val;
+    };
+
+    /**
+
+    Returns the frame number for the given dataset. Purely virtual func.
+    \param buff pointer to the dataset
+    \returns frame number
+
+    */
+
+    int getFrameNumber(char *buff) {
+#ifdef ALDO                                      // VH
+        return ((jf_header *)buff)->bunchNumber; // VH
+#endif                                           // VH
+        return ((header *)buff)->detHeader.frameNumber;
+    };
+
+    /**
+
+       Returns the packet number for the given dataset. purely virtual func
+       \param buff pointer to the dataset
+       \returns packet number number
+
+
+
+    */
+    int getPacketNumber(char *buff) {
+#ifdef ALDO // VH
+        // uint32_t fakePacketNumber = 1000;
+        // return fakePacketNumber; //VH //TODO: Keep in mind in case of bugs!
+        // //This is definitely bad!
+        return 1000;
+#endif // VH
+        return ((header *)buff)->detHeader.packetNumber;
+    };
+
+    char *readNextFrame(std::ifstream &filebin) {
+        int ff = -1, np = -1;
+        return readNextFrame(filebin, ff, np);
+    };
+
+    char *readNextFrame(std::ifstream &filebin, int &ff) {
+        int np = -1;
+        return readNextFrame(filebin, ff, np);
+    };
+
+    char *readNextFrame(std::ifstream &filebin, int &ff, int &np) {
+        char *data = new char[dataSize];
+        char *d = readNextFrame(filebin, ff, np, data);
+        if (d == NULL) {
+            delete[] data;
+            data = NULL;
+        }
+        return data;
+    };
+
+    char *readNextFrame(std::ifstream &filebin, int &ff, int &np, char *data) {
+      //char *retval = 0;
+      //int nd;
+      //int fnum = -1;
+        np = 0;
+        //int pn;
+
+        //std::cout << dataSize << std::endl;
+        //if (ff >= 0) {
+          //  fnum = ff; }
+
+        if (filebin.is_open()) {
+            if (filebin.read(data, dataSize)) {
+                std::cout << "*";
+                ff = getFrameNumber(data);
+                np = getPacketNumber(data);
+                return data;
+            }
+	    std::cout << "#";
+        } else {
+	  std::cout << "File not open" << std::endl;
+	}
+        return NULL;
+    };
+
+    /*    Loops over a memory slot until a complete frame is found (i.e. all */
+    /*    packets 0 to nPackets, same frame number). purely virtual func \param
+     */
+    /*    data pointer to the memory to be analyzed \param ndata reference to
+     * the */
+    /*    amount of data found for the frame, in case the frame is incomplete at
+     */
+    /*    the end of the memory slot \param dsize size of the memory slot to be
+     */
+    /*    analyzed \returns pointer to the beginning of the last good frame
+     * (might */
+    /*    be incomplete if ndata smaller than dataSize), or NULL if no frame is
+     */
+    /*    found */
+
+    /* *\/ */
+    virtual char *findNextFrame(char *data, int &ndata, int dsize) {
+        if (dsize < dataSize)
+            ndata = dsize;
+        else
+            ndata = dataSize;
+        return data;
+    };
+
+    // int getPacketNumber(int x, int y) {return dataMap[y][x]/packetSize;};
+};
+
+#endif

slsDetectorCalibration/dataStructures/jungfrauLGADStrixelsDataQuadH5.h

@@ -0,0 +1,432 @@
+// SPDX-License-Identifier: LGPL-3.0-or-other
+// Copyright (C) 2021 Contributors to the SLS Detector Package
+#ifndef JUNGFRAULGADSTRIXELSDATAQUADH5_H
+#define JUNGFRAULGADSTRIXELSDATAQUADH5_H
+#ifdef CINT
+#include "sls/sls_detector_defs_CINT.h"
+#else
+#include "sls/sls_detector_defs.h"
+#endif
+#include "slsDetectorData.h"
+
+// This needs to be linked correctly
+#include "HDF5File.cpp"
+#include "HDF5File.h" //this includes hdf5.h and hdf5_hl.h
+
+// #define VERSION_V2
+/**
+    @short  structure for a Detector Packet or Image Header
+    @li frameNumber is the frame number
+    @li expLength is the subframe number (32 bit eiger) or real time exposure
+   time in 100ns (others)
+    @li packetNumber is the packet number
+    @li bunchId is the bunch id from beamline
+    @li timestamp is the time stamp with 10 MHz clock
+    @li modId is the unique module id (unique even for left, right, top, bottom)
+    @li xCoord is the x coordinate in the complete detector system
+    @li yCoord is the y coordinate in the complete detector system
+    @li zCoord is the z coordinate in the complete detector system
+    @li debug is for debugging purposes
+    @li roundRNumber is the round robin set number
+    @li detType is the detector type see :: detectorType
+    @li version is the version number of this structure format
+*/
+
+// #include <algorithm>
+#include <numeric>
+#include <tuple>
+
+namespace strixelQuad {
+constexpr int nc_rawimg = 1024; // for full images  //256;
+constexpr int nc_quad = 512;
+constexpr int nr_rawimg = 512;
+constexpr int nr_chip = 256;
+constexpr int gr = 9;
+
+// shift due to extra pixels
+constexpr int shift_x = 2; // left
+
+constexpr int nc_strixel = (nc_quad - shift_x - 2 * gr) / 3; // 164
+constexpr int nr_strixel =
+    (nr_chip - 1 - gr) * 3;   // one half (-1 because double sided pixel) //738
+constexpr int nr_center = 12; // double sided pixels to be skipped
+
+// boundaries in ASIC coordinates (pixels at both bounds are included)
+constexpr int xstart = 256 + gr;               // 265
+constexpr int xend = 255 + nc_quad - gr;       // 758
+constexpr int bottom_ystart = gr;              //   9
+constexpr int bottom_yend = nr_chip - 2;       // 254
+constexpr int top_ystart = nr_chip + 1;        // 257
+constexpr int top_yend = nr_chip * 2 - gr - 1; // 502
+
+// x shift because of 2-pixel strixels on one side
+constexpr int shift = 2;
+
+} // namespace strixelQuad
+
+// to account for module rotation
+enum rotation { NORMAL = 0, INVERSE = 1 };
+
+const int rota = NORMAL;
+
+typedef struct {
+    uint64_t bunchNumber; /**< is the frame number */
+    uint64_t pre;         /**< something */
+
+} jf_header; // Aldo's header
+
+using namespace strixelQuad;
+
+class jungfrauLGADStrixelsDataQuadH5 : public slsDetectorData<uint16_t> {
+
+  private:
+    int iframe;
+    int x0, y0, x1, y1, shifty;
+    struct {
+        uint16_t xmin;
+        uint16_t xmax;
+        uint16_t ymin;
+        uint16_t ymax;
+        int nc;
+    } globalROI;
+
+    // to account for the inverted routing of the two different quad halfs
+    enum location { BOTTOM = 0, TOP = 1 };
+
+    int multiplicator = 3;
+    std::vector<int> mods{0, 1, 2};
+
+    void reverseVector(std::vector<int> &v) {
+        std::reverse(v.begin(), v.end());
+        std::cout << "mods reversed ";
+        for (auto i : v)
+            std::cout << i << " ";
+        std::cout << '\n';
+    }
+
+    void setMappingShifts(const int rot, const int half) {
+
+        x0 = xstart;
+        x1 = xend;
+
+        if (rot == NORMAL) {
+            x0 += shift;
+        } else {
+            x1 -= shift;
+            reverseVector(mods);
+        }
+
+        if (half == BOTTOM) {
+            y0 = bottom_ystart;
+            y1 = bottom_yend;
+            shifty = 0;
+        } else {
+            y0 = top_ystart;
+            y1 = top_yend;
+            reverseVector(mods);
+            shifty = nr_strixel + nr_center; // double-sided pixels in the
+                                             // center have to be jumped
+        }
+    }
+
+    void remap(int xmin = 0, int xmax = 0, int ymin = 0, int ymax = 0) {
+
+        int ix, iy = 0;
+        // remapping loop
+        for (int ipy = y0; ipy <= y1; ++ipy) {
+            for (int ipx = x0; ipx <= x1; ++ipx) {
+
+                ix = int((ipx - x0) / multiplicator);
+                for (int m = 0; m < multiplicator; ++m) {
+                    if ((ipx - x0) % multiplicator == m)
+                        iy = (ipy - y0) * multiplicator + mods[m] + shifty;
+                }
+
+                //	  if (iy< 40)	  cout << iy << "  " << ix <<endl;
+                if (xmin < xmax && ymin < ymax) { // if ROI
+                    if (ipx >= xmin && ipx <= xmax && ipy >= ymin &&
+                        ipy <= ymax)
+                        dataMap[iy][ix] =
+                            (globalROI.nc * (ipy - globalROI.ymin) +
+                             (ipx - globalROI.xmin)) *
+                            2;
+                } else { // if full Quad
+                    dataMap[iy][ix] = (nc_rawimg * ipy + ipx) * 2;
+                }
+            }
+        }
+    }
+
+    void remapQuad(const int rot) {
+
+        setMappingShifts(rot, BOTTOM);
+        remap();
+        setMappingShifts(rot, TOP);
+        remap();
+    }
+
+    std::tuple<uint16_t, uint16_t, uint16_t, uint16_t>
+    adjustROItoLimits(uint16_t xmin, uint16_t xmax, uint16_t ymin,
+                      uint16_t ymax, uint16_t lim_roi_xmin,
+                      uint16_t lim_roi_xmax, uint16_t lim_roi_ymin,
+                      uint16_t lim_roi_ymax) {
+        uint16_t xmin_roi, xmax_roi, ymin_roi, ymax_roi;
+        if (xmin < lim_roi_xmin)
+            xmin_roi = lim_roi_xmin;
+        else
+            xmin_roi = xmin;
+        if (xmax > lim_roi_xmax)
+            xmax_roi = lim_roi_xmax;
+        else
+            xmax_roi = xmax;
+        if (ymin < lim_roi_ymin)
+            ymin_roi = lim_roi_ymin;
+        else
+            ymin_roi = ymin;
+        if (ymax > lim_roi_ymax)
+            ymax_roi = lim_roi_ymax;
+        else
+            ymax_roi = ymax;
+        return std::make_tuple(xmin_roi, xmax_roi, ymin_roi, ymax_roi);
+    }
+
+    // The strixel Quad has a mirrored symmetry from the center axis
+    // So we need to distinguish between bottom and top half for remapping
+    std::vector<std::tuple<int, uint16_t, uint16_t, uint16_t, uint16_t>>
+    mapSubROIs(uint16_t xmin, uint16_t xmax, uint16_t ymin, uint16_t ymax) {
+        bool bottom = false;
+        bool top = false;
+
+        for (int x = xmin; x != xmax + 1; ++x) {
+            for (int y = ymin; y != ymax; ++y) {
+                if (xstart <= x && x <= xend && bottom_ystart <= y &&
+                    y <= bottom_yend)
+                    bottom = true;
+                if (xstart <= x && x <= xend && top_ystart <= y &&
+                    y <= top_yend)
+                    top = true;
+            }
+        }
+
+        uint16_t xmin_roi{}, xmax_roi{}, ymin_roi{}, ymax_roi{};
+        std::vector<std::tuple<int, uint16_t, uint16_t, uint16_t, uint16_t>>
+            rois{};
+
+        if (bottom) {
+            std::tie(xmin_roi, xmax_roi, ymin_roi, ymax_roi) =
+                adjustROItoLimits(xmin, xmax, ymin, ymax, xstart, xend,
+                                  bottom_ystart, bottom_yend);
+            rois.push_back(std::make_tuple(BOTTOM, xmin_roi, xmax_roi, ymin_roi,
+                                           ymax_roi));
+        }
+        if (top) {
+            std::tie(xmin_roi, xmax_roi, ymin_roi, ymax_roi) =
+                adjustROItoLimits(xmin, xmax, ymin, ymax, xstart, xend,
+                                  top_ystart, top_yend);
+            rois.push_back(
+                std::make_tuple(TOP, xmin_roi, xmax_roi, ymin_roi, ymax_roi));
+        }
+
+        return rois;
+    }
+
+    void remapROI(std::tuple<int, uint16_t, uint16_t, uint16_t, uint16_t> roi,
+                  const int rot) {
+
+        int half, xmin, xmax, ymin, ymax;
+        std::tie(half, xmin, xmax, ymin, ymax) = roi;
+
+        setMappingShifts(rot, half);
+
+        std::cout << "remapping roi: "
+                  << ", x0: " << x0 << ", x1: " << x1 << ", y0: " << y0
+                  << ", y1: " << y1 << std::endl;
+        std::cout << "Adjusted roi: [" << xmin << ", " << xmax << ", " << ymin
+                  << ", " << ymax << "]" << std::endl;
+
+        remap(xmin, xmax, ymin, ymax);
+    }
+
+    // The following functions are pure virtual in the base class. But I don't
+    // want them to be accessible here! Implement the functions as private (to
+    // satisfy the linker) int getFrameNumber(char* buff){return 0;} //This is
+    // actually needed because the cluster finder writes the framenumber
+    int getPacketNumber(char *buff) { return 0; } // Not provided
+
+    // Mark overwritten functions as override final
+    char *readNextFrame(std::ifstream &filebin) override final {
+        return nullptr;
+    }
+
+  public:
+    using header = sls::defs::sls_receiver_header;
+
+    jungfrauLGADStrixelsDataQuadH5(uint16_t xmin = 0, uint16_t xmax = 0,
+                                   uint16_t ymin = 0, uint16_t ymax = 0)
+        : slsDetectorData<uint16_t>(
+              // nc_strixel,
+              // nr_strixel * 2 + nr_center,
+              // nc_strixel * ( nr_strixel * 2 + nr_center ) * 2
+              512 / 2, 1024 * 2, 512 * 1024 * 2) {
+        std::cout << "Jungfrau strixels quad with full module data "
+                  << std::endl;
+
+        // Fill all strixels with dummy values
+        // for (int ix = 0; ix != nc_strixel; ++ix) {
+        //    for (int iy = 0; iy != nr_strixel * 2 + nr_center; ++iy) {
+        for (int ix = 0; ix != 512 / 2; ++ix) {
+            for (int iy = 0; iy != 1024 * 2; ++iy) {
+                // Set everything to dummy value
+                dataMap[iy][ix] = sizeof(header);
+            }
+        }
+
+        globalROI.xmin = xmin;
+        globalROI.xmax = xmax;
+        globalROI.ymin = ymin;
+        globalROI.ymax = ymax;
+
+        // std::cout << "sizeofheader = " << sizeof(header) << std::endl;
+        std::cout << "Jungfrau strixels quad with full module data "
+                  << std::endl;
+
+        if (xmin < xmax && ymin < ymax) {
+
+            // get ROI raw image number of columns
+            globalROI.nc = xmax - xmin + 1;
+            std::cout << "nc_roi = " << globalROI.nc << std::endl;
+
+            dataSize = (xmax - xmin + 1) * (ymax - ymin + 1) * 2;
+            std::cout << "datasize " << dataSize << std::endl;
+
+            auto rois = mapSubROIs(xmin, xmax, ymin, ymax);
+            // function to fill vector of rois from globalROI
+
+            for (auto roi : rois)
+                remapROI(roi, rota);
+
+        } else {
+            remapQuad(rota);
+        }
+
+        iframe = 0;
+        std::cout << "data struct created" << std::endl;
+    };
+
+    /**
+       Returns the value of the selected channel for the given dataset as
+       double. \param data pointer to the dataset (including headers etc) \param
+       ix pixel number in the x direction \param iy pixel number in the y
+       direction \returns data for the selected channel, with inversion if
+       required as double
+
+    */
+    virtual double getValue(char *data, int ix, int iy = 0) {
+
+        uint16_t val = getChannel(data, ix, iy) & 0x3fff;
+        return val;
+    };
+
+    char *readNextFrame(HDF5File &hfile) {
+        int fn = 0;
+        std::vector<hsize_t> h5offset(1);
+        return readNextFrame(hfile, fn, h5offset);
+    };
+
+    char *readNextFrame(HDF5File &hfile, int &fn) {
+        std::vector<hsize_t> h5offset(1);
+        return readNextFrame(hfile, fn, h5offset);
+    };
+
+    char *readNextFrame(HDF5File &hfile, int &fn,
+                        std::vector<hsize_t> &h5offset) {
+
+        // Ensure dataSize is a valid size for allocation
+        if (dataSize <= 0) {
+            // Handle error case appropriately, e.g., log an error message
+            return nullptr;
+        }
+
+        char *data = new char[dataSize];
+        char *readResult = readNextFrame(hfile, fn, h5offset, data);
+
+        // Check if reading failed
+        if (readResult == nullptr) {
+            delete[] data;  // Free allocated memory
+            data = nullptr; // Set to nullptr to avoid dangling pointer
+        }
+
+        return data; // returning data is equivalent to returning
+                     // reinterpret_cast<char*>(data_ptr) as they both point to
+                     // the same memory
+    };
+
+    /*
+     * This is the most recent function. This is used in the cluster finder!
+     * The overloads are legacy!
+     * Note that caller has to allocate and deallocate memory for data!
+     * \param hfile object of type HDF5File (reader class)
+     * \param framenumber frame number as read from the HDF5 file
+     * \param h5offset vector defining offset parameters for HDF5 hyperslab
+     * selection (dimensions Z and S), incremented automatially
+     * \param data pointer to image buffer (converted to hold uint16_t by
+     * definition of HDF5File)
+     */
+    char *readNextFrame(HDF5File &hfile, int &framenumber,
+                        std::vector<hsize_t> &h5offset, char *data) {
+
+        if (framenumber >= 0) {
+            if (h5offset[0] % 10 == 0)
+                std::cout << "*";
+
+            // Storing the reinterpret_cast in the variable data_ptr ensures
+            // that I can pass it to a function that expects at uint16_t*
+            uint16_t *data_ptr = reinterpret_cast<uint16_t *>(
+                data); // now data_ptr points where data points (thus modifies
+                       // the same memory)
+
+            framenumber = hfile.ReadImage(data_ptr, h5offset);
+            iframe = h5offset[0]; // iframe is a class member!
+            return data; // return reinterpret_cast<char*>(data_ptr); //
+                         // Equivalent
+        }
+
+        std::cout << "#";
+        return nullptr;
+    };
+
+    int getFrameNumber(char *buff) {
+        return iframe;
+    } // Provided via public method readNextFrame
+    // It is debatable if one might not instead want to provide the "real" frame
+    // number as read from the file here For now, this is the frame offset
+    // counter (that always has to start at 0 for each new file)
+
+    /*    Loops over a memory slot until a complete frame is found (i.e. all */
+    /*    packets 0 to nPackets, same frame number). purely virtual func \param
+     */
+    /*    data pointer to the memory to be analyzed \param ndata reference to
+     * the */
+    /*    amount of data found for the frame, in case the frame is incomplete at
+     */
+    /*    the end of the memory slot \param dsize size of the memory slot to be
+     */
+    /*    analyzed \returns pointer to the beginning of the last good frame
+     * (might */
+    /*    be incomplete if ndata smaller than dataSize), or NULL if no frame is
+     */
+    /*    found */
+
+    /* *\/ */
+    virtual char *findNextFrame(char *data, int &ndata, int dsize) {
+        if (dsize < dataSize)
+            ndata = dsize;
+        else
+            ndata = dataSize;
+        return data;
+    };
+
+    // int getPacketNumber(int x, int y) {return dataMap[y][x]/packetSize;};
+};
+
+#endif

slsDetectorCalibration/dataStructures/jungfrauLGADStrixelsData_new.h

@@ -28,6 +28,10 @@
     @li version is the version number of this structure format
 */
 
+#include <algorithm>
+#include <numeric>
+#include <tuple>
+
 namespace strixelSingleChip {
 constexpr int nc_rawimg = 1024; // for full images  //256;
 constexpr int nr_rawimg = 512;
@@ -77,7 +81,7 @@ constexpr int c6g1_ystart = c6g2_yend + 1;        // 448
 constexpr int c6g1_yend = c6g2_yend + 64 - gr;    // 502
 
 // y shift due to faulty bonding (relevant for M408)
-constexpr int bond_shift_y = 1; // CHANGE IF YOU CHANGE MODULE!
+constexpr int bond_shift_y = 0; // CHANGE IF YOU CHANGE MODULE!
 
 } // namespace strixelSingleChip
 
@@ -97,6 +101,13 @@ class jungfrauLGADStrixelsData : public slsDetectorData<uint16_t> {
     int chip_x0;
     int chip_y0;
     int x0, y0, x1, y1, shifty;
+    struct {
+      uint16_t xmin;
+      uint16_t xmax;
+      uint16_t ymin;
+      uint16_t ymax;
+      int nc;
+    } globalROI;
 
     int getMultiplicator(const int group) {
         int multiplicator;
@@ -215,9 +226,113 @@ class jungfrauLGADStrixelsData : public slsDetectorData<uint16_t> {
         }
     }
 
-    void remapROI(uint16_t xmin, uint16_t xmax, uint16_t ymin, uint16_t ymax) {
+
+    std::tuple< uint16_t, uint16_t, uint16_t, uint16_t > adjustROItoLimits(uint16_t xmin,
+									   uint16_t xmax,
+									   uint16_t ymin,
+									   uint16_t ymax,
+									   uint16_t lim_roi_xmin,
+									   uint16_t lim_roi_xmax,
+									   uint16_t lim_roi_ymin,
+									   uint16_t lim_roi_ymax) {
+      uint16_t xmin_roi, xmax_roi, ymin_roi, ymax_roi;
+      if ( xmin < lim_roi_xmin)
+        xmin_roi = lim_roi_xmin;
+      else
+        xmin_roi = xmin;
+      if ( xmax > lim_roi_xmax )
+        xmax_roi = lim_roi_xmax;
+      else
+        xmax_roi = xmax;
+      if ( ymin < lim_roi_ymin )
+        ymin_roi = lim_roi_ymin;
+      else
+        ymin_roi = ymin;
+      if ( ymax > lim_roi_ymax )
+        ymax_roi = lim_roi_ymax;
+      else
+        ymax_roi = ymax;
+      return std::make_tuple(xmin_roi, xmax_roi, ymin_roi, ymax_roi);
+    }
+
+    std::vector < std::tuple< int, int, uint16_t, uint16_t, uint16_t, uint16_t > > mapSubROIs(uint16_t xmin,
+											      uint16_t xmax,
+											      uint16_t ymin,
+											      uint16_t ymax) {
+      bool chip_1_1 = false;
+      bool chip_1_2 = false;
+      bool chip_1_3 = false;
+      bool chip_6_1 = false;
+      bool chip_6_2 = false;
+      bool chip_6_3 = false;
+
+      for ( int x=xmin; x!=xmax+1; ++x ) {
+	for ( int y=ymin; y!=ymax; ++y ) {
+	  if ( c1g1_xstart<=x && x<=c1_xend && (c1g1_ystart+bond_shift_y)<=y && y<=(c1g1_yend+bond_shift_y) )
+	    chip_1_1 = true;
+	  if ( c1g2_xstart<=x && x<=c1_xend && (c1g2_ystart+bond_shift_y)<=y && y<=(c1g2_yend+bond_shift_y) )
+	    chip_1_2 = true;
+	  if ( c1g3_xstart<=x && x<=c1_xend && (c1g3_ystart+bond_shift_y)<=y && y<=(c1g3_yend+bond_shift_y) )
+	    chip_1_3 = true;
+	  if ( c6_xstart<=x && x<=c6g1_xend && (c6g1_ystart-bond_shift_y)<=y && y<=(c6g1_yend-bond_shift_y) )
+	    chip_6_1 = true;
+	  if ( c6_xstart<=x && x<=c6g2_xend && (c6g2_ystart-bond_shift_y)<=y && y<=(c6g2_yend-bond_shift_y) )
+	    chip_6_2 = true;
+	  if ( c6_xstart<=x && x<=c6g3_xend && (c6g3_ystart-bond_shift_y)<=y && y<=(c6g3_yend-bond_shift_y) )
+	    chip_6_3 = true;
+	}
+      }
+
+      uint16_t xmin_roi{}, xmax_roi{}, ymin_roi{}, ymax_roi{};
+       //[ chip, group, xmin, xmax, ymin, ymax ]
+      std::vector < std::tuple< int, int, uint16_t, uint16_t, uint16_t, uint16_t > > rois{};
+
+      if (chip_1_1) {
+	std::tie( xmin_roi, xmax_roi, ymin_roi, ymax_roi ) =
+	  adjustROItoLimits( xmin, xmax, ymin, ymax,
+			     c1g1_xstart, c1_xend, 0, c1g1_yend+bond_shift_y );
+	rois.push_back( std::make_tuple( 1, 1, xmin_roi, xmax_roi, ymin_roi, ymax_roi ) );
+      }
+      if (chip_1_2) {
+	std::tie( xmin_roi, xmax_roi, ymin_roi, ymax_roi ) =
+	  adjustROItoLimits( xmin, xmax, ymin, ymax,
+			     c1g2_xstart, c1_xend, c1g2_ystart+bond_shift_y, c1g2_yend+bond_shift_y );
+	rois.push_back( std::make_tuple( 1, 2, xmin_roi, xmax_roi, ymin_roi, ymax_roi ) );
+      }
+      if (chip_1_3) {
+	std::tie( xmin_roi, xmax_roi, ymin_roi, ymax_roi ) =
+	  adjustROItoLimits( xmin, xmax, ymin, ymax,
+			     c1g3_xstart, c1_xend, c1g3_ystart+bond_shift_y, c1g3_yend+bond_shift_y );
+	rois.push_back( std::make_tuple( 1, 3, xmin_roi, xmax_roi, ymin_roi, ymax_roi ) );
+      }
+      if (chip_6_3) {
+	std::tie( xmin_roi, xmax_roi, ymin_roi, ymax_roi ) =
+	  adjustROItoLimits( xmin, xmax, ymin, ymax,
+			     c6_xstart, c6g3_xend, c6g3_ystart-bond_shift_y, c6g3_yend-bond_shift_y );
+	rois.push_back( std::make_tuple( 6, 3, xmin_roi, xmax_roi, ymin_roi, ymax_roi ) );
+      }
+      if (chip_6_2) {
+	std::tie( xmin_roi, xmax_roi, ymin_roi, ymax_roi ) =
+	  adjustROItoLimits( xmin, xmax, ymin, ymax,
+			     c6_xstart, c6g2_xend, c6g2_ystart-bond_shift_y, c6g2_yend-bond_shift_y );
+	rois.push_back( std::make_tuple( 6, 2, xmin_roi, xmax_roi, ymin_roi, ymax_roi ) );
+      }
+      if (chip_6_1) {
+	std::tie( xmin_roi, xmax_roi, ymin_roi, ymax_roi ) =
+	  adjustROItoLimits( xmin, xmax, ymin, ymax,
+			     c6_xstart, c6g1_xend, c6g1_ystart-bond_shift_y, 511 );
+	rois.push_back( std::make_tuple( 6, 1, xmin_roi, xmax_roi, ymin_roi, ymax_roi ) );
+      }
+
+      return rois;
+
+    }
+
+    void remapROI(std::tuple< int, int, uint16_t, uint16_t, uint16_t, uint16_t > roi) {
         // determine group and chip selected by ROI
-        int group;
+      int group, xmin, xmax, ymin, ymax;
+      std::tie( mchip, group, xmin, xmax, ymin, ymax ) = roi;
+	/*
         if (ymax <= c1g1_yend + bond_shift_y) {
             group = 1;
             mchip = 1;
@@ -243,18 +358,17 @@ class jungfrauLGADStrixelsData : public slsDetectorData<uint16_t> {
 	  group = -1;
 	  mchip = -1;
 	}
-        int multiplicator = getMultiplicator(group);
-        setMappingShifts(group);
-
-        std::cout << "chip: " << mchip << ", group: " << group << ", m: " << multiplicator
-                  << ", x0: " << x0 << ", x1: " << x1 << ", y0: " << y0
-                  << ", y1: " << y1 << std::endl;
-
-        // get ROI raw image number of columns
-        int nc_roi = xmax - xmin + 1;
-        std::cout << "nc_roi = " << nc_roi << std::endl;
+	*/
+      int multiplicator = getMultiplicator(group);
+      setMappingShifts(group);
+      
+      std::cout << "remapping chip: " << mchip << ", group: " << group << ", m: " << multiplicator
+		<< ", x0: " << x0 << ", x1: " << x1 << ", y0: " << y0
+		<< ", y1: " << y1 << std::endl;
+      std::cout << "Adjusted roi: [" << xmin << ", " << xmax << ", " << ymin << ", " << ymax << "]" << std::endl;
 
         // make sure loop bounds are correct
+	/*
         if (y0 < ymin)
             std::cout << "Error ymin" << std::endl;
         if (y1 > ymax)
@@ -264,24 +378,27 @@ class jungfrauLGADStrixelsData : public slsDetectorData<uint16_t> {
             std::cout << "Error xmin" << std::endl;
         if (x1 > xmax)
             std::cout << "Error xmax" << std::endl;
+	*/
 
         // remapping loop
-        int ix, iy = 0;
-        for (int ipy = y0; ipy <= y1; ++ipy) {
-            for (int ipx = x0; ipx <= x1; ++ipx) {
+      int ix, iy = 0;
+      for (int ipy = y0; ipy <= y1; ++ipy) {
+	for (int ipx = x0; ipx <= x1; ++ipx) {
 
-                ix = int((ipx - x0 /*-xmin*/) / multiplicator);
-                for (int m = 0; m < multiplicator; m++) {
-                    if ((ipx - x0 /*-xmin*/) % multiplicator == m)
-                        iy = (ipy - y0 /*-ymin*/) * multiplicator + m + shifty;
-                }
-
-                //	  if (iy< 40)	  cout << iy << "  " << ix <<endl;
-                dataMap[iy][ix] =
-                    sizeof(header) + (nc_roi * (ipy - ymin) + (ipx - xmin)) * 2;
-                groupmap[iy][ix] = group - 1;
-            }
-        }
+	  ix = int((ipx - x0 /*-xmin*/) / multiplicator);
+	  for (int m = 0; m < multiplicator; m++) {
+	    if ((ipx - x0 /*-xmin*/) % multiplicator == m)
+	      iy = (ipy - y0 /*-ymin*/) * multiplicator + m + shifty;
+	  }
+	  
+	  //	  if (iy< 40)	  cout << iy << "  " << ix <<endl;
+	  if ( ipx>=xmin && ipx<=xmax && ipy>=ymin && ipy <=ymax )
+	    dataMap[iy][ix] =
+	      sizeof(header) + (globalROI.nc * (ipy - globalROI.ymin) + (ipx - globalROI.xmin)) * 2;
+	  else dataMap[iy][ix] = sizeof(header);
+	  groupmap[iy][ix] = group - 1;
+	}
+      }
     }
 
   public:
@@ -307,16 +424,31 @@ class jungfrauLGADStrixelsData : public slsDetectorData<uint16_t> {
             }
         }
 
+	globalROI.xmin = xmin;
+	globalROI.xmax = xmax;
+	globalROI.ymin = ymin;
+	globalROI.ymax = ymax;
+
         std::cout << "sizeofheader = " << sizeof(header) << std::endl;
         std::cout << "Jungfrau strixels 2X single chip with full module data "
                   << std::endl;
 
         if (xmin < xmax && ymin < ymax) {
 
-            dataSize =
-                (xmax - xmin + 1) * (ymax - ymin + 1) * 2 + sizeof(header);
-            std::cout << "datasize " << dataSize << std::endl;
-            remapROI(xmin, xmax, ymin, ymax);
+	  // get ROI raw image number of columns
+	  globalROI.nc = xmax - xmin + 1;
+	  std::cout << "nc_roi = " << globalROI.nc << std::endl;
+	  
+	  dataSize =
+	    (xmax - xmin + 1) * (ymax - ymin + 1) * 2 + sizeof(header);
+	  std::cout << "datasize " << dataSize << std::endl;
+
+	  //[ chip, group, xmin, xmax, ymin, ymax ]
+	  auto rois = mapSubROIs(xmin, xmax, ymin, ymax);
+	  //function to fill vector of rois from globalROI
+
+	  for ( auto roi : rois )
+	    remapROI(roi);
 
         } else {
 

slsDetectorCalibration/dataStructures/jungfrauModuleData.h

@@ -3,6 +3,7 @@
 #ifndef JUNGFRAUMODULEDATA_H
 #define JUNGFRAUMODULEDATA_H
 #include <cstdint>
+#include "sls/sls_detector_defs.h"
 #include "slsDetectorData.h"
 
 //#define VERSION_V2
@@ -27,7 +28,7 @@ typedef struct {
     uint64_t bunchNumber; /**< is the frame number */
     uint64_t pre;         /**< something */
 
-} jf_header;
+} jf_header; //Aldo's header!
 
 using namespace std;  
 class jungfrauModuleData : public slsDetectorData<uint16_t> {
@@ -42,20 +43,56 @@ class jungfrauModuleData : public slsDetectorData<uint16_t> {
        1286 large etc.) \param c crosstalk parameter for the output buffer
 
    */
+
+#ifdef ALDO
+    using header = jf_header;
+#else
+    using header = sls::defs::sls_receiver_header;
+#endif
+
 #ifndef ZMQ
-#define off sizeof(jf_header)
+#define off sizeof(header)
 #endif
 #ifdef ZMQ
 #define off 0
 #endif
 
 
-    jungfrauModuleData()
+ jungfrauModuleData(uint16_t xmin=0, uint16_t xmax=0,
+		    uint16_t ymin=0, uint16_t ymax=0)
         : slsDetectorData<uint16_t>(1024, 512,
                                     1024* 512 * 2 + off) {
-  
-        for (int ix = 0; ix < 1024; ix++) {
-            for (int iy = 0; iy < 512; iy++) {
+
+      for (int ix = 0; ix != 1024; ++ix) {
+	for (int iy = 0; iy != 512; ++iy) {
+	  dataMap[iy][ix] = off;
+	}
+      }
+
+      if (xmin < xmax && ymin < ymax) {
+
+	int nc_roi = xmax - xmin + 1;
+	int nr_roi = ymax - ymin + 1;
+	std::cout << "nc_roi = " << nc_roi << std::endl;
+	std::cout << "nr_roi = " << nr_roi << std::endl;
+	  
+	dataSize =
+	  (xmax - xmin + 1) * (ymax - ymin + 1) * 2 + off;
+	std::cout << "datasize " << dataSize << std::endl;
+
+	for (int ix = xmin; ix < xmax+1; ++ix) {
+            for (int iy = ymin; iy < ymax+1; ++iy) {
+                dataMap[iy][ix] = off + (nc_roi * iy + ix) * 2;
+#ifdef HIGHZ
+                dataMask[iy][ix] = 0x3fff;
+#endif
+            }
+        }
+
+      } else {
+	  
+        for (int ix = 0; ix < 1024; ++ix) {
+            for (int iy = 0; iy < 512; ++iy) {
                 dataMap[iy][ix] = off + (1024 * iy + ix) * 2;
 #ifdef HIGHZ
                 dataMask[iy][ix] = 0x3fff;
@@ -63,7 +100,7 @@ class jungfrauModuleData : public slsDetectorData<uint16_t> {
             }
         }
 
-
+      }
 
 
         iframe = 0;
@@ -150,7 +187,7 @@ class jungfrauModuleData : public slsDetectorData<uint16_t> {
         //int pn;
 
         //  cout << dataSize << endl;
-        if (ff >= 0)
+        //if (ff >= 0)
 	  //fnum = ff;
 
         if (filebin.is_open()) {

slsDetectorCalibration/dataStructures/jungfrauSingleChipData.h

@@ -0,0 +1,226 @@
+// SPDX-License-Identifier: LGPL-3.0-or-other
+// Copyright (C) 2021 Contributors to the SLS Detector Package
+#ifndef JUNGFRAUSINGLECHIPDATA_H
+#define JUNGFRAUSINGLECHIPDATA_H
+#include <cstdint>
+#include "sls/sls_detector_defs.h"
+#include "slsDetectorData.h"
+
+//#define VERSION_V2
+/**
+    @short  structure for a Detector Packet or Image Header
+    @li frameNumber is the frame number
+    @li expLength is the subframe number (32 bit eiger) or real time exposure
+   time in 100ns (others)
+    @li packetNumber is the packet number
+    @li bunchId is the bunch id from beamline
+    @li timestamp is the time stamp with 10 MHz clock
+    @li modId is the unique module id (unique even for left, right, top, bottom)
+    @li xCoord is the x coordinate in the complete detector system
+    @li yCoord is the y coordinate in the complete detector system
+    @li zCoord is the z coordinate in the complete detector system
+    @li debug is for debugging purposes
+    @li roundRNumber is the round robin set number
+    @li detType is the detector type see :: detectorType
+    @li version is the version number of this structure format
+*/
+typedef struct {
+    uint64_t bunchNumber; /**< is the frame number */
+    uint64_t pre;         /**< something */
+
+} jf_header; //Aldo's header!
+
+using namespace std;  
+class jungfrauSingleChipData : public slsDetectorData<uint16_t> {
+
+  private:
+    int iframe;
+
+  public:
+    /**
+       Implements the slsReceiverData structure for the moench02 prototype read
+       out by a module i.e. using the slsReceiver (160x160 pixels, 40 packets
+       1286 large etc.) \param c crosstalk parameter for the output buffer
+
+   */
+
+#ifdef ALDO
+    using header = jf_header;
+#else
+    using header = sls::defs::sls_receiver_header;
+#endif
+
+#ifndef ZMQ
+#define off sizeof(header)
+#endif
+#ifdef ZMQ
+#define off 0
+#endif
+
+
+ jungfrauSingleChipData(uint16_t xmin=0, uint16_t xmax=0,
+		    uint16_t ymin=0, uint16_t ymax=0)
+        : slsDetectorData<uint16_t>(256, 256,
+                                    256* 256 * 2 + off) {
+
+      for (int ix = 0; ix != 256; ++ix) {
+	for (int iy = 0; iy != 256; ++iy) {
+	  dataMap[iy][ix] = off;
+	}
+      }
+
+      if (xmin < xmax && ymin < ymax) {
+
+	int nc_roi = xmax - xmin + 1;
+	int nr_roi = ymax - ymin + 1;
+	std::cout << "nc_roi = " << nc_roi << std::endl;
+	std::cout << "nr_roi = " << nr_roi << std::endl;
+	  
+	dataSize =
+	  (xmax - xmin + 1) * (ymax - ymin + 1) * 2 + off;
+	std::cout << "datasize " << dataSize << std::endl;
+
+	for (int ix = xmin; ix < xmax+1; ++ix) {
+            for (int iy = ymin; iy < ymax+1; ++iy) {
+                dataMap[iy][ix] = off + (nc_roi * iy + ix) * 2;
+#ifdef HIGHZ
+                dataMask[iy][ix] = 0x3fff;
+#endif
+            }
+        }
+
+      } else {
+	  
+        for (int ix = 0; ix < 256; ++ix) {
+            for (int iy = 0; iy < 256; ++iy) {
+                dataMap[iy][ix] = off + (256 * iy + ix) * 2;
+#ifdef HIGHZ
+                dataMask[iy][ix] = 0x3fff;
+#endif
+            }
+        }
+
+      }
+
+
+        iframe = 0;
+        //  cout << "data struct created" << endl;
+    };
+
+    /**
+       Returns the value of the selected channel for the given dataset as
+       double. \param data pointer to the dataset (including headers etc) \param
+       ix pixel number in the x direction \param iy pixel number in the y
+       direction \returns data for the selected channel, with inversion if
+       required as double
+
+    */
+    virtual double getValue(char *data, int ix, int iy = 0) {
+
+        uint16_t val = getChannel(data, ix, iy) & 0x3fff;
+	/* if (ix==0 && iy==0) */
+	/*   cout << val << endl; */
+        return val;
+    };
+
+   
+
+    /**
+
+    Returns the frame number for the given dataset. Purely virtual func.
+    \param buff pointer to the dataset
+    \returns frame number
+
+ */
+
+    /* class jfrau_packet_header_t { */
+    /*  public: */
+    /* 	unsigned char reserved[4]; */
+    /* 	unsigned char packetNumber[1]; */
+    /* 	unsigned char frameNumber[3]; */
+    /* 	unsigned char bunchid[8]; */
+    /* }; */
+
+    int getFrameNumber(char *buff) {
+        return ((jf_header *)buff)->bunchNumber;
+    };
+  
+    /**
+
+       Returns the packet number for the given dataset. purely virtual func
+       \param buff pointer to the dataset
+       \returns packet number number
+
+
+
+    */
+    int getPacketNumber(char *buff) {
+        return 0;
+    };
+
+   
+    char *readNextFrame(ifstream &filebin) {
+        int ff = -1, np = -1;
+        return readNextFrame(filebin, ff, np);
+    };
+
+    char *readNextFrame(ifstream &filebin, int &ff) {
+        int np = -1;
+        return readNextFrame(filebin, ff, np);
+    };
+
+    char *readNextFrame(ifstream &filebin, int &ff, int &np) {
+        char *data = new char[dataSize];
+        char *d = readNextFrame(filebin, ff, np, data);
+        if (d == NULL) {
+            delete[] data;
+            data = NULL;
+        }
+        return data;
+    };
+
+    char *readNextFrame(ifstream &filebin, int &ff, int &np,char *data) {
+      //char *retval = 0;
+      //int nd;
+      //int fnum = -1;
+        np = 0;
+        //int pn;
+
+        //  cout << dataSize << endl;
+        //if (ff >= 0)
+	  //fnum = ff;
+
+        if (filebin.is_open()) {
+            if (filebin.read(data, dataSize)) {
+                ff = getFrameNumber(data);
+                np = getPacketNumber(data);
+                return data;
+            }
+        }
+        return NULL;
+    };
+
+    /**
+
+       Loops over a memory slot until a complete frame is found (i.e. all
+       packets 0 to nPackets, same frame number). purely virtual func \param
+       data pointer to the memory to be analyzed \param ndata reference to the
+       amount of data found for the frame, in case the frame is incomplete at
+       the end of the memory slot \param dsize size of the memory slot to be
+       analyzed \returns pointer to the beginning of the last good frame (might
+       be incomplete if ndata smaller than dataSize), or NULL if no frame is
+       found
+
+    */
+    virtual char *findNextFrame(char *data, int &ndata, int dsize) {
+        if (dsize < dataSize)
+            ndata = dsize;
+        else
+            ndata = dataSize;
+        return data;
+    };
+
+    // int getPacketNumber(int x, int y) {return dataMap[y][x]/packetSize;};
+};
+
+#endif

slsDetectorCalibration/interpolatingDetector.h

@@ -54,6 +54,10 @@ class interpolatingDetector : public singlePhotonDetector {
         pthread_mutex_init(fi, NULL);
     };
 
+    /**
+       pointer-based copy constructor (cloner)
+       \param orig detector to be copied
+    */
     interpolatingDetector(interpolatingDetector *orig)
         : singlePhotonDetector(orig) {
         // if (orig->interp)
@@ -66,10 +70,28 @@ class interpolatingDetector : public singlePhotonDetector {
         fi = orig->fi;
     }
 
+    /**
+     * copy constructor (deep copy)
+     * stricly, TODO: Implement Rule of Five!
+     * (copy op=, move ctor, and move op= would need to be defined)
+    */
+    interpolatingDetector(interpolatingDetector const& other)
+        : singlePhotonDetector(other) {
+
+        interp = other.interp;
+
+        id = other.id;
+        fi = other.fi;
+    }
+
     virtual interpolatingDetector *Clone() {
         return new interpolatingDetector(this);
     }
 
+    virtual interpolatingDetector *Copy() {
+        return new interpolatingDetector(*this);
+    }
+
     virtual int setId(int i) {
         id = i;
         //  interp->setId(id);

slsDetectorCalibration/jungfrauExecutables/CMakeLists.txt

@@ -6,18 +6,86 @@
 
 set(JUNGFRAU_EXECUTABLES)
 
-find_package(fmt REQUIRED)
+#find_package(fmt REQUIRED)
+
+#nlohmann_json
+#If the library was INSTALLED (i.e. sudo dnf install nlohmann-json3-dev), do stuff described in the repo under CMake -> External
+#find_package(nlohmann_json 3.2.0 REQUIRED)
+#find_package(nlohmann_json 3.11.2 REQUIRED)
+#find_package(nlohmann_json 3.11.3 REQUIRED)
+#
+#If the library was not installed but just cloned from https://github.com/nlohmann/json.git (possible, because it is a header-only library),
+# do stuff described in the repo under CMake -> Embedded
+#set(JSON_BuildTests OFF CACHE INTERNAL "")
+# If you only include this third party in PRIVATE source files, you do not
+# need to install it when your main project gets installed.
+# set(JSON_Install OFF CACHE INTERNAL "")
+#add_subdirectory(nlohmann_json) #Put the actual path to json
+#
+#Alternative: Fetch and install it from remote
+include(FetchContent)
+FetchContent_Declare(
+  json
+  GIT_REPOSITORY https://github.com/nlohmann/json.git
+  GIT_TAG v3.11.3  # Replace with the version you need
+)
+#FetchContent_Declare(json URL https://github.com/nlohmann/json/releases/download/v3.11.3/json.tar.xz) #Alternative (from the repo documentation)
+FetchContent_MakeAvailable(json)
+
+FetchContent_Declare(
+  fmt
+  GIT_REPOSITORY https://github.com/fmtlib/fmt.git
+  GIT_TAG master
+)
+FetchContent_MakeAvailable(fmt)
+
+# HDF5 file writing 
+if (SLS_USE_HDF5)
+    find_package(HDF5 1.10 COMPONENTS CXX REQUIRED)
+	    list (APPEND SOURCES 
+	        HDF5File.cpp 
+	    )
+endif (SLS_USE_HDF5)
+
+
 
 # jungfrauRawDataProcess
 add_executable(jungfrauRawDataProcess jungfrauRawDataProcess.cpp)
 target_compile_definitions(jungfrauRawDataProcess PRIVATE MODULE)
 list(APPEND JUNGFRAU_EXECUTABLES jungfrauRawDataProcess)
 
+# jungfrauRawDataProcessChipAldo
+add_executable(jungfrauRawDataProcessChipAldo jungfrauRawDataProcess_filetxt.cpp)
+target_compile_definitions(jungfrauRawDataProcessChipAldo PRIVATE CHIP ALDO)
+list(APPEND JUNGFRAU_EXECUTABLES jungfrauRawDataProcessChipAldo)
+
 # jungfrauRawDataProcessStrx
 add_executable(jungfrauRawDataProcessStrx jungfrauRawDataProcess_filetxt.cpp)
 target_compile_definitions(jungfrauRawDataProcessStrx PRIVATE JFSTRX)
 list(APPEND JUNGFRAU_EXECUTABLES jungfrauRawDataProcessStrx)
 
+# jungfrauRawDataProcessStrxQuad
+add_executable(jungfrauRawDataProcessStrxQuad jungfrauRawDataProcess_filetxt.cpp)
+target_compile_definitions(jungfrauRawDataProcessStrxQuad PRIVATE JFSTRXQ)
+list(APPEND JUNGFRAU_EXECUTABLES jungfrauRawDataProcessStrxQuad)
+
+# jungfrauRawDataProcessStrxQuadH5
+# HDF5
+if (SLS_USE_HDF5)
+    if (HDF5_FOUND)
+       add_executable(jungfrauRawDataProcessStrxQuadH5 jungfrauRawDataProcess_filetxtH5.cpp)
+       #target_compile_definitions(jungfrauRawDataProcessStrxQuadH5 PRIVATE JFSTRXQH5)
+       list(APPEND JUNGFRAU_EXECUTABLES jungfrauRawDataProcessStrxQuadH5)
+       target_include_directories(jungfrauRawDataProcessStrxQuadH5 PRIVATE ${HDF5_INCLUDE_DIRS}  ${CMAKE_INSTALL_PREFIX}/include)
+       target_link_libraries(jungfrauRawDataProcessStrxQuadH5 PRIVATE ${HDF5_LIBRARIES})
+
+       add_executable(jungfrauRawDataProcessStrxQuadH5SC jungfrauRawDataProcess_filetxtH5_SC.cpp)
+       #target_compile_definitions(jungfrauRawDataProcessStrxQuadH5 PRIVATE JFSTRXQH5)
+       list(APPEND JUNGFRAU_EXECUTABLES jungfrauRawDataProcessStrxQuadH5SC)
+       target_include_directories(jungfrauRawDataProcessStrxQuadH5SC PRIVATE ${HDF5_INCLUDE_DIRS}  ${CMAKE_INSTALL_PREFIX}/include)
+       target_link_libraries(jungfrauRawDataProcessStrxQuadH5SC PRIVATE ${HDF5_LIBRARIES})
+    endif ()
+endif (SLS_USE_HDF5)
 
 # jungfrauRawDataProcessStrxChip1
 add_executable(jungfrauRawDataProcessStrxChip1 jungfrauRawDataProcess.cpp)
@@ -58,6 +126,22 @@ list(APPEND JUNGFRAU_EXECUTABLES jungfrauRawDataProcessStrxOldAldo)
 # others to be added if needed (might already be there in Makefile.cluster_finder TO BE CHECKED)
 
 
+if (SLS_USE_HDF5)
+   if (HDF5_FOUND)
+        target_include_directories(jungfrauRawDataProcessStrxQuadH5 PRIVATE 
+								    ${HDF5_INCLUDE_DIRS}
+								    ${CMAKE_INSTALL_PREFIX}/include
+								    )
+	    target_link_libraries(jungfrauRawDataProcessStrxQuadH5 PRIVATE ${HDF5_LIBRARIES})
+
+        target_include_directories(jungfrauRawDataProcessStrxQuadH5SC PRIVATE 
+								    ${HDF5_INCLUDE_DIRS}
+								    ${CMAKE_INSTALL_PREFIX}/include
+								    )
+	    target_link_libraries(jungfrauRawDataProcessStrxQuadH5SC PRIVATE ${HDF5_LIBRARIES})
+
+   endif ()
+endif (SLS_USE_HDF5)
 
 
 foreach(exe ${JUNGFRAU_EXECUTABLES})
@@ -67,11 +151,16 @@ foreach(exe ${JUNGFRAU_EXECUTABLES})
         ../interpolations
         ../dataStructures 
         ../interpolations/etaVEL
-        ../../slsSupportLib/include/
+        ../../slsSupportLib/include/sls/
         ../../slsReceiverSoftware/include/
         ../tiffio/include
         ${fmt_INCLUDE_DIRS}
     )
+#    if (SLS_USE_HDF5)
+#       if (HDF5_FOUND)
+#       	  target_include_directories(${exe} PRIVATE ${HDF5_INCLUDE_DIRS}  ${CMAKE_INSTALL_PREFIX}/include)
+#       endif ()
+#    endif (SLS_USE_HDF5)
 
     target_link_libraries(${exe} 
         PUBLIC
@@ -79,6 +168,7 @@ foreach(exe ${JUNGFRAU_EXECUTABLES})
         pthread 
         tiffio
         fmt::fmt
+	    nlohmann_json::nlohmann_json
           #-L/usr/lib64/
         #-lm  -lstdc++  -lrt
 
@@ -86,7 +176,11 @@ foreach(exe ${JUNGFRAU_EXECUTABLES})
         slsProjectWarnings
         slsProjectOptions
     )
-
+#    if (SLS_USE_HDF5)
+#       if (HDF5_FOUND)
+#	target_link_libraries(${exe} PRIVATE ${HDF5_LIBRARIES})
+#       endif ()
+#    endif (SLS_USE_HDF5)
 
     set_target_properties(${exe}  PROPERTIES 
         RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin
@@ -98,4 +192,8 @@ foreach(exe ${JUNGFRAU_EXECUTABLES})
     
 endforeach(exe ${JUNGFRAU_EXECUTABLES})
 
-install(TARGETS ${JUNGFRAU_EXECUTABLES} DESTINATION bin)
+
+
+
+
+install(TARGETS ${JUNGFRAU_EXECUTABLES} DESTINATION bin)

slsDetectorCalibration/jungfrauExecutables/jungfrauRawDataProcess.cpp

@@ -11,7 +11,7 @@
 #define RAWDATA
 
 #if !defined JFSTRX && !defined JFSTRXOLD && !defined JFSTRXCHIP1 &&           \
-    !defined JFSTRXCHIP6
+    !defined JFSTRXCHIP6 && !defined CHIP
 #ifndef MODULE
 #include "jungfrauHighZSingleChipData.h"
 #endif
@@ -20,6 +20,10 @@
 #endif
 #endif
 
+#ifdef CHIP
+#include "jungfrauSingleChipData.h"
+#endif
+
 #ifdef JFSTRX
 #include "jungfrauLGADStrixelsData_new.h"
 #endif
@@ -41,6 +45,9 @@
 #include <ctime>
 #include <fmt/core.h>
 
+#include <nlohmann/json.hpp>
+using json = nlohmann::json;
+
 
 std::string getRootString( const std::string& filepath ) {
   size_t pos1 = filepath.find_last_of("/");
@@ -61,20 +68,24 @@ std::string getRootString( const std::string& filepath ) {
 std::string createFileName( const std::string& dir, const std::string& fprefix="run", const std::string& fsuffix="", const std::string& fext="raw", int aindex=0, int mindex=0, int findex=0, int outfilecounter=-1 ) {
   if (outfilecounter >= 0)
     return fmt::format("{:s}/{:s}_d{:d}_f{:d}_{:d}_f{:05d}.{:s}", dir, fprefix, mindex, findex, aindex, outfilecounter, fext);
-  else if (fsuffix.length()!=0)
-    return fmt::format("{:s}/{:s}_{:s}.{:s}", dir, fprefix, fsuffix, fext);
+  else if (fsuffix.length()!=0) {
+    if (fsuffix == "master")
+      return fmt::format("{:s}/{:s}_master_{:d}.{:s}", dir, fprefix, aindex, fext);
+    else
+      return fmt::format("{:s}/{:s}_{:s}.{:s}", dir, fprefix, fsuffix, fext);
+  }
   else
     return fmt::format("{:s}/{:s}_d{:d}_f{:d}_{:d}.{:s}", dir, fprefix, mindex, findex, aindex, fext);
 }
 
 int main(int argc, char *argv[]) {
 
-    if (argc < 5) {
+    if (argc < 6) {
         std::cout
             << "Usage is " << argv[0]
-            << "indir outdir fprefix(excluding slsDetector standard suffixes and extension) fextension "
-               "[runmin] [runmax] [pedfile (raw or tiff)] [threshold] "
-               "[nframes] [xmin xmax ymin ymax] [gainmap]"
+            << "indir outdir [fprefix(excluding slsDetector standard suffixes and extension)] [fextension] "
+	       "[fmin] [fmax] [runmin] [runmax] [pedfile (raw or tiff)] [threshold] "
+               "[nframes] [xmin xmax ymin ymax] [optional: bool read rxroi from data file header] [gainmap]"
             << std::endl;
         std::cout
             << "threshold <0 means analog; threshold=0 means cluster finder; "
@@ -106,94 +117,131 @@ int main(int argc, char *argv[]) {
     std::string outdir(argv[2]);
     std::string fprefix(argv[3]);
     std::string fext(argv[4]);
+
+    int fmin = 0;
+    if (argc >= 6)
+        fmin = atoi(argv[5]);
+    int fmax = fmin;
+    if (argc >= 7)
+        fmax = atoi(argv[6]);
+
     int runmin = 0;
-
     // cout << "argc is " << argc << endl;
-    if (argc >= 6) {
-        runmin = atoi(argv[5]);
+    if (argc >= 8) {
+        runmin = atoi(argv[7]);
     }
-
     int runmax = runmin;
-
-    if (argc >= 7) {
-        runmax = atoi(argv[6]);
+    if (argc >= 9) {
+        runmax = atoi(argv[8]);
     }
 
     std::string pedfilename{};
-    if (argc >= 8) {
-        pedfilename = argv[7];
+    if (argc >= 10) {
+        pedfilename = argv[9];
     }
     double thr = 0;
     double thr1 = 1;
 
-    if (argc >= 9) {
-        thr = atof(argv[8]);
+    if (argc >= 11) {
+        thr = atof(argv[10]);
     }
 
     int nframes = 0;
 
-    if (argc >= 10) {
-        nframes = atoi(argv[9]);
+    if (argc >= 12) {
+        nframes = atoi(argv[11]);
     }
+
+    bool readrxroifromdatafile = false;
+    if (argc >= 17)
+      readrxroifromdatafile = atoi(argv[16]);
     
+    // Receiver ROI
+    uint16_t rxroi_xmin = 0;
+    uint16_t rxroi_xmax = 0;
+    uint16_t rxroi_ymin = 0;
+    uint16_t rxroi_ymax = 0;
+
+    { //protective scope so ifstream gets destroyed properly
+
+      auto jsonmastername = createFileName( indir, fprefix, "master", "json", runmin );
+      std::cout << "json master file " << jsonmastername << std::endl;
+      std::ifstream masterfile(jsonmastername); //, ios::in | ios::binary);
+      if (masterfile.is_open()) {
+	json j;
+	masterfile >> j;
+	rxroi_xmin = j["Receiver Roi"]["xmin"];
+	rxroi_xmax = j["Receiver Roi"]["xmax"];
+	rxroi_ymin = j["Receiver Roi"]["ymin"];
+	rxroi_ymax = j["Receiver Roi"]["ymax"];
+	masterfile.close();
+	std::cout << "Read Receiver ROI [" << rxroi_xmin << ", " << rxroi_xmax << ", "
+		  << rxroi_ymin << ", " << rxroi_ymax << "] from json master file" << std::endl;
+      } else 
+	std::cout << "Could not open master file " << jsonmastername << std::endl;
+    
+    }
+
     // Define decoders...
 #if !defined JFSTRX && !defined JFSTRXOLD && !defined JFSTRXCHIP1 &&           \
-    !defined JFSTRXCHIP6
+    !defined JFSTRXCHIP6 && !defined CHIP
 #ifndef MODULE
     jungfrauHighZSingleChipData *decoder = new jungfrauHighZSingleChipData();
     int nx = 256, ny = 256;
 #endif
 #ifdef MODULE
-    jungfrauModuleData *decoder = new jungfrauModuleData();
+    jungfrauModuleData *decoder = new jungfrauModuleData(rxroi_xmin, rxroi_xmax, rxroi_ymin, rxroi_ymax);
     int nx = 1024, ny = 512;
 #endif
 #endif
 
+#ifdef CHIP
+    std::cout << "Jungfrau pixel module single chip readout" << std::endl;
+    jungfrauSingleChipData *decoder = new jungfrauSingleChipData();
+    int nx = 256, ny = 256;
+#endif
+
 #ifdef JFSTRX
-    cout << "Jungfrau strixel full module readout" << endl;
-    // ROI
-    uint16_t xxmin = 0;
-    uint16_t xxmax = 0;
-    uint16_t yymin = 0;
-    uint16_t yymax = 0;
+    std::cout << "Jungfrau strixel full module readout" << std::endl;
 
 #ifndef ALDO
-    { //THIS SCOPE IS IMPORTANT! (To ensure proper destruction of ifstream)
-      using header = sls::defs::sls_receiver_header;
-      // check if there is a roi in the header
-      typedef struct {
-        uint16_t xmin;
-        uint16_t xmax;
-        uint16_t ymin;
-        uint16_t ymax;
-      } receiverRoi_compact;
-      receiverRoi_compact croi;
-      std::string fsuffix{};
-      auto filename = createFileName( indir, fprefix, fsuffix, fext, runmin );
-      std::cout << "Reading header of file " << filename << " to check for ROI "
-		<< std::endl;
-      ifstream firstfile(filename, ios::in | ios::binary);
-      if (firstfile.is_open()) {
-        header hbuffer;
-        std::cout << "sizeof(header) = " << sizeof(header) << std::endl;
-        if (firstfile.read((char *)&hbuffer, sizeof(header))) {
-	  memcpy(&croi, &hbuffer.detHeader.detSpec1, 8);
-	  std::cout << "Read ROI [" << croi.xmin << ", " << croi.xmax << ", "
-		    << croi.ymin << ", " << croi.ymax << "]" << std::endl;
-	  xxmin = croi.xmin;
-	  xxmax = croi.xmax;
-	  yymin = croi.ymin;
-	  yymax = croi.ymax;
-        } else
-	  std::cout << "reading error" << std::endl;
-        firstfile.close();
-      } else
-        std::cout << "Could not open " << filename << " for reading " << std::endl;
-    } //end of protective scope
+    if (readrxroifromdatafile)
+      { //THIS SCOPE IS IMPORTANT! (To ensure proper destruction of ifstream)
+	using header = sls::defs::sls_receiver_header;
+	// check if there is a roi in the header
+	typedef struct {
+	  uint16_t xmin;
+	  uint16_t xmax;
+	  uint16_t ymin;
+	  uint16_t ymax;
+	} receiverRoi_compact;
+	receiverRoi_compact croi;
+	std::string fsuffix{};
+	auto filename = createFileName( indir, fprefix, fsuffix, fext, runmin );
+	std::cout << "Reading header of file " << filename << " to check for ROI "
+		  << std::endl;
+	ifstream firstfile(filename, ios::in | ios::binary);
+	if (firstfile.is_open()) {
+	  header hbuffer;
+	  std::cout << "sizeof(header) = " << sizeof(header) << std::endl;
+	  if (firstfile.read((char *)&hbuffer, sizeof(header))) {
+	    memcpy(&croi, &hbuffer.detHeader.detSpec1, 8);
+	    std::cout << "Read ROI [" << croi.xmin << ", " << croi.xmax << ", "
+		      << croi.ymin << ", " << croi.ymax << "]" << std::endl;
+	    rxroi_xmin = croi.xmin;
+	    rxroi_xmax = croi.xmax;
+	    rxroi_ymin = croi.ymin;
+	    rxroi_ymax = croi.ymax;
+	  } else
+	    std::cout << "reading error" << std::endl;
+	  firstfile.close();
+	} else
+	  std::cout << "Could not open " << filename << " for reading " << std::endl;
+      } //end of protective scope
 #endif
 
     jungfrauLGADStrixelsData *decoder =
-        new jungfrauLGADStrixelsData(xxmin, xxmax, yymin, yymax);
+      new jungfrauLGADStrixelsData(rxroi_xmin, rxroi_xmax, rxroi_ymin, rxroi_ymax);
     int nx = 1024 / 3, ny = 512 * 5;
 #endif
 #ifdef JFSTRXCHIP1
@@ -218,19 +266,20 @@ int main(int argc, char *argv[]) {
     decoder->getDetectorSize(nx, ny);
     std::cout << "Detector size is " << nx << " " << ny << std::endl;
 
-    int xmin = 0, xmax = nx, ymin = 0, ymax = ny;
-    if (argc >= 14) {
-        xmin = atoi(argv[10]);
-        xmax = atoi(argv[11]);
-        ymin = atoi(argv[12]);
-        ymax = atoi(argv[13]);
+    //Cluster finder ROI
+    int xmin = 0, xmax = nx-1, ymin = 0, ymax = ny-1;
+    if (argc >= 16) {
+        xmin = atoi(argv[12]);
+        xmax = atoi(argv[13]);
+        ymin = atoi(argv[14]);
+        ymax = atoi(argv[15]);
     }
-    std::cout << xmin << " " << xmax << " " << ymin << " " << ymax << " "
+    std::cout << "Cluster finder ROI: [" << xmin << ", " << xmax << ", " << ymin << ", " << ymax << "]"
               << std::endl;
 
     char *gainfname = NULL;
-    if (argc > 14) {
-        gainfname = argv[14];
+    if (argc > 17) {
+        gainfname = argv[17];
         std::cout << "Gain map file name is: " << gainfname << std::endl;
     }
 
@@ -239,6 +288,8 @@ int main(int argc, char *argv[]) {
     std::cout << "input directory is " << indir << std::endl;
     std::cout << "output directory is " << outdir << std::endl;
     std::cout << "input file prefix is " << fprefix << std::endl;
+    std::cout << "fmin is " << fmin << std::endl;
+    std::cout << "fmax is " << fmax << std::endl;
     std::cout << "runmin is " << runmin << std::endl;
     std::cout << "runmax is " << runmax << std::endl;
     if (pedfilename.length()!=0)
@@ -319,7 +370,7 @@ int main(int argc, char *argv[]) {
 
 	  mt->setFrameMode(ePedestal);
 
-	  ifstream pedefile(fname, ios::in | ios::binary);
+	  std::ifstream pedefile(fname, ios::in | ios::binary);
 	  //      //open file
 	  if (pedefile.is_open()) {
 	    std::cout << "bbbb " << std::ctime(&end_time) << std::endl;
@@ -380,26 +431,27 @@ int main(int argc, char *argv[]) {
     }
 
     ifr = 0;
-    int ifile = 0;
+    int ioutfile = 0;
 
     mt->setFrameMode(eFrame);
 
     FILE *of = NULL;
     
-    for (int irun = runmin; irun <= runmax; irun++) {
+    for (int irun = runmin; irun <= runmax; ++irun) {
+      for (int ifile = fmin; ifile <= fmax; ++ifile) {
         std::cout << "DATA ";
 	std::string fsuffix{};
-	auto fname = createFileName( indir, fprefix, fsuffix, fext, irun );
-	auto imgfname = createFileName( outdir, fprefix, fsuffix, "tiff", irun );
-	auto cfname = createFileName( outdir, fprefix, fsuffix, "clust", irun );
+	auto fname = createFileName( indir, fprefix, fsuffix, fext, irun, 0, ifile );
+	auto imgfname = createFileName( outdir, fprefix, fsuffix, "tiff", irun, 0, ifile );
+	auto cfname = createFileName( outdir, fprefix, fsuffix, "clust", irun, 0, ifile );
         std::cout << fname << " ";
         std::cout << imgfname << std::endl;
         std::time(&end_time);
         std::cout << std::ctime(&end_time) << std::endl;
         //  std::cout <<  fname << " " << outfname << " " << imgfname <<  std::endl;
-        ifstream filebin(fname, ios::in | ios::binary);
+	std::ifstream filebin(fname, ios::in | ios::binary);
         //      //open file
-        ifile = 0;
+        ioutfile = 0;
         if (filebin.is_open()) {
             if (thr <= 0 && cf != 0) { // cluster finder
                 if (of == NULL) {
@@ -436,10 +488,10 @@ int main(int argc, char *argv[]) {
                     std::cout << " " << ifr << " " << ff << std::endl;
                 if (nframes > 0) {
                     if (ifr % nframes == 0) {
-		      imgfname = createFileName( outdir, fprefix, fsuffix, "tiff", irun, 0, 0, ifile );
+		      imgfname = createFileName( outdir, fprefix, fsuffix, "tiff", irun, 0, 0, ioutfile );
 		      mt->writeImage(imgfname.c_str(), thr1);
                         mt->clearImage();
-                        ifile++;
+                        ioutfile++;
                     }
                 }
                 // } else
@@ -453,7 +505,7 @@ int main(int argc, char *argv[]) {
             }
             if (nframes >= 0) {
                 if (nframes > 0)
-		  imgfname = createFileName( outdir, fprefix, fsuffix, "tiff", irun, 0, 0, ifile );
+		  imgfname = createFileName( outdir, fprefix, fsuffix, "tiff", irun, 0, 0, ioutfile );
                 std::cout << "Writing tiff to " << imgfname << " " << thr1
                           << std::endl;
                 mt->writeImage(imgfname.c_str(), thr1);
@@ -469,9 +521,10 @@ int main(int argc, char *argv[]) {
         } else
             std::cout << "Could not open " << fname << " for reading "
                       << std::endl;
+      }
     }
     if (nframes < 0) {
-      auto imgfname = createFileName( outdir, fprefix, "sum", "tiff", -1, 0, 0, -1 );
+      auto imgfname = createFileName( outdir, fprefix, "sum", "tiff", runmin, 0, fmin, -1 );
         std::cout << "Writing tiff to " << imgfname << " " << thr1 << std::endl;
         mt->writeImage(imgfname.c_str(), thr1);
     }

slsDetectorCalibration/jungfrauExecutables/jungfrauRawDataProcess_filetxt.cpp

@@ -10,8 +10,8 @@
 
 #define RAWDATA
 
-#if !defined JFSTRX && !defined JFSTRXOLD && !defined JFSTRXCHIP1 &&           \
-    !defined JFSTRXCHIP6
+#if !defined JFSTRX && !defined JFSTRXQ && !defined JFSTRXOLD && !defined JFSTRXCHIP1 &&           \
+    !defined JFSTRXCHIP6 && !defined CHIP
 #ifndef MODULE
 #include "jungfrauHighZSingleChipData.h"
 #endif
@@ -20,9 +20,16 @@
 #endif
 #endif
 
+#ifdef CHIP
+#include "jungfrauSingleChipData.h"
+#endif
+
 #ifdef JFSTRX
 #include "jungfrauLGADStrixelsData_new.h"
 #endif
+#ifdef JFSTRXQ
+#include "jungfrauLGADStrixelsDataQuad.h"
+#endif
 #if defined JFSTRXCHIP1 || defined JFSTRXCHIP6
 #include "jungfrauLGADStrixelsDataSingleChip.h"
 #endif
@@ -41,6 +48,9 @@
 #include <ctime>
 #include <fmt/core.h>
 
+#include <nlohmann/json.hpp>
+using json = nlohmann::json;
+
 
 std::string getRootString( const std::string& filepath ) {
   size_t pos1;
@@ -71,11 +81,11 @@ std::string createFileName( const std::string& dir, const std::string& fprefix="
 //NOTE THAT THE DATA FILES HAVE TO BE IN THE RIGHT ORDER SO THAT PEDESTAL TRACKING WORKS!
 int main(int argc, char *argv[]) {
 
-    if (argc < 10) {
+    if (argc < 11) {
         std::cout
             << "Usage is " << argv[0]
-            << " filestxt outdir [pedfile (raw or tiff)] [xmin xmax ymin ymax] "
-               "[threshold] [nframes] "
+            << " filestxt outdir [json master] [pedfile (raw or tiff)] [xmin xmax ymin ymax] "
+               "[threshold] [nframes] [optional: bool read rxroi from data file header]"
 	       "NOTE THAT THE DATA FILES HAVE TO BE IN THE RIGHT ORDER SO THAT PEDESTAL TRACKING WORKS! "
             << std::endl;
         std::cout
@@ -105,19 +115,19 @@ int main(int argc, char *argv[]) {
 
     const std::string txtfilename(argv[1]);
     const std::string outdir(argv[2]);
-    const std::string pedfilename(argv[3]);
-
-    int xmin = atoi(argv[4]);
-    int xmax = atoi(argv[5]);
-    int ymin = atoi(argv[6]);
-    int ymax = atoi(argv[7]);
+    const std::string jsonmastername(argv[3]);
+    const std::string pedfilename(argv[4]);
 
     double thr = 0;
     double thr1 = 1;
-    thr = atof(argv[8]);
+    thr = atof(argv[9]);
 
     int nframes = 0;
-    nframes = atoi(argv[9]);
+    nframes = atoi(argv[10]);
+
+    bool readrxroifromdatafile = false;
+    if (argc > 11)
+      readrxroifromdatafile = atoi(argv[11]);
 
     //Get vector of filenames from input txt-file
     std::vector<std::string> filenames{};
@@ -146,10 +156,34 @@ int main(int argc, char *argv[]) {
     }
 
     std::cout << "###############" << std::endl;
-    
+
+    // Receiver ROI
+    uint16_t rxroi_xmin = 0;
+    uint16_t rxroi_xmax = 0;
+    uint16_t rxroi_ymin = 0;
+    uint16_t rxroi_ymax = 0;
+
+    { //protective scope so ifstream gets destroyed properly
+
+      std::ifstream masterfile(jsonmastername); //, ios::in | ios::binary);
+      if (masterfile.is_open()) {
+	json j;
+	masterfile >> j;
+	rxroi_xmin = j["Receiver Roi"]["xmin"];
+	rxroi_xmax = j["Receiver Roi"]["xmax"];
+	rxroi_ymin = j["Receiver Roi"]["ymin"];
+	rxroi_ymax = j["Receiver Roi"]["ymax"];
+	masterfile.close();
+	std::cout << "Read rxROI [" << rxroi_xmin << ", " << rxroi_xmax << ", "
+		      << rxroi_ymin << ", " << rxroi_ymax << "]" << std::endl;
+      } else 
+	std::cout << "Could not open master file " << jsonmastername << std::endl;
+      
+    }
+
     // Define decoders...
-#if !defined JFSTRX && !defined JFSTRXOLD && !defined JFSTRXCHIP1 &&           \
-    !defined JFSTRXCHIP6
+#if !defined JFSTRX && !defined JFSTRXQ && !defined JFSTRXOLD && !defined JFSTRXCHIP1 &&           \
+    !defined JFSTRXCHIP6 && !defined CHIP
 #ifndef MODULE
     jungfrauHighZSingleChipData *decoder = new jungfrauHighZSingleChipData();
     int nx = 256, ny = 256;
@@ -160,52 +194,60 @@ int main(int argc, char *argv[]) {
 #endif
 #endif
 
+#ifdef CHIP
+    std::cout << "Jungfrau pixel module single chip readout" << std::endl;
+    jungfrauSingleChipData *decoder = new jungfrauSingleChipData();
+    int nx = 256, ny = 256;
+#endif
+
 #ifdef JFSTRX
-    cout << "Jungfrau strixel full module readout" << endl;
-    // ROI
-    uint16_t xxmin = 0;
-    uint16_t xxmax = 0;
-    uint16_t yymin = 0;
-    uint16_t yymax = 0;
+    std::cout << "Jungfrau strixel full module readout" << std::endl;
 
 #ifndef ALDO
-    { //THIS SCOPE IS IMPORTANT! (To ensure proper destruction of ifstream)
-      using header = sls::defs::sls_receiver_header;
-      // check if there is a roi in the header
-      typedef struct {
-        uint16_t xmin;
-        uint16_t xmax;
-        uint16_t ymin;
-        uint16_t ymax;
-      } receiverRoi_compact;
-      receiverRoi_compact croi;
-      //std::string filepath(argv[9]); //This is a problem if the input files have different ROIs!
-      std::cout << "Reading header of file " << filenames[0] << " to check for ROI "
-		<< std::endl;
-      ifstream firstfile(filenames[0], ios::in | ios::binary);
-      if (firstfile.is_open()) {
-        header hbuffer;
-        std::cout << "sizeof(header) = " << sizeof(header) << std::endl;
-        if (firstfile.read((char *)&hbuffer, sizeof(header))) {
-	  memcpy(&croi, &hbuffer.detHeader.detSpec1, 8);
-	  std::cout << "Read ROI [" << croi.xmin << ", " << croi.xmax << ", "
-		    << croi.ymin << ", " << croi.ymax << "]" << std::endl;
-	  xxmin = croi.xmin;
-	  xxmax = croi.xmax;
-	  yymin = croi.ymin;
-	  yymax = croi.ymax;
-        } else
-	  std::cout << "reading error" << std::endl;
-        firstfile.close();
-      } else
-        std::cout << "Could not open " << filenames[0] << " for reading " << std::endl;
-    } //end of protective scope
+    if (readrxroifromdatafile)
+      { //THIS SCOPE IS IMPORTANT! (To ensure proper destruction of ifstream)
+	using header = sls::defs::sls_receiver_header;
+	// check if there is a roi in the header
+	typedef struct {
+	  uint16_t xmin;
+	  uint16_t xmax;
+	  uint16_t ymin;
+	  uint16_t ymax;
+	} receiverRoi_compact;
+	receiverRoi_compact croi;
+	//std::string filepath(argv[9]); //This is a problem if the input files have different ROIs!
+	std::cout << "Reading header of file " << filenames[0] << " to check for ROI "
+		  << std::endl;
+	std::ifstream firstfile( filenames[0], ios::in | ios::binary);
+	if (firstfile.is_open()) {
+	  header hbuffer;
+	  std::cout << "sizeof(header) = " << sizeof(header) << std::endl;
+	  if (firstfile.read((char *)&hbuffer, sizeof(header))) {
+	    memcpy(&croi, &hbuffer.detHeader.detSpec1, 8);
+	    std::cout << "Read ROI [" << croi.xmin << ", " << croi.xmax << ", "
+		      << croi.ymin << ", " << croi.ymax << "]" << std::endl;
+	    rxroi_xmin = croi.xmin;
+	    rxroi_xmax = croi.xmax;
+	    rxroi_ymin = croi.ymin;
+	    rxroi_ymax = croi.ymax;
+	  } else
+	    std::cout << "reading error" << std::endl;
+	  firstfile.close();
+	} else
+	  std::cout << "Could not open " << filenames[0] << " for reading " << std::endl;
+      } //end of protective scope
 #endif
 
     jungfrauLGADStrixelsData *decoder =
-        new jungfrauLGADStrixelsData(xxmin, xxmax, yymin, yymax);
+        new jungfrauLGADStrixelsData(rxroi_xmin, rxroi_xmax, rxroi_ymin, rxroi_ymax);
     int nx = 1024 / 3, ny = 512 * 5;
 #endif
+#ifdef JFSTRXQ
+    std::cout << "Jungfrau strixel quad" << std::endl;
+    jungfrauLGADStrixelsDataQuad *decoder =
+        new jungfrauLGADStrixelsDataQuad(rxroi_xmin, rxroi_xmax, rxroi_ymin, rxroi_ymax);
+    int nx = 1024 / 3, ny = 512 * 3;
+#endif
 #ifdef JFSTRXCHIP1
     std::cout << "Jungfrau strixel LGAD single chip 1" << std::endl;
     jungfrauLGADStrixelsDataSingleChip *decoder =
@@ -228,7 +270,16 @@ int main(int argc, char *argv[]) {
     decoder->getDetectorSize(nx, ny);
     std::cout << "Detector size is " << nx << " " << ny << std::endl;
 
+    //Cluster finder ROI
+    int xmin = 0, xmax = nx-1, ymin = 0, ymax = ny-1;
+    xmin = atoi(argv[5]);
+    xmax = atoi(argv[6]);
+    ymin = atoi(argv[7]);
+    ymax = atoi(argv[8]);
+    std::cout << "Cluster finder ROI: [" << xmin << ", " << xmax << ", " << ymin << ", " << ymax << "]"
+              << std::endl;
 
+    /* old
     if ( xmin == xmax ) {
       xmin = 0;
       xmax = nx;
@@ -239,6 +290,7 @@ int main(int argc, char *argv[]) {
     }
     std::cout << xmin << " " << xmax << " " << ymin << " " << ymax << " "
               << std::endl;
+    */
 
     /*
     char *gainfname = NULL;
@@ -410,7 +462,7 @@ int main(int argc, char *argv[]) {
       std::time(&end_time);
       std::cout << std::ctime(&end_time) << std::endl;
 
-      ifstream filebin(filenames[ifile], ios::in | ios::binary);
+      std::ifstream filebin(filenames[ifile], ios::in | ios::binary);
       //      //open file
       ioutfile = 0;
       if (filebin.is_open()) {

slsDetectorCalibration/jungfrauExecutables/jungfrauRawDataProcess_filetxtH5.cpp

@@ -0,0 +1,457 @@
+// SPDX-License-Identifier: LGPL-3.0-or-other
+// Copyright (C) 2021 Contributors to the SLS Detector Package
+// #include "sls/ansi.h"
+//#include <iostream>
+#undef CORR
+
+#define C_GHOST 0.0004
+
+#define CM_ROWS 50
+
+#define RAWDATA
+
+
+#include "jungfrauLGADStrixelsDataQuadH5.h"
+
+#include "multiThreadedCountingDetector.h"
+#include "singlePhotonDetector.h"
+
+#include <fstream>
+#include <map>
+#include <memory>
+#include <stdio.h>
+#include <sys/stat.h>
+
+#include <ctime>
+#include <fmt/core.h>
+
+/*
+#include <nlohmann/json.hpp>
+using json = nlohmann::json;
+*/
+
+/*Dataset paths according to different beamlines*/
+std::string const data_datasetname_furka("/data/JF18T01V01/data");
+std::string const index_datasetname_furka("/data/JF18T01V01/frame_index");
+std::string const data_datasetname_xfelSCS("/INSTRUMENT/SCS_HRIXS_JUNGF/DET/JNGFR01:daqOutput/data/adc");
+std::string const index_datasetname_xfelSCS("/INSTRUMENT/SCS_HRIXS_JUNGF/DET/JNGFR01:daqOutput/data/frameNumber");
+
+std::string getRootString( std::string const& filepath ) {
+  size_t pos1;
+  if (filepath.find("/") == std::string::npos )
+    pos1 = 0;
+  else
+    pos1 = filepath.find_last_of("/")+1;
+  size_t pos2 = filepath.find_last_of(".");
+  //std::cout << "pos1 " << pos1 << " pos2 " << pos2 << " size " << filepath.length() << std::endl;
+  return filepath.substr( pos1, pos2-pos1 );
+}
+
+//Create file name string
+//   dir:     directory
+//   fprefix:   fileprefix (without extension)
+//   fsuffix:   filesuffix (for output files, e.g. "ped")
+//   fext:    file extension (e.g. "raw")
+std::string createFileName( std::string const& dir, std::string const& fprefix="run",
+			    std::string const& fsuffix="", std::string const& fext="raw", int const outfilecounter=-1 ) {
+  std::string return_string;
+  if (outfilecounter >= 0)
+    return_string = fmt::format("{:s}/{:s}_{:s}_f{:05d}", dir, fprefix, fsuffix, outfilecounter);
+  else if (fsuffix.length()!=0)
+    return_string = fmt::format("{:s}/{:s}_{:s}", dir, fprefix, fsuffix);
+  else
+    return_string = fmt::format("{:s}/{:s}", dir, fprefix);
+
+  if (fext.length()!=0)
+    return_string += "." + fext;
+
+  return return_string;
+}
+
+
+//NOTE THAT THE DATA FILES HAVE TO BE IN THE RIGHT ORDER SO THAT PEDESTAL TRACKING WORKS!
+int main(int argc, char *argv[]) {
+
+    if (argc < 4) {
+        std::cout
+            << "Usage is " << argv[0]
+            << " filestxt outdir [pedfile (h5)] "
+               " optional: [int dataset path; 0 means Furka, 1 means XFEL; overwrites default given in HDF5File.h] "
+               " [bool validate h5 rank] "
+               " [xmin xmax ymin ymax] [threshold] [nframes] "
+	             " NOTE THAT THE DATA FILES HAVE TO BE IN THE RIGHT ORDER SO THAT PEDESTAL TRACKING WORKS! "
+            << std::endl;
+        std::cout
+            << "threshold <0 means analog; threshold=0 means cluster finder; "
+               "threshold>0 means photon counting"
+            << std::endl;
+        std::cout
+            << "nframes <0 means sum everything; nframes=0 means one file per "
+               "run; nframes>0 means one file every nframes"
+            << std::endl;
+        return 1;
+    }
+
+    int const fifosize = 100; //1000;
+    int const nthreads = 10;
+    int const csize = 3; // 3
+    int const nsigma = 5;
+    int const nped = 10000;
+
+    int cf = 0;
+
+    std::string const txtfilename(argv[1]);
+    std::string const outdir(argv[2]);
+    std::string const pedfilename(argv[3]);
+
+    std::string datasetpath{};
+    std::string frameindexpath{};
+    if (argc > 4) {
+      switch (atoi(argv[4])) {
+      case 0:
+        datasetpath = data_datasetname_furka;
+        frameindexpath = index_datasetname_furka;
+        break;
+      case 1:
+        datasetpath = data_datasetname_xfelSCS;
+        frameindexpath = index_datasetname_xfelSCS;
+        break;
+      default:
+        break;
+      }
+    }
+
+    bool validate_rank=true;
+    if (argc > 5)
+      validate_rank = atoi(argv[5]);
+
+    double thr = 0;
+    double thr1 = 1;
+    if (argc > 9)
+      thr = atof(argv[9]);
+
+    int nframes = 0;
+    if (argc > 10)
+      nframes = atoi(argv[10]);
+
+    //Get vector of filenames from input txt-file
+    std::vector<std::string> filenames{};
+    //filenames.reserve(512);
+    { //Safety scope for ifstream
+      ifstream inputs( txtfilename, std::ios::in );
+      if (inputs.is_open()) {
+	    std::cout << "Reading imput filenames from txt-file ..." << std::endl;
+	    std::string line{};
+	    while (!inputs.eof()) {
+	      std::getline(inputs, line);
+        if(line.find(".h5") != std::string::npos) {
+	        filenames.emplace_back(line);
+	        std::cout << line << std::endl; //" line.max_size() " << line.max_size() << " filenames.capacity() " << filenames.capacity() << '\n';
+        }
+	    }
+	    inputs.close();
+	    std::cout << "---- Reached end of txt-file. ----" << std::endl;
+      } else
+	      std::cout << "Could not open " << txtfilename << std::endl;
+      if (filenames.size()>0) {
+	      std::cout << filenames.size() << " filenames found in " << txtfilename << std::endl;
+	      std::cout << "The files will be processed in the same order as found in the txt-file." << std::endl;
+      } else {
+	      std::cout << "No files found in txt-file!" << std::endl;
+	      return 1;
+      }
+    }
+
+    std::cout << "###############" << std::endl;
+
+    // Define decoder
+    std::cout << "Jungfrau strixel quad h5" << std::endl;
+    jungfrauLGADStrixelsDataQuadH5* decoder = new jungfrauLGADStrixelsDataQuadH5();
+    //auto decoder = std::make_unique<jungfrauLGADStrixelsDataQuadH5>();
+    int nx = 1024 / 3, ny = 512 * 3;
+
+    //Cluster finder ROI
+    int xmin = 0, xmax = nx-1, ymin = 0, ymax = ny-1;
+    if (argc > 9) {
+      xmin = atoi(argv[6]);
+      xmax = atoi(argv[7]);
+      ymin = atoi(argv[8]);
+      ymax = atoi(argv[9]);
+    }
+    std::cout << "Cluster finder ROI: [" << xmin << ", " << xmax << ", " << ymin << ", " << ymax << "]"
+              << std::endl;
+
+    decoder->getDetectorSize(nx, ny);
+    std::cout << "Detector size is " << nx << " " << ny << std::endl;
+
+    std::time_t end_time;
+
+    std::cout << "output directory is " << outdir << std::endl;
+    if (pedfilename.length()!=0)
+        std::cout << "pedestal file is " << pedfilename << std::endl;
+    if (thr > 0)
+        std::cout << "threshold is " << thr << std::endl;
+    std::cout << "Nframes is " << nframes << std::endl;
+
+    uint32_t nnx, nny;
+
+    singlePhotonDetector* filter =
+      new singlePhotonDetector(decoder, 3, nsigma, 1, NULL, nped, 200, -1, -1, NULL, NULL);
+    //auto filter = std::make_unique<singlePhotonDetector>(decoder.get(), 3, nsigma, 1, nullptr, nped, 200, -1, -1, nullptr, nullptr);
+
+    thr = 0.15 * thr;
+    //filter->newDataSet(); //This only initializes the dataset for the first thread (the other threads are created via cloning)
+    // int dsize = decoder->getDataSize();
+
+    if (thr > 0) {
+        std::cout << "threshold is " << thr << std::endl;
+        filter->setThreshold(thr);
+        cf = 0;
+
+    } else
+        cf = 1;
+
+    filter->setROI(xmin, xmax, ymin, ymax);
+    std::time(&end_time);
+    std::cout << std::ctime(&end_time) << std::endl;
+
+    char* buff;
+
+    multiThreadedCountingDetector* mt =
+        new multiThreadedCountingDetector(filter, nthreads, fifosize);
+    //auto mt = std::make_unique<multiThreadedCountingDetector>(filter.get(), nthreads, fifosize);
+    mt->setClusterSize(csize, csize);
+    mt->newDataSet(); //Initialize new dataset for each thread
+
+#ifndef ANALOG
+    mt->setDetectorMode(ePhotonCounting);
+    std::cout << "Counting!" << std::endl;
+    if (thr > 0) {
+        cf = 0;
+    }
+#endif
+//{
+#ifdef ANALOG
+    mt->setDetectorMode(eAnalog);
+    std::cout << "Analog!" << std::endl;
+    cf = 0;
+    // thr1=thr;
+#endif
+    //  }
+
+    mt->StartThreads();
+    mt->popFree(buff);
+
+    int ifr = 0; //frame counter of while loop
+    int framenumber = 0; //framenumber as read from file (detector)
+    std::vector<hsize_t> h5offset(1,0); //frame counter internal to HDF5File::ReadImage (provided for sanity check/debugging)
+
+    if (pedfilename.length()>1) {
+
+      std::cout << "PEDESTAL " << std::endl;
+
+      if (pedfilename.find(".tif") == std::string::npos) { //not a tiff file
+	      std::string const fname(pedfilename);
+	      std::cout << fname << std::endl;
+	      std::time(&end_time);
+	      std::cout << "aaa " << std::ctime(&end_time) << std::endl;
+
+	      mt->setFrameMode(ePedestal);
+
+	      //HDF5File pedefile;
+	      auto pedefile = std::make_unique<HDF5File>();
+        pedefile->SetFrameIndexPath(frameindexpath);
+        pedefile->SetImageDataPath(datasetpath);
+	      //      //open file
+	      if ( pedefile->OpenResources(fname.c_str(),validate_rank) ) {
+	        std::cout << "bbbb " << std::ctime(&end_time) << std::endl;
+	    
+	        framenumber = 0;
+
+	        while ( decoder->readNextFrame(*pedefile, framenumber, h5offset, buff) ) {
+
+	          if ((ifr + 1) % 100 == 0) {
+	            std::cout
+		            << " ****"
+		            << decoder->getValue(buff, 20, 20); // << std::endl;
+	          }
+	          mt->pushData(buff);
+	          mt->nextThread();
+	          mt->popFree(buff);
+	          ++ifr;
+	          if (ifr % 100 == 0) {
+	            std::cout << " ****" << ifr << " " << framenumber << " " << h5offset[0]
+			                  << std::endl;
+	          } // else
+	      
+	          if (ifr >= 1000)
+	            break;
+	          //framenumber = 0;
+	        }
+	    
+	        pedefile->CloseResources();
+	        while (mt->isBusy()) {
+	          ;
+	        }
+
+	        std::cout << "Writing pedestal to " << getRootString(pedfilename) << "_ped.tiff" << std::endl;
+	        auto imgfname = createFileName( outdir, getRootString(pedfilename), "ped", "");
+	        mt->writePedestal(imgfname.c_str());
+          std::cout << "Writing pedestal rms to " << getRootString(pedfilename) << "_rms.tiff" << std::endl;
+	        imgfname = createFileName( outdir, getRootString(pedfilename), "rms", "");
+	        mt->writePedestalRMS(imgfname.c_str());
+
+	      } else
+	        std::cout << "Could not open pedestal file " << fname
+		                << " for reading " << std::endl;
+
+      } else { //is a tiff file
+
+	      std::vector<double> ped(nx * ny);
+	      float* pp = ReadFromTiff(pedfilename.c_str(), nny, nnx);
+	      if (pp && (int)nnx == nx && (int)nny == ny) {
+	        for (int i = 0; i < nx * ny; i++) {
+	          ped[i] = pp[i];
+	        }
+	        delete[] pp;
+	        mt->setPedestal(ped.data());
+	        std::cout << "Pedestal set from tiff file " << pedfilename
+		                << std::endl;
+	      } else {
+	        std::cout << "Could not open pedestal tiff file " << pedfilename
+		                << " for reading " << std::endl;
+	      }
+      }
+      std::time(&end_time);
+      std::cout << std::ctime(&end_time) << std::endl;
+    }
+    
+    ifr = 0;
+    int ioutfile = 0;
+
+    mt->setFrameMode(eFrame);
+
+    FILE* of = nullptr;
+    
+    //NOTE THAT THE DATA FILES HAVE TO BE IN THE RIGHT ORDER SO THAT PEDESTAL TRACKING WORKS!
+    for (unsigned int ifile = 0; ifile != filenames.size(); ++ifile) {
+      std::cout << "DATA " << filenames[ifile] << " " << std::endl;
+      auto imgfname( createFileName( outdir, getRootString(filenames[ifile]), "", "" ) );
+      std::string const cfname( createFileName( outdir, getRootString(filenames[ifile]), "", "clust" ) );
+      
+      std::time(&end_time);
+      std::cout << std::ctime(&end_time) << std::endl;
+
+      //HDF5File fileh5;
+      auto fileh5 = std::make_unique<HDF5File>();
+      fileh5->SetFrameIndexPath(frameindexpath);
+      fileh5->SetImageDataPath(datasetpath);
+      //      //open file
+      ioutfile = 0;
+      if ( fileh5->OpenResources(filenames[ifile].c_str(), validate_rank) ) {
+	      if (thr <= 0 && cf != 0) { // cluster finder
+	        if (of == nullptr) {
+	          of = fopen(cfname.c_str(), "w");
+	          if (of) {
+	            if (mt) {
+		            mt->setFilePointer(of);
+		            std::cout << "file pointer set " << std::endl;
+	            } else {
+		            std::cerr << "Error: mt is null." << std::endl;
+		            return 1;
+	            }
+	            //mt->setFilePointer(of);
+	            //std::cout << "file pointer set " << std::endl;
+	            //std::cout << "Here! " << framenumber << " ";
+	          } else {
+	            std::cout << "Could not open " << cfname
+			                  << " for writing " << std::endl;
+	            mt->setFilePointer(nullptr);
+	            return 1;
+	          }
+	        }
+	      }
+	
+	//     //while read frame
+	      framenumber = 0;
+	      h5offset[0] = 0;
+	      ifr = 0;
+	      //std::cout << "Here! " << framenumber << " ";
+	      while ( decoder->readNextFrame(*fileh5, framenumber, h5offset, buff) ) {
+	        //std::cout << "Here! " << framenumber << " ";
+	        //         //push
+	        if ((ifr + 1) % 1000 == 0) {
+	          std::cout << " ****"
+		                  << decoder->getValue(buff, 20, 20); // << std::endl;
+	        }
+	        mt->pushData(buff);
+
+	        // 	//         //pop
+	        mt->nextThread();
+	        mt->popFree(buff); /* In the last execution of the loop,
+                              * this leaves buff outside of the Fifo!
+                              * Free explicitely at the end! */
+	  
+	        ++ifr;
+	        if (ifr % 1000 == 0)
+	          std::cout << " " << ifr << " " << framenumber << " " << h5offset[0]
+                      << std::endl;
+	        if (nframes > 0) {
+	          if (ifr % nframes == 0) {
+	            imgfname = createFileName( outdir, getRootString(filenames[ifile]), "", "", ioutfile );
+	            mt->writeImage(imgfname.c_str(), thr1);
+	            mt->clearImage();
+	            ++ioutfile;
+	          }
+	        }
+
+	        //framenumber = 0;
+	      }
+
+	      //std::cout << "aa --" << std::endl;
+	      fileh5->CloseResources();
+
+	      //std::cout << "bb --" << std::endl;
+	      while (mt->isBusy()) {
+	        ;
+	      }
+
+	      //std::cout << "cc --" << std::endl;
+	      if (nframes >= 0) {
+	        if (nframes > 0)
+	          imgfname = createFileName( outdir, getRootString(filenames[ifile]), "", "", ioutfile );
+	        std::cout << "Writing tiff to " << imgfname << " " << thr1 << ".tiff"
+		                << std::endl;
+	        mt->writeImage(imgfname.c_str(), thr1);
+	        mt->clearImage();
+	        if (of) {
+	          fclose(of);
+	          of = nullptr;
+	          mt->setFilePointer(nullptr);
+	        }
+	      }
+	      std::time(&end_time);
+	      std::cout << std::ctime(&end_time) << std::endl;
+      } else
+	      std::cout << "Could not open " << filenames[ifile] << " for reading "
+		              << std::endl;
+    }
+    if (nframes < 0) {     
+      //std::string fprefix( getRootString(filenames[0]) ); //Possibly, non-ideal name choice for file
+      auto imgfname( createFileName( outdir, getRootString(filenames[0]), "sum", "" ) );
+      std::cout << "Writing tiff to " << imgfname << " " << thr1 << ".tiff" << std::endl;
+      mt->writeImage(imgfname.c_str(), thr1);
+    }
+
+
+    //std::cout << "Calling delete..." << std::endl;
+    /* Info: Previously, 'delete mt' caused crash
+     (double calls of StopThread() in both destructors of
+     multiThreadedAnalogDetector and threadedAnalogDetector)
+     Now fixed! */
+    delete mt; // triggers cleanup of all threads and singlePhotonDetector instances (delete filter is obsolete)
+    delete decoder;
+    free(buff); // Free explicitly as it gets popped out of the Fifo at termination of while(readNextFrame)
+
+    return 0;
+}

slsDetectorCalibration/jungfrauExecutables/jungfrauRawDataProcess_filetxtH5_SC.cpp

@@ -0,0 +1,468 @@
+// SPDX-License-Identifier: LGPL-3.0-or-other
+// Copyright (C) 2021 Contributors to the SLS Detector Package
+// #include "sls/ansi.h"
+//#include <iostream>
+#undef CORR
+
+#define RAWDATA
+
+
+#include "jungfrauLGADStrixelsDataQuadH5.h"
+
+#include "multiThreadedCountingDetector.h"
+#include "singlePhotonDetector.h"
+
+#include <fstream>
+#include <csignal>
+#include <map>
+#include <memory>
+#include <stdio.h>
+#include <sys/stat.h>
+
+#include <ctime>
+#include <fmt/core.h>
+
+
+/*Dataset paths according to different beamlines*/
+std::string const data_datasetname_furka("/data/JF18T01V01/data");
+std::string const index_datasetname_furka("/data/JF18T01V01/frame_index");
+std::string const data_datasetname_xfelSCS("/INSTRUMENT/SCS_HRIXS_JUNGF/DET/JNGFR01:daqOutput/data/adc");
+std::string const index_datasetname_xfelSCS("/INSTRUMENT/SCS_HRIXS_JUNGF/DET/JNGFR01:daqOutput/data/frameNumber");
+
+std::string getRootString( std::string const& filepath ) {
+  size_t pos1;
+  if (filepath.find("/") == std::string::npos )
+    pos1 = 0;
+  else
+    pos1 = filepath.find_last_of("/")+1;
+  size_t pos2 = filepath.find_last_of(".");
+  //std::cout << "pos1 " << pos1 << " pos2 " << pos2 << " size " << filepath.length() << std::endl;
+  return filepath.substr( pos1, pos2-pos1 );
+}
+
+/* Create file name string
+ * \param dir directory
+ * \param fprefix fileprefix (without extension)
+ * \param fsuffix filesuffix (for output files, e.g. "ped")
+ * \param fext    file extension (e.g. "raw")
+ */
+std::string createFileName( std::string const& dir, std::string const& fprefix="run",
+			                std::string const& fsuffix="", std::string const& fext="raw", int const outfilecounter=-1 ) {
+    std::string return_string;
+    if (outfilecounter >= 0)
+        return_string = fmt::format("{:s}/{:s}_{:s}_f{:05d}", dir, fprefix, fsuffix, outfilecounter);
+    else if (fsuffix.length()!=0)
+        return_string = fmt::format("{:s}/{:s}_{:s}", dir, fprefix, fsuffix);
+    else
+        return_string = fmt::format("{:s}/{:s}", dir, fprefix);
+
+    if (fext.length()!=0)
+        return_string += "." + fext;
+
+    return return_string;
+}
+
+/* Adjusts number of threads to be a multiple of number of storage cells
+ * \param requestedThreads number of threads requested by the user
+ * \param nSC number of storage cells
+ */
+int adjustThreads(int requestedThreads, int nSC) {
+    if (nSC <= 0) {
+        std::cerr << "Error: Number of S values must be greater than zero!" << std::endl;
+        return requestedThreads; // Return the original value as a fallback
+    }
+
+    // Calculate the remainder
+    int remainder = requestedThreads % nSC;
+
+    // If remainder is non-zero, round up by adding the difference
+    int adjustedThreads = (remainder == 0) ? requestedThreads : requestedThreads + (nSC - remainder);
+
+    // Ensure at least `nSC` threads are used
+    if (adjustedThreads < nSC) {
+        adjustedThreads = nSC;
+    }
+
+    std::cout << "Adjusted thread count (rounded up): " << adjustedThreads << " (nearest multiple of " 
+              << nSC << ")" << std::endl;
+
+    return adjustedThreads;
+}
+
+// Signal handler for segmentation faults
+void signal_handler(int signum) {
+    std::cerr << "Caught signal " << signum << ": Segmentation fault (core dump)" << std::endl;
+    // Handle the error (e.g., clean up, abort, etc.)
+    exit(signum);  // Exit program with the signal code
+}
+
+
+int main(int argc, char *argv[]) {
+
+    if (argc < 4) {
+        std::cout
+            << "Usage is " << argv[0]
+            << " filestxt outdir [pedfile (h5)] "
+               " optional: [int dataset path; 0 means Furka, 1 means XFEL; overwrites default given in HDF5File.h] "
+               " [bool validate h5 rank] "
+               " [xmin xmax ymin ymax] [nframes] "
+	             " NOTE THAT THE DATA FILES HAVE TO BE IN THE RIGHT ORDER SO THAT PEDESTAL TRACKING WORKS! "
+            << std::endl;
+        return 1;
+    }
+
+    // Set up the signal handler for segmentation faults
+    signal(SIGSEGV, signal_handler);
+
+    int const fifosize = 100; //1000;
+    int const nthreads = 10;
+    int const csize = 3; // 3
+    int const nsigma = 5;
+    int const nped = 10000;
+
+    //int cf = 0;
+
+    std::string const txtfilename(argv[1]);
+    std::string const outdir(argv[2]);
+    std::string const pedfilename(argv[3]);
+
+    std::string datasetpath{};
+    std::string frameindexpath{};
+    if (argc > 4) {
+      switch (atoi(argv[4])) {
+      case 0:
+        datasetpath = data_datasetname_furka;
+        frameindexpath = index_datasetname_furka;
+        break;
+      case 1:
+        datasetpath = data_datasetname_xfelSCS;
+        frameindexpath = index_datasetname_xfelSCS;
+        break;
+      default:
+        break;
+      }
+    }
+
+    bool validate_rank=true;
+    if (argc > 5)
+      validate_rank = atoi(argv[5]);
+
+    //Get vector of filenames from input txt-file
+    std::vector<std::string> filenames{};
+
+    { //Safety scope for ifstream
+      ifstream inputs( txtfilename, std::ios::in );
+      if (inputs.is_open()) {
+	    std::cout << "Reading imput filenames from txt-file ..." << std::endl;
+	    std::string line{};
+	    while (!inputs.eof()) {
+	      std::getline(inputs, line);
+            if(line.find(".h5") != std::string::npos) {
+	            filenames.emplace_back(line);
+	            std::cout << line << std::endl;
+            }
+	    }
+	    inputs.close();
+	    std::cout << "---- Reached end of txt-file. ----" << std::endl;
+      } else
+	      std::cout << "Could not open " << txtfilename << std::endl;
+      if (filenames.size()>0) {
+	      std::cout << filenames.size() << " filenames found in " << txtfilename << std::endl;
+	      std::cout << "The files will be processed in the same order as found in the txt-file." << std::endl;
+      } else {
+	      std::cout << "No files found in txt-file!" << std::endl;
+	      return 1;
+      }
+    }
+
+    std::cout << "###############" << std::endl;
+
+    // Define decoder
+    std::cout << "Jungfrau strixel quad h5" << std::endl;
+    jungfrauLGADStrixelsDataQuadH5* decoder = new jungfrauLGADStrixelsDataQuadH5();
+    //auto decoder = std::make_unique<jungfrauLGADStrixelsDataQuadH5>();
+    int nx = 1024 / 3, ny = 512 * 3;
+
+    //Cluster finder ROI
+    int xmin = 0, xmax = nx-1, ymin = 0, ymax = ny-1;
+    if (argc > 9) {
+      xmin = atoi(argv[6]);
+      xmax = atoi(argv[7]);
+      ymin = atoi(argv[8]);
+      ymax = atoi(argv[9]);
+    }
+    std::cout << "Cluster finder ROI: [" << xmin << ", " << xmax << ", " << ymin << ", " << ymax << "]"
+              << std::endl;
+
+    decoder->getDetectorSize(nx, ny);
+    std::cout << "Detector size is " << nx << " " << ny << std::endl;
+
+    std::time_t end_time;
+
+    std::cout << "Output directory is " << outdir << std::endl;
+    if (pedfilename.length()!=0)
+        std::cout << "Pedestal file is " << pedfilename << std::endl;
+
+    //uint32_t nnx, nny;
+
+    singlePhotonDetector* filter =
+      new singlePhotonDetector(decoder, 3, nsigma, 1, NULL, nped, 200, -1, -1, NULL, NULL);
+    //auto filter = std::make_unique<singlePhotonDetector>(decoder.get(), 3, nsigma, 1, nullptr, nped, 200, -1, -1, nullptr, nullptr);
+
+    filter->setROI(xmin, xmax, ymin, ymax);
+
+    std::time(&end_time);
+    std::cout << std::ctime(&end_time) << std::endl;
+
+    // Validate number of threads for number of storage cells (if applicable)
+    int nThreads = nthreads;
+    int nSC = 1;
+    // Determine the dimensions of the dataset from the first datafile
+    auto firstfileh5 = std::make_unique<HDF5File>();
+    firstfileh5->SetFrameIndexPath(frameindexpath);
+    firstfileh5->SetImageDataPath(datasetpath);
+    //std::cout << "Debug: Attempting to open file " << filenames[0].c_str() << std::endl;
+    if ( firstfileh5->OpenResources(filenames[0].c_str(), validate_rank) ) {
+
+        // Validate number of threads
+        if( firstfileh5->GetRank() == 4 ) {
+            auto h5dims = firstfileh5->GetDatasetDimensions();
+            nSC = h5dims[1];
+            nThreads = adjustThreads(nthreads,nSC);
+        }
+
+        firstfileh5->CloseResources();
+    } else {
+        std::cerr << "Error: Could not open data file " << filenames[0]
+		          << " for validating rank " << std::endl;
+    }
+
+    multiThreadedCountingDetector* mt =
+        new multiThreadedCountingDetector(filter, nThreads, fifosize, nSC);
+    //auto mt = std::make_unique<multiThreadedCountingDetector>(filter.get(), nthreads, fifosize);
+    mt->setClusterSize(csize, csize);
+    mt->newDataSet(); //Initialize new dataset for each thread
+    mt->setDetectorMode(ePhotonCounting);
+
+    char* buff;
+
+    mt->StartThreads();
+    mt->popFree(buff); // Get the first pointer to write image to
+
+    size_t ifr = 0; //frame counter of while loop
+    int framenumber = 0; //framenumber as read from file (detector)
+    std::vector<hsize_t> h5offset; //hyperslab offset internal to HDF5File::ReadImage
+    hsize_t h5rank;
+
+    if (pedfilename.length()>1) {
+
+        std::string froot = getRootString(pedfilename);
+
+        std::cout << "PEDESTAL " << pedfilename << std::endl;
+	    std::time(&end_time);
+	    std::cout << "aaa " << std::ctime(&end_time) << std::endl;
+
+	    mt->setFrameMode(ePedestal);
+
+	    //HDF5File pedefile;
+	    auto pedefile = std::make_unique<HDF5File>();
+        pedefile->SetFrameIndexPath(frameindexpath);
+        pedefile->SetImageDataPath(datasetpath);
+	    //      //open file
+	    if ( pedefile->OpenResources(pedfilename.c_str(),validate_rank) ) {
+
+            // Initialize offset vector to 0
+            h5rank = pedefile->GetRank();
+	        h5offset.resize(h5rank-2, 0);
+	        framenumber = 0;
+
+	        while ( decoder->readNextFrame(*pedefile, framenumber, h5offset, buff) ) {
+
+	            if ((ifr + 1) % 100 == 0) {
+	                std::cout
+		                << " ****"
+		                << decoder->getValue(buff, 20, 20); // << std::endl;
+	            }
+
+                int storageCell = 0;
+                hsize_t n_storageCells = 1;
+                if (h5rank == 4) {
+                    storageCell = h5offset[1];
+                    n_storageCells = pedefile->GetDatasetDimensions()[1];
+                }
+
+                // push buff into fifoData for a thread corresponding to the active storage cell
+	            mt->pushData(buff, storageCell);
+                // increment (round-robin) the internal thread counter for that storage cell
+	            mt->nextThread(storageCell);
+                // get a free memory address from fifoFree of the active storage cell for the next read operation
+	            mt->popFree(buff, storageCell);
+                /* NOTE: the buff that was popped free from the current thread, will be (likely) pushed into
+                 * the fifoData of a different thread in the next iteration of the loop! */
+
+	            ++ifr;
+	            if (ifr % 100 == 0) {
+	                std::cout << " ****" << ifr << " " << framenumber << " " << h5offset[0];
+                    if (n_storageCells>1)
+			            std::cout << " sc " << storageCell;
+                    std::cout << "\n";
+	            } // else
+	      
+	            if (ifr >= 1000*n_storageCells)
+	                break;
+	          //framenumber = 0;
+	        }
+	    
+	        pedefile->CloseResources();
+	        while (mt->isBusy()) {
+	          ;
+	        }
+
+	        std::cout << "Writing pedestal to " << getRootString(pedfilename) << "_ped_SCxx.tiff" << std::endl;
+	        auto imgfname = createFileName( outdir, getRootString(pedfilename), "ped", "" );
+	        mt->writePedestal(imgfname.c_str());
+            std::cout << "Writing pedestal rms to " << getRootString(pedfilename) << "_rms_SCxx.tiff" << std::endl;
+	        imgfname = createFileName( outdir, getRootString(pedfilename), "rms", "");
+	        mt->writePedestalRMS(imgfname.c_str());
+
+	      } else {
+	        std::cerr << "Error: Could not open pedestal file " << pedfilename
+		                << " for reading " << std::endl;
+          }
+
+      std::time(&end_time);
+      std::cout << std::ctime(&end_time) << std::endl;
+    }
+    
+    ifr = 0;
+    //int ioutfile = 0;
+
+    mt->setFrameMode(eFrame);
+
+    std::vector<FILE*> of(nSC, nullptr);
+    
+    for (unsigned int ifile = 0; ifile != filenames.size(); ++ifile) {
+        std::cout << "DATA " << filenames[ifile] << " " << std::endl;
+        std::time(&end_time);
+        std::cout << std::ctime(&end_time) << std::endl;
+
+        //HDF5File fileh5;
+        auto fileh5 = std::make_unique<HDF5File>();
+        fileh5->SetFrameIndexPath(frameindexpath);
+        fileh5->SetImageDataPath(datasetpath);
+
+        // Open HDF5 file
+        if ( fileh5->OpenResources(filenames[ifile].c_str(), validate_rank) ) {
+
+            std::vector<std::string> cfnames(nSC);
+            for ( int s = 0; s < nSC; ++s ) {
+                std::string fsuffix = "SC" + std::to_string(s);
+                cfnames[s] = createFileName( outdir, getRootString(filenames[ifile]), fsuffix, "clust" );
+            }
+	      
+            //Open output files and set file pointers according to storage cells
+            for ( size_t f = 0; f < of.size(); ++f ) {
+                if (!of[f]) {
+                    of[f] = fopen(cfnames[f].c_str(), "w");
+                    if (of[f])
+                    {
+                        if (mt) {
+		                    mt->setFilePointer(of[f],f); // assumes f == sc
+		                    std::cout << "File pointer set for storage cell " << f << std::endl;
+	                    } else {
+		                    std::cerr << "Error: mt is null." << std::endl;
+		                    return 1;
+	                    }
+                    } else {
+	                    std::cerr << "Error: could not open " << cfnames[f]
+			                      << " for writing " << std::endl;
+	                    mt->setFilePointer(nullptr,f);
+	                    return 1;
+                
+                    }
+                }
+            }
+	      
+	
+	        // Read frames
+	        framenumber = 0;
+	        std::fill(h5offset.begin(), h5offset.end(), 0);
+	        ifr = 0;
+	        
+	        while ( decoder->readNextFrame(*fileh5, framenumber, h5offset, buff) ) {
+	        
+	            if ((ifr + 1) % 1000 == 0) {
+	                std::cout << " ****"
+		                      << decoder->getValue(buff, 20, 20); // << std::endl;
+	            }
+
+                int storageCell = 0;
+                hsize_t n_storageCells = 1;
+                if (h5rank == 4) {
+                    storageCell = h5offset[1];
+                    n_storageCells = fileh5->GetDatasetDimensions()[1];
+                }
+
+	            // push buff into fifoData for a thread corresponding to the active storage cell
+	            mt->pushData(buff, storageCell);
+                // increment (round-robin) the internal thread counter for that storage cell
+	            mt->nextThread(storageCell);
+                // get a free memory address from fifoFree of the active storage cell for the next read operation
+	            mt->popFree(buff, storageCell); /* In the last execution of the loop,
+                                                * this leaves buff outside of the Fifo!
+                                                * Free explicitely at the end! */
+	  
+	            ++ifr;
+	            if (ifr % 1000 == 0) {
+	                std::cout << " " << ifr << " " << framenumber << " " << h5offset[0];
+                    if (n_storageCells>1) std::cout << " sc " << storageCell;
+                    std::cout << "\n";
+                }
+
+	            //framenumber = 0;
+	        }
+
+	        //std::cout << "aa --" << std::endl;
+	        fileh5->CloseResources();
+
+	        //std::cout << "bb --" << std::endl;
+	        while (mt->isBusy()) {
+	            ;
+	        }
+
+	        //std::cout << "cc --" << std::endl;
+	      
+	        auto imgfname = createFileName( outdir, getRootString(filenames[ifile]), "", "" );
+	        std::cout << "Writing tiff to " << imgfname << "_SCxx.tiff" << std::endl;
+	        mt->writeImage(imgfname.c_str());
+	        mt->clearImage();
+
+            // Close output files
+            for ( size_t f = 0; f < of.size(); ++f ) {
+	            if (of[f]) {
+	                fclose(of[f]);
+	                mt->setFilePointer(nullptr,f);
+                    of[f] = nullptr;
+	            }
+            }
+	      
+	        std::time(&end_time);
+	        std::cout << std::ctime(&end_time) << std::endl;
+        } else {
+	        std::cerr << "Error: Could not open " << filenames[ifile] << " for reading "
+		              << std::endl;
+        }
+    }
+    
+
+    //std::cout << "Calling delete..." << std::endl;
+    /* Info: Previously, 'delete mt' caused crash
+     (double calls of StopThread() in both destructors of
+     multiThreadedAnalogDetector and threadedAnalogDetector)
+     Now fixed! */
+    delete mt; // triggers cleanup of all threads and singlePhotonDetector instances (delete filter is obsolete)
+    delete decoder;
+    free(buff); // Free explicitly as it gets popped out of the Fifo at termination of while(readNextFrame)
+
+    return 0;
+}

slsDetectorCalibration/moenchExecutables/moenchZmqProcess.cpp

@@ -336,7 +336,7 @@ int main(int argc, char *argv[]) {
 
     string fname;
     //	int length;
-    int *detimage = NULL;
+    int *detimage = nullptr;
     int nnx, nny, nnsx, nnsy;
     // uint32_t imageSize = 0, nPixelsX = 0, nPixelsY = 0,
     // uint32_t  dynamicRange = 0;

slsDetectorCalibration/multiThreadedAnalogDetector.h

@@ -23,7 +23,9 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <unistd.h>
-//#include <mutex>
+#include <mutex>
+#include <algorithm>
+#include<unordered_map>
 
 using namespace std;
 
@@ -47,12 +49,20 @@ class threadedAnalogDetector {
 
             if (mm) {
                 // memset(mm,0, det->getDataSize());
+                /*
+                if (i == 0) { // debug
+                    first_mm = mm;
+                }
+                */
+                
                 fifoFree->push(mm);
+                //std::cout << "Allocated memory at: " << static_cast<void*>(mm) << " (fifoslot " << i << ")" << std::endl;
             } else
                 break;
         }
+
         if (i < fs)
-            cout << "Could allocate only " << i << " frames";
+            std::cout << "Could allocate only " << i << " frames";
 
         busy = 0;
         stop = 1;
@@ -102,24 +112,50 @@ class threadedAnalogDetector {
     };
 
     virtual ~threadedAnalogDetector() {
+
+        // std::cout << "#### Debug: Destructing threadedAnalogDetector! ####" << std::endl;
+
         StopThread();
-        delete fifoFree;
-        delete fifoData;
+
+        if (fifoFree) { delete fifoFree; fifoFree = nullptr; }
+        if (fifoData) { delete fifoData; fifoData = nullptr; }
+        if (det) { 
+            delete det; // Call destructor for singlePhotonDetector
+            det = nullptr; 
+        } 
     }
 
     /** Returns true if the thread was successfully started, false if there was
      * an error starting the thread */
     virtual bool StartThread() {
         stop = 0;
-        cout << "Detector number " << det->getId() << endl;
-        cout << "common mode is " << det->getCommonModeSubtraction() << endl;
-        cout << "ghos summation is " << det->getGhostSummation() << endl;
+        std::cout << "Detector number " << det->getId() << std::endl;
+        std::cout << "common mode is " << det->getCommonModeSubtraction() << std::endl;
+        std::cout << "ghos summation is " << det->getGhostSummation() << std::endl;
+
         return (pthread_create(&_thread, NULL, processData, this) == 0);
     }
 
     virtual void StopThread() {
         stop = 1;
-        (void)pthread_join(_thread, NULL);
+        //std::cout << "Attempting to stop thread..." << std::endl;
+
+        // Free all remaining allocated memory in fifoFree
+        char *mm = nullptr;
+        while (fifoFree->pop(mm, true)) {  // Use no_block to avoid waiting
+            //std::cout << "fifo Free: Freeing memory at: " << static_cast<void*>(mm) << std::endl;
+            free(mm);  // Free the allocated memory
+        }
+
+        if (_thread) {
+            //(void)pthread_join(_thread, NULL);
+            //std::cout << "Calling pthread_join for thread: " << det->getId() << std::endl;
+            pthread_join(_thread, NULL);
+            _thread = 0;
+            std::cout << "Thread " << det->getId() << " stopped and joined." << std::endl;
+        } else {
+            std::cout << "No thread to join." << std::endl;
+        }
     }
 
     virtual bool pushData(char *&ptr) { return fifoData->push(ptr); }
@@ -266,6 +302,8 @@ class threadedAnalogDetector {
     char *data;
     int *ff;
 
+    //char* first_mm = nullptr; // For debug; to track first allocated block
+
     static void *processData(void *ptr) {
         threadedAnalogDetector *This = ((threadedAnalogDetector *)ptr);
         return This->processData();
@@ -278,55 +316,118 @@ class threadedAnalogDetector {
                 usleep(100);
                 if (fifoData->isEmpty()) {
                     busy = 0;
-                } else
+                } else {
                     busy = 1;
-            } else
+                }
+            } else {
                 busy = 1;
+            }
 
             if (busy == 1) {
+                // Check stop flag before making a blocking call
+                //if (stop) {
+                //    break;
+                //}
+
+                // Blocking call
                 fifoData->pop(data); // blocking!
+
+                // Process data if not stopping
+                //if (!stop) {
                 det->processData(data);
                 fifoFree->push(data);
+                    
+                //}
                 // busy=0;
             }
         }
+
         return NULL;
     }
 };
 
 class multiThreadedAnalogDetector {
   public:
-    multiThreadedAnalogDetector(analogDetector<uint16_t> *d, int n,
-                                int fs = 1000)
-        : stop(0), nThreads(n), ithread(0) {
+    multiThreadedAnalogDetector(analogDetector<uint16_t> *d, int num_threads,
+                                int fs = 1000, int num_sc = 1)
+        : stop(0), nThreads(num_threads), nSC(num_sc) {
+
+        /*
         dd[0] = d;
         if (nThreads == 1)
             dd[0]->setId(100);
         else
             dd[0]->setId(0);
-        for (int i = 1; i < nThreads; i++) {
-            dd[i] = d->Clone();
-            dd[i]->setId(i);
+        */
+
+        // Create separate detectorObjects for each SC (each owns its mutex)
+        std::vector< analogDetector<uint16_t>* > sc_detectors(nSC, nullptr);
+        sc_detectors[0] = d; // First storage cell uses the given detector
+        // std::cout << "#### Debug: Copied analogDetector object for storage cell 0! ####" << std::endl;
+
+        for (int sc = 1; sc < nSC; ++sc) {
+            sc_detectors[sc] = d->Copy(); // Ensure unique mutex for each SC
+            // std::cout << "#### Debug: Copied analogDetector object for storage cell " << sc << "! ####" << std::endl;
         }
 
+        // Distribute threads among storage cells
+        int threads_per_sc = nThreads / nSC;
+        int remaining_threads = nThreads % nSC; // additional safety measure if nThreads is not divisible by nSC
+
+        sc_to_threads.clear();
+        for (int s = 0, thread_idx = 0; s < nSC; ++s) {
+            // Remaining threads (if any) are assigned to the first storage cells
+            int current_sc_threads = threads_per_sc + (s < remaining_threads ? 1 : 0);
+
+            for (int t = 0; t < current_sc_threads; ++t, ++thread_idx) {
+
+                if (t == 0) {
+                    dd[thread_idx] = sc_detectors[s]; // First thread gets main SC detector
+                } else {
+                    dd[thread_idx] = sc_detectors[s]->Clone(); // Other threads get clones
+                }
+
+                std::cout << "Assigned thread " << thread_idx << " to storage cell " << s << std::endl;
+
+                dd[thread_idx]->setId(thread_idx);
+                
+                // Store which threads belong to which SC
+                sc_to_threads[s].push_back(thread_idx);
+            }
+        }
+
+        if (nSC == 1 && nThreads == 1) {
+            dd[0]->setId(100);
+        }
+
+        // Initialize threadedAnalogDetector objects
         for (int i = 0; i < nThreads; i++) {
             cout << "**" << i << endl;
             dets[i] = new threadedAnalogDetector(dd[i], fs);
         }
 
-        image = NULL;
+        // Set all thread counters to zero for each storage cell
+        thread_counters_by_sc.resize(nSC,0);
+
+        image = nullptr;
         ff = NULL;
         ped = NULL;
-        cout << "Ithread is " << ithread << endl;
+        //std::cout << "Ithread is " << ithread << std::endl;
     }
 
     virtual ~multiThreadedAnalogDetector() {
-        StopThreads();
-        for (int i = 0; i < nThreads; i++)
-            delete dets[i];
-        /* for (int i=1; i<nThreads; i++)  */
-        /*   delete dd[i]; */
-        // delete [] image;
+        // std::cout << "#### Debug: Destructing multiThreadedAnalogDetector! ####" << std::endl;
+        //StopThreads(); // Superfluous, leads to double delete
+        
+        /* Reverse loop for destruction.
+         * Deletes clones first, then root object, which owns the mutex
+         * (ensure shared mutex is deleted last).
+         * Optional solution: reference counting (safer but more complex) */
+        for (int i = nThreads - 1; i >= 0; --i) {
+            delete dets[i]; //StopThread() called by each ~threadedAnalogDetector()
+            dets[i] = nullptr;
+        }
+        
     }
 
     virtual int setFrameMode(int fm) {
@@ -359,36 +460,52 @@ class multiThreadedAnalogDetector {
             dets[i]->newDataSet();
     };
 
-    virtual int *getImage(int &nnx, int &nny, int &ns, int &nsy) {
-        int *img;
+    // Storage cell sensitive
+    virtual int *getImage(int &nnx, int &nny, int &ns, int &nsy, int sc = 0) {
+        //int *img;
         // int nnx, nny, ns;
         // int nnx, nny, ns;
         int nn = dets[0]->getImageSize(nnx, nny, ns, nsy);
-        if (image) {
-            delete[] image;
-            image = NULL;
+
+        if (sc_images[sc]) {
+            delete[] sc_images[sc];
+            sc_images[sc] = nullptr;
         }
-        image = new int[nn];
+
+        // Allocate memory for image and zero-initialize
+        sc_images[sc] = new int[nn]();
         // int nn=dets[0]->getImageSize(nnx, nny, ns);
         // for (i=0; i<nn; i++) image[i]=0;
 
-        for (int ii = 0; ii < nThreads; ii++) {
-            // cout << ii << " " << nn << " " << nnx << " " << nny << " " << ns
-            // << endl;
-            img = dets[ii]->getImage();
+        // Get the threads assigned to this storage cell
+        auto const& assigned_threads = sc_to_threads[sc];
+
+        // Only iterate over threads assigned to this storage cell
+        for (int thread_id : assigned_threads) {
+            
+            int* tmp_img = dets[thread_id]->getImage();
+            if (!tmp_img) continue; // Skip if null
+
+            /* std::cout << "## Thread " << ii
+                      << " # image size " << nn
+                      << " # nnx " << nnx
+                      << " # nny " << nny
+                      << " # ns " << ns; */
+
+            // Sum images across threads
             for (int i = 0; i < nn; i++) {
-                if (ii == 0)
-                    //	  if (img[i]>0)
-                    image[i] = img[i];
-                // else
-                //  image[i]=0;
-                else // if (img[i]>0)
-                    image[i] += img[i];
+
+                /* std::cout << " # pixel " << i
+                          << " # value " << tmp_img[i]
+                          << " ## " << std::endl; */
+
+                
+                sc_images[sc][i] += tmp_img[i];
                 // if (img[i])	  cout << "det " << ii << " pix " << i << " val
                 // " <<  img[i] << " " << image[i] << endl;
             }
         }
-        return image;
+        return sc_images[sc];
     }
 
     virtual void clearImage() {
@@ -398,7 +515,7 @@ class multiThreadedAnalogDetector {
         }
     }
 
-    virtual void *writeImage(const char *imgname, double t = 1) {
+    virtual void *writeImage(char const* base_imgname, double t = 1) {
         /* #ifdef SAVE_ALL   */
         /*      for (int ii=0; ii<nThreads; ii++) { */
         /*        char tit[10000];cout << "m" <<endl; */
@@ -407,11 +524,30 @@ class multiThreadedAnalogDetector {
         /*      } */
         /* #endif */
         int nnx, nny, ns, nsy;
-        getImage(nnx, nny, ns, nsy);
-        // int nnx, nny, ns;
         int nn = dets[0]->getImageSize(nnx, nny, ns, nsy);
-        float *gm = new float[nn];
-        if (gm) {
+
+        // Allocate teporary float buffer and zero-initialize
+        std::vector<float> gm(nn);
+
+        // Lambda for pixel conversion
+        auto convert_pixel = [t](int pixel) -> float {
+            return (t > 0) ? static_cast<float>(std::max(0, pixel)) / static_cast<float>(t) : pixel;
+        }; // t ... threshold
+
+        // Loop over each storage cell
+        for (auto const& [sc, _] : sc_to_threads) { // structured bindings [sc, _] only available with -std=c++17
+            std::string imgname(base_imgname);
+            if (nSC > 1) imgname += "_SC" + std::to_string(sc);
+            imgname += ".tiff";
+
+            //Retrieve the image for this storage cell
+            int *image = getImage(nnx, nny, ns, nsy, sc);
+            if (!image) continue; // Skip if null
+
+            // Convert image data to float
+            std::transform(image, image + nn, gm.begin(), convert_pixel);
+
+            /* old loop implementing same logic as convert_pixel
             for (int ix = 0; ix < nn; ix++) {
                 if (t) {
                     if (image[ix] < 0)
@@ -424,11 +560,18 @@ class multiThreadedAnalogDetector {
                 // if (image[ix]>0 && ix/nnx<350) cout << ix/nnx << " " <<
                 // ix%nnx << " " << image[ix]<< " " << gm[ix] << endl;
             }
-            // cout << "image " << nnx << " " << nny << endl;
-            WriteToTiff(gm, imgname, nnx, nny);
-            delete[] gm;
-        } else
-            cout << "Could not allocate float image " << endl;
+            */
+
+            WriteToTiff(gm.data(), imgname.c_str(), nnx, nny);
+
+            // Clean up memory for this storage cell
+            if (sc_images[sc]) {
+                delete[] sc_images[sc];
+                sc_images[sc] = nullptr;
+            }
+            
+        }
+        
         return NULL;
     }
 
@@ -439,6 +582,7 @@ class multiThreadedAnalogDetector {
     }
 
     virtual void StopThreads() {
+        std::cout << "Stopping all threads ..." << std::endl;
         for (int i = 0; i < nThreads; i++)
             dets[i]->StopThread();
     }
@@ -453,76 +597,133 @@ class multiThreadedAnalogDetector {
         return ret;
     }
 
-    virtual bool pushData(char *&ptr) { return dets[ithread]->pushData(ptr); }
+/*
+    virtual std::vector<int> getThreadsForSc(int sc) {
+        return sc_to_threads[sc];
+    }
+*/
+    virtual bool pushData(char *&ptr, int sc=0) {
+        //Additional logic implemented to accommodate storage cells
 
-    virtual bool popFree(char *&ptr) {
-        //  cout << ithread << endl;
-        return dets[ithread]->popFree(ptr);
+        std::unique_lock<std::mutex> lock(map_mutex);
+
+        // Get assigned threads for this storage cell
+        auto& assigned_threads = sc_to_threads[sc];
+        auto& counter = thread_counters_by_sc[sc];
+        
+        // Distribute workload among threads using round-robin
+        int selected_thread = assigned_threads[counter % assigned_threads.size()];
+
+        return dets[selected_thread]->pushData(ptr); 
     }
 
-    virtual int nextThread() {
-        ithread++;
-        if (ithread == nThreads)
-            ithread = 0;
-        return ithread;
+    virtual bool popFree(char *&ptr, int sc=0) {
+        //Additional logic implemented to accommodate storage cells
+
+        std::unique_lock<std::mutex> lock(map_mutex);
+
+        // Get assigned threads for this storage cell
+        auto& assigned_threads = sc_to_threads[sc];
+        auto& counter = thread_counters_by_sc[sc];
+        
+        // Distribute workload among threads using round-robin
+        int selected_thread = assigned_threads[counter % assigned_threads.size()];
+
+        return dets[selected_thread]->popFree(ptr);
     }
 
-    virtual double *getPedestal() {
+    virtual int nextThread(int sc=0) {
+        //Additional logic implemented to accommodate storage cells
+
+        auto& counter = thread_counters_by_sc[sc];
+        //counter++;
+        if (++counter == nThreads/nSC)
+            counter = 0;
+        return counter;
+    }
+
+    // Storage cell sensitive
+    virtual double *getPedestal(int sc = 0) {
         int nx, ny;
         dets[0]->getDetectorSize(nx, ny);
-        if (ped)
-            delete[] ped;
-        ped = new double[nx * ny];
+
+        if (sc_pedestals.count(sc) && sc_pedestals[sc]) {
+            delete[] sc_pedestals[sc];
+            sc_pedestals[sc] = nullptr;
+        }
+
+        // allocate memory and initialize all values to zero
+        sc_pedestals[sc] = new double[nx * ny](); // parentheses initialize elements to zero
+        //std::fill(sc_pedestals[sc], sc_pedestals[sc] + (nx * ny), 0.0); // explicit zero initialization
         double *p0 = new double[nx * ny];
 
-        for (int i = 0; i < nThreads; i++) {
+        // Get the threads assigned to this storage cell
+        auto const& assigned_threads = sc_to_threads[sc];
+        int num_threads = assigned_threads.size();
+
+        // Only iterate over threads assigned to this storage cell
+        for ( int thread_id : assigned_threads ) {
             // inte=(slsInterpolation*)dets[i]->getInterpolation(nb,emi,ema);
             //   cout << i << endl;
-            p0 = dets[i]->getPedestal(p0);
-            if (p0) {
-                if (i == 0) {
+            //p0 = dets[thread_id]->getPedestal(p0);
+            dets[thread_id]->getPedestal(p0);
 
+            if (p0) { /*
+                if (i == 0) {
+                    // If first thread, initialize ped with first thread's values
                     for (int ib = 0; ib < nx * ny; ib++) {
                         ped[ib] = p0[ib] / ((double)nThreads);
                         //  cout << p0[ib] << " ";
                     }
                 } else {
-                    for (int ib = 0; ib < nx * ny; ib++) {
-                        ped[ib] += p0[ib] / ((double)nThreads);
-                        //  cout << p0[ib] << " ";
-                    }
+                */
+                    // For subsequent threads, accumulate pedestal values
+                // if ( i == 0 ) becomes superfluous if we zero-initialize earlier
+                for (int ib = 0; ib < nx * ny; ib++) {
+                    sc_pedestals[sc][ib] += p0[ib] / ((double)num_threads);
+                    //  cout << p0[ib] << " ";
                 }
+                //}
             }
         }
         delete[] p0;
-        return ped;
+        return sc_pedestals[sc];
     };
 
-    virtual double *getPedestalRMS() {
+    // Storage cell sensitive
+    virtual double *getPedestalRMS(int sc = 0) {
         int nx, ny;
         dets[0]->getDetectorSize(nx, ny);
-        // if (ped) delete [] ped;
-        double *rms = new double[nx * ny];
+        
+        if (sc_pedestals_rms.count(sc) && sc_pedestals_rms[sc]) {
+            delete[] sc_pedestals_rms[sc];
+            sc_pedestals_rms[sc] = nullptr;
+        }
+
+        // allocate memory and initialize all values to zero
+        sc_pedestals_rms[sc] = new double[nx * ny](); // Zero-initialize
+        //std::fill(sc_pedestals_rms[sc], sc_pedestals_rms[sc] + (nx * ny), 0.0); // explicit zero initialization
+        //double *rms = sc_pedestals_rms[sc];
         double *p0 = new double[nx * ny];
 
-        for (int i = 0; i < nThreads; i++) {
+        // Get the threads assigned to this storage cell
+        auto const& assigned_threads = sc_to_threads[sc];
+        int num_threads = assigned_threads.size();
+
+        // Only iterate over threads assigned to this storage cell
+        for (int thread_id : assigned_threads) {
             // inte=(slsInterpolation*)dets[i]->getInterpolation(nb,emi,ema);
             //   cout << i << endl;
-            p0 = dets[i]->getPedestalRMS(p0);
-            if (p0) {
-                if (i == 0) {
+            //p0 = dets[thread_id]->getPedestalRMS(p0);
+            dets[thread_id]->getPedestalRMS(p0);
 
-                    for (int ib = 0; ib < nx * ny; ib++) {
-                        rms[ib] = p0[ib] * p0[ib] / ((double)nThreads);
+            if (p0) {
+                for (int ib = 0; ib < nx * ny; ib++) {
+                    sc_pedestals_rms[sc][ib] += (p0[ib] * p0[ib]) / ((double)num_threads);
                         //  cout << p0[ib] << " ";
-                    }
-                } else {
-                    for (int ib = 0; ib < nx * ny; ib++) {
-                        rms[ib] += p0[ib] * p0[ib] / ((double)nThreads);
-                        //  cout << p0[ib] << " ";
-                    }
                 }
             }
+            
         }
         delete[] p0;
 
@@ -533,64 +734,117 @@ class multiThreadedAnalogDetector {
         /* 	rms[ib]=0; */
         /* } */
 
-        return rms;
+        return sc_pedestals_rms[sc];
     };
 
-    virtual double *setPedestal(double *h = NULL) {
+    /**
+     * Sets pedestal for given storage cell
+     * \param h pedestal
+     * \param sc storage cell
+     * \returns NULL
+     */
+    virtual double *setPedestal(double *h = NULL, int sc = 0) {
         // int nb=0;
 
         int nx, ny;
         dets[0]->getDetectorSize(nx, ny);
         if (h == NULL)
             h = ped;
-        for (int i = 0; i < nThreads; i++) {
+        for (auto i : sc_to_threads[sc]) {
             dets[i]->setPedestal(h);
         }
 
         return NULL;
     };
 
-    virtual void *writePedestal(const char *imgname) {
+    // Storage cell sensitive
+    virtual void *writePedestal(char const* base_imgname) {
 
         int nx, ny;
         dets[0]->getDetectorSize(nx, ny);
 
-        getPedestal();
-        float *gm = new float[nx * ny];
-        if (gm) {
+        //float *gm = new float[nx * ny];
+
+        // Loop over each storage cell
+        for ( auto const& entry : sc_to_threads ) {
+            int sc = entry.first;
+            std::string imgname = std::string(base_imgname);
+            if (nSC>1)
+                imgname += "_SC" + std::to_string(sc);
+            imgname += ".tiff";
+
+            getPedestal(sc); // Compute pedestal for this storage cell
+            std::vector<float> gm(nx * ny);
+
+            // Copy pedestal data into the float array
+            /*
             for (int ix = 0; ix < nx * ny; ix++) {
-                gm[ix] = ped[ix];
+                gm[ix] = sc_pedestals[sc][ix];
             }
-            WriteToTiff(gm, imgname, nx, ny);
-            delete[] gm;
-        } else
-            cout << "Could not allocate float image " << endl;
+            */
+            std::copy(sc_pedestals[sc], sc_pedestals[sc] + (nx * ny), gm.data());
+            WriteToTiff(gm.data(), imgname.c_str(), nx, ny);
+            
+            // Clean up memory
+            if(sc_pedestals[sc]) {
+                delete[] sc_pedestals[sc];
+                sc_pedestals[sc] = nullptr;
+            }
+        }
 
         return NULL;
     };
 
-    virtual void *writePedestalRMS(const char *imgname) {
+    // Storage cell sensitive
+    virtual void *writePedestalRMS(char const* base_imgname) {
 
         int nx, ny;
         dets[0]->getDetectorSize(nx, ny);
 
-        double *rms = getPedestalRMS();
-        float *gm = new float[nx * ny];
-        if (gm) {
+        //float *gm = new float[nx * ny];
+
+        // Loop over each stoarge cell
+        for ( auto const& entry : sc_to_threads ) {
+            int sc = entry.first;
+            std::string imgname = std::string(base_imgname);
+            if (nSC>1)
+                imgname += "_SC" + std::to_string(sc);
+            imgname += ".tiff";
+
+            double *rms = getPedestalRMS(sc); // Compute pedestal RMS for this storage cell
+            std::vector<float> gm(nx * ny);
+
+            // Copy rms data into the float array
+            /*
             for (int ix = 0; ix < nx * ny; ix++) {
-                gm[ix] = rms[ix];
+                    gm[ix] = rms[ix];
             }
-            WriteToTiff(gm, imgname, nx, ny);
-            delete[] gm;
-            delete[] rms;
-        } else
-            cout << "Could not allocate float image " << endl;
+            */
+            std::copy(rms, rms + (nx * ny), gm.data());
+            WriteToTiff(gm.data(), imgname.c_str(), nx, ny);
+
+            // Clean up memory
+            //delete[] rms; //This would cause double-free since the pointer is already handled by sc_pedestals_rms[sc]
+            if(sc_pedestals_rms[sc]) {
+                delete[] sc_pedestals_rms[sc]; 
+                sc_pedestals_rms[sc] = nullptr;
+            }  
+        }
 
         return NULL;
     };
 
+    /** 
+     * Reads pedestal and sets it for given storage cell
+     * \param imgname name of pedestal file
+     * \param nb obsolete
+     * \param emin obsolete
+     * \param emax obsolete
+     * \param sc storage cell
+     * \returns setPedestal(NULL, sc)
+     * */
     virtual void *readPedestal(const char *imgname, int nb = -1,
-                               double emin = 1, double emax = 0) {
+                               double emin = 1, double emax = 0, int sc = 0) {
 
         int nx, ny;
         dets[0]->getDetectorSize(nx, ny);
@@ -610,36 +864,68 @@ class multiThreadedAnalogDetector {
         }
         delete[] gm;
 
-        return setPedestal();
+        return setPedestal(NULL, sc);
     };
 
-    /** sets file pointer where to write the clusters to
+    /** Sets file pointer where to write the clusters to
         \param f file pointer
-        \returns current file pointer
+        \param sc storage cell index
+        \returns current file pointer, or nullptr if invalid
     */
-    virtual FILE *setFilePointer(FILE *f) {
-        for (int i = 0; i < nThreads; i++) {
-            dets[i]->setFilePointer(f);
-            // dets[i]->setMutex(&fmutex);
+    virtual FILE *setFilePointer(FILE *f, int sc = 0) {
+        // Check if the storage cell exists
+        if (sc_to_threads.find(sc) == sc_to_threads.end() || sc_to_threads[sc].empty()) {
+            std::cerr << "Error: Invalid storage cell index " << sc << std::endl;
+            return nullptr;
         }
-        return dets[0]->getFilePointer();
+
+        // Assign file pointer to all threads belonging to this storage cell
+        for (auto i : sc_to_threads[sc]) {
+            if(dets[i]) {
+                dets[i]->setFilePointer(f);
+                // dets[i]->setMutex(&fmutex);
+            } else {
+                std::cerr << "Warning: dets[" << i << "] is null, skipping file pointer set." << std::endl;
+            }
+        }
+        // Return file pointer of the first thread in this storage cell
+        return dets[sc_to_threads[sc][0]] ? dets[sc_to_threads[sc][0]]->getFilePointer() : nullptr;
     };
 
-    /** gets file pointer where to write the clusters to
-        \returns current file pointer
+    /** Gets file pointer where to write the clusters to
+        \param sc storage cell index
+        \returns current file pointer, or nullptr if invalid
     */
-    virtual FILE *getFilePointer() { return dets[0]->getFilePointer(); };
+    virtual FILE *getFilePointer(int sc = 0) {
+        // Ensure storage cell index is valid
+        if (sc_to_threads.find(sc) == sc_to_threads.end() || sc_to_threads[sc].empty()) {
+            std::cerr << "Error: Invalid storage cell index " << sc << std::endl;
+            return nullptr;
+        }
+
+        // Return file pointer of the first thread in this storage cell
+        return dets[sc_to_threads[sc][0]] ? dets[sc_to_threads[sc][0]]->getFilePointer() : nullptr;
+    };
 
   protected:
     bool stop;
     const int nThreads;
     threadedAnalogDetector *dets[MAXTHREADS];
     analogDetector<uint16_t> *dd[MAXTHREADS];
-    int ithread;
-    int *image;
-    int *ff;
-    double *ped;
-    pthread_mutex_t fmutex;
+    //int ithread{0}; // Thread index
+    std::vector<int> thread_counters_by_sc{}; // Counters for threads for each storage cell 
+    int* image;
+    int* ff;
+    double* ped;
+    //pthread_mutex_t fmutex; //unused
+    std::unordered_map<int,std::vector<int>> sc_to_threads; // Maps storage cell -> vector of assigned thread ids
+    std::mutex map_mutex; // Ensure thread-safe access to the map
+    int nSC{1}; // Number of storage cells
+    
+    std::unordered_map<int,int*> sc_images; // Store images per storage cell
+    std::unordered_map<int,double*> sc_pedestals; // Store pedestal arrays per storage cell
+    std::unordered_map<int,double*> sc_pedestals_rms; // Store pedestal RMS arrays per storage cell
+    // at the moment, these maps could be avoided, but this implementation is more robust in allowing future changes
 };
 
 #endif

slsDetectorCalibration/multiThreadedCountingDetector.h

@@ -19,8 +19,8 @@ using namespace std;
 
 class multiThreadedCountingDetector : public multiThreadedAnalogDetector {
   public:
- multiThreadedCountingDetector(singlePhotonDetector *d, int n, int fs = 1000)
-        : multiThreadedAnalogDetector(d, n, fs){};
+ multiThreadedCountingDetector(singlePhotonDetector *d, int num_threads, int fs = 1000, int num_sc = 1)
+        : multiThreadedAnalogDetector(d, num_threads, fs, num_sc){};
     // virtual
     // ~multiThreadedCountingDetector{multiThreadedAnalogDetector::~multiThreadedAnalogDetector();};
     virtual double setNSigma(double n) {

slsDetectorCalibration/multiThreadedInterpolatingDetector.h

@@ -70,7 +70,6 @@ class multiThreadedInterpolatingDetector
         if (getInterpolation() == NULL)
             return multiThreadedAnalogDetector::getImage(nnx, nny, nsx, nsy);
         // if one interpolates, the whole image is stored in detector 0;
-        int *img;
         // int nnx, nny, ns;
         // int nnx, nny, ns;
         int nn = dets[0]->getImageSize(nnx, nny, nsx, nsy);
@@ -79,9 +78,9 @@ class multiThreadedInterpolatingDetector
             image = NULL;
         }
         image = new int[nn];
-        img = dets[0]->getImage();
+        int* tmp_img = dets[0]->getImage();
         for (int i = 0; i < nn; i++) {
-            image[i] = img[i];
+            image[i] = tmp_img[i];
         }
         return image;
     };

slsDetectorCalibration/singlePhotonDetector.h

@@ -51,7 +51,6 @@ class singlePhotonDetector : public analogDetector<uint16_t> {
 
 
     */
-
     singlePhotonDetector(slsDetectorData<uint16_t> *d, int csize = 3,
                          double nsigma = 5, int sign = 1,
                          commonModeSubtraction *cm = NULL, int nped = 1000,
@@ -60,10 +59,11 @@ class singlePhotonDetector : public analogDetector<uint16_t> {
         : analogDetector<uint16_t>(d, sign, cm, nped, nnx, nny, gm, gs),
           nDark(nd), eventMask(NULL), nSigma(nsigma), eMin(-1), eMax(-1),
           clusterSize(csize), clusterSizeY(csize), c2(1), c3(1), clusters(NULL),
-          quad(UNDEFINED_QUADRANT), tot(0), quadTot(0) {
+          quad(UNDEFINED_QUADRANT), tot(0), quadTot(0), ownsMutex(true) {    // The original object owns the mutex {
 
-        fm = new pthread_mutex_t;
-	pthread_mutex_init(fm, NULL);
+        //fm = new pthread_mutex_t;
+	    //pthread_mutex_init(fm, NULL);
+        fm = new std::mutex();
 
         eventMask = new eventType *[ny];
         //  val=new double*[ny];
@@ -86,24 +86,31 @@ class singlePhotonDetector : public analogDetector<uint16_t> {
         nphFrame = 0;
     };
     /**
-       destructor. Deletes the cluster structure, the pdestalSubtraction and the
-       image array
+       Destructor. Deletes the cluster structure, event mask, and destroys the mutex.
     */
     virtual ~singlePhotonDetector() {
-        delete[] clusters;
-        for (int i = 0; i < ny; i++)
-            delete[] eventMask[i];
-        delete[] eventMask;
+        // std::cout << "#### Debug: Destructing singlePhotonDetector! ####" << std::endl;
+        if (clusters) { delete[] clusters; clusters = nullptr; }
+        for (int i = 0; i < ny; i++) {
+            if (eventMask[i]) { delete[] eventMask[i]; eventMask[i] = nullptr; }
+        }
+        if (eventMask) { delete[] eventMask; eventMask = nullptr; }
+        if (ownsMutex) {
+            if (fm) {
+                //pthread_mutex_destroy(fm); // Destroy the mutex (not necessary with std::mutex)
+                delete fm; // Free the memory allocated for the mutex
+                fm = nullptr;  // Set the pointer to nullptr to avoid dangling pointer
+            }
+        }
+        
     };
 
     /**
-       copy constructor
+       pointer-based copy constructor (cloner)
        \param orig detector to be copied
-
     */
-
     singlePhotonDetector(singlePhotonDetector *orig)
-        : analogDetector<uint16_t>(orig) {
+        : analogDetector<uint16_t>(orig), fm(orig->fm), ownsMutex(false) {
 
         nDark = orig->nDark;
         myFile = orig->myFile;
@@ -126,11 +133,12 @@ class singlePhotonDetector : public analogDetector<uint16_t> {
         c3 = sqrt(clusterSizeY * clusterSize);
 
         clusters = new single_photon_hit[nx * ny];
+        //std::copy(orig->clusters, orig->clusters + (nx * ny), clusters); //possibly superfluous
 
         // cluster=clusters;
 
         setClusterSize(clusterSize);
-        fm = orig->fm;
+        //fm = orig->fm;
 
         quad = UNDEFINED_QUADRANT;
         tot = 0;
@@ -143,13 +151,67 @@ class singlePhotonDetector : public analogDetector<uint16_t> {
     }
 
     /**
-       duplicates the detector structure
-       \returns new single photon detector with same parameters
-
+     * copy constructor (deep copy), creates a new mutex
+     * stricly, TODO: Implement Rule of Five!
+     * (copy op=, move ctor, and move op= would need to be defined)
     */
-    virtual singlePhotonDetector *Clone() {
+    singlePhotonDetector(singlePhotonDetector const& other)
+        : analogDetector<uint16_t>(other), ownsMutex(true) {
+
+        //fm = new pthread_mutex_t; // create a new mutex
+	    //pthread_mutex_init(fm, NULL);
+        fm = new std::mutex();  // New unique mutex per copy
+
+        nDark = other.nDark;
+        myFile = other.myFile;
+
+        eventMask = new eventType *[ny];
+        for (int i = 0; i < ny; i++) {
+            eventMask[i] = new eventType[nx];
+            //std::copy(other.eventMask[i], other.eventMask[i] + nx, eventMask[i]);
+        }
+
+        eMin = other.eMin;
+        eMax = other.eMax;
+
+        nSigma = other.nSigma;
+        clusterSize = other.clusterSize;
+        clusterSizeY = other.clusterSizeY;
+
+        c2 = sqrt((clusterSizeY + 1) / 2 * (clusterSize + 1) / 2);
+        c3 = sqrt(clusterSizeY * clusterSize);
+
+        clusters = new single_photon_hit[nx * ny];
+        //std::copy(other.clusters, other.clusters + (nx * ny), clusters);
+
+        setClusterSize(clusterSize);
+
+        quad = other.quad;
+        tot = other.tot;
+        quadTot = other.quadTot;
+        gmap = other.gmap;
+        nphTot = other.nphTot;
+        nphFrame = other.nphFrame;
+   }
+
+    /**
+       Clones the detector structure
+       \returns new single photon detector with same parameters
+       that shares the mutex of the original
+    */
+    virtual singlePhotonDetector* Clone() {
         return new singlePhotonDetector(this);
     }
+    /**
+       Copies the detector structure
+       \returns new single photon detector with same parameters
+       that owns a new mutex
+    */
+    virtual singlePhotonDetector* Copy() {
+        return new singlePhotonDetector(*this);  // Calls the copy constructor
+    }
+    
+
     /** sets/gets number of rms threshold to detect photons
         \param n number of sigma to be set (0 or negative gets)
         \returns actual number of sigma parameter
@@ -381,7 +443,7 @@ class singlePhotonDetector : public analogDetector<uint16_t> {
 
         // int ir, ic;
         eventType ee;
-        double max = 0, tl = 0, tr = 0, bl = 0, br = 0, *v;
+        double max = 0, tl = 0, tr = 0, bl = 0, br = 0, v = 0;//, *v;
         int cm = 0;
         int good = 1;
         int ir, ic;
@@ -403,7 +465,8 @@ class singlePhotonDetector : public analogDetector<uint16_t> {
             cm = 1;
         }
 
-        double *val = new double[ny * nx];
+        //double *val = new double[ny * nx];
+        std::vector<double> val( ny * nx );
 
         for (int iy = ymin; iy < ymax; ++iy) {
             for (int ix = xmin; ix < xmax; ++ix) {
@@ -435,26 +498,34 @@ class singlePhotonDetector : public analogDetector<uint16_t> {
                             (ix + ic) >= 0 && (ix + ic) < nx) {
 
 
-			  if ((iy + ir) > iy && (ix + ic) > ix ) {
+			                if ((iy + ir) > iy && (ix + ic) > ix ) {
                             
-			    val[(iy + ir) * nx + ix + ic] =
-			      subtractPedestal(data, ix + ic, iy + ir, cm);
+			                    val[(iy + ir) * nx + ix + ic] =
+			                        subtractPedestal(data, ix + ic, iy + ir, cm);
 			    
 			    
-			  }
-			  v = &(val[(iy + ir) * nx + ix + ic]);
-			  tot += *v;
-			  if (ir <= 0 && ic <= 0)
-			    bl += *v;
-			  if (ir <= 0 && ic >= 0)
-			    br += *v;
-			  if (ir >= 0 && ic <= 0)
-			    tl += *v;
-			  if (ir >= 0 && ic >= 0)
-			  tr += *v;
-			  if (*v > max) //{
-			    max = *v;
-			  //}
+			                }
+			                //v = &(val[(iy + ir) * nx + ix + ic]);
+                            v = val[(iy + ir) * nx + ix + ic];
+			                //tot += *v;
+                            tot += v;
+			                if (ir <= 0 && ic <= 0)
+			                    bl += v;
+                                //bl += *v;
+			                if (ir <= 0 && ic >= 0)
+			                    br += v;
+                                //br += *v;
+			                if (ir >= 0 && ic <= 0)
+			                    tl += v;
+                                //tl += *v;
+			                if (ir >= 0 && ic >= 0)
+			                    tr += v;
+                                //tr += *v;
+			                //if (*v > max) //{
+                                //max = *v;
+                            if (v > max)
+			                    max = v;
+			                //}
                         }
                     }
                 }
@@ -524,19 +595,19 @@ class singlePhotonDetector : public analogDetector<uint16_t> {
                              ic < (clusterSize / 2) + 1; ic++) {
                             if ((iy + ir) >= 0 && (iy + ir) < ny &&
                                 (ix + ic) >= 0 && (ix + ic) < nx) {
-			      (clusters + nph)
+			                    (clusters + nph)
                                     ->set_data(val[(iy + ir) * nx + ix + ic],
                                                ic, ir);
-			      if (val[(iy + ir) * nx + ix + ic]>max) 
-				good=0;
-			    }
+			                    if (val[(iy + ir) * nx + ix + ic]>max) 
+				                    good=0;
+			                }
                         }
                     }
-		    if (good==0) {
-		      (clusters + nph)->print();
-		      cout << max << " " <<  val[iy * nx + ix] << endl;
-		    }
-		    //else (clusters + nph)->print();
+		            if (good==0) {
+		                (clusters + nph)->print();
+		                cout << max << " " <<  val[iy * nx + ix] << endl;
+		            }
+		            //else (clusters + nph)->print();
                     if (eMin > 0 && tot < eMin)
                         good = 0;
                     if (eMax > 0 && tot > eMax)
@@ -561,7 +632,7 @@ class singlePhotonDetector : public analogDetector<uint16_t> {
         // cout <<id << " **********************************"<< iframe << " " <<
         // det->getFrameNumber(data) << " " << nphFrame << endl;
         writeClusters(det->getFrameNumber(data));
-        delete[] val;
+        //delete[] val;
         return image;
     };
 
@@ -667,13 +738,14 @@ class singlePhotonDetector : public analogDetector<uint16_t> {
     void writeClusters(int fn) {
         if (myFile) {
             // cout << "++" << endl;
-            pthread_mutex_lock(fm);
+            //pthread_mutex_lock(fm); // This is dangerous! What if writeClusters() throws? Then the mutex is never unlocked!
+            std::lock_guard<std::mutex> lock(*fm); // safer, RAII-based locking
             //   cout <<"**********************************"<< fn << " " <<
             //   nphFrame << endl;
             writeClusters(myFile, clusters, nphFrame, fn);
             // for (int i=0; i<nphFrame; i++)
             //  (clusters+i)->write(myFile);
-            pthread_mutex_unlock(fm);
+            //pthread_mutex_unlock(fm); // unsafe
             // cout << "--" << endl;
         }
     };
@@ -711,7 +783,14 @@ class singlePhotonDetector : public analogDetector<uint16_t> {
         ema = eMax;
     };
 
-    void setMutex(pthread_mutex_t *m) { fm = m; };
+    //void setMutex(pthread_mutex_t *m) { fm = m; };
+    void setMutex(std::mutex* m) { 
+        if (ownsMutex && fm) {
+            delete fm;  // Cleanup old mutex
+        }
+        fm = m;
+        ownsMutex = false;
+    };
 
   protected:
     int nDark; /**< number of frames to be used at the beginning of the dataset
@@ -733,7 +812,9 @@ class singlePhotonDetector : public analogDetector<uint16_t> {
     int nphFrame;
 
     //    double **val;
-    pthread_mutex_t *fm;
+    //pthread_mutex_t* fm;      // Pointer to the shared mutex
+    std::mutex* fm;           // Pointer to the shared (or unique) mutex (safer version)
+    bool ownsMutex;           // Flag to indicate ownership
 };
 
 #endif