first draft

slsDetectorCalibration/dataStructures/HDF5File.cpp

@@ -0,0 +1,291 @@
+#include "HDF5File.h"
+
+#include "ansi.h"
+
+#include <string.h>
+
+HDF5File::HDF5File () {
+  //InitializeParameters();
+}
+
+HDF5File::~HDF5File () {
+	if(frame_index_list)
+		delete [] frame_index_list;
+	if(current_image)
+		delete [] current_image;
+}
+
+void HDF5File::InitializeParameters () {
+  /*
+    memset(file_name, 0, MAX_STR_LENGTH); //awkward, initializes all file_name characters to 0
+    file = -1;
+    dataspace = -1;
+    //memspace = -1;
+    dataset = -1;
+    number_of_frames = 0;
+    frame_index_list = NULL;
+    current_image = NULL;
+  */
+  for (int i = 0; i < RANK; ++i) {
+    file_dims[i] = 0; //also awkward
+    chunk_dims[i] = 0;
+    frame_offset[i] = 0;
+  }
+}
+
+
+int HDF5File::OpenResources (const char* const fname, bool validate) {
+	// Initialize
+	//InitializeParameters();
+
+
+	// Open File
+	file = H5Fopen (fname, H5F_ACC_RDONLY, H5P_DEFAULT);
+	if (file < 0) {
+	  cprintf(RED,"could not open hdf5 file\n");
+	  return 0;
+	}
+	cprintf(BLUE, "Opened File: %s\n", fname);
+
+
+	// Open Dataset
+	dataset = H5Dopen2 (file, DATA_DATASETNAME, H5P_DEFAULT);
+	if (dataset < 0){
+	  cprintf(RED,"could not open dataset\n");
+	  CloseResources ();
+	  return 0;
+	}
+	cprintf(BLUE, "Opened Dataset: %s\n", DATA_DATASETNAME);
+	
+
+	// Create Dataspace
+	dataspace = H5Dget_space (dataset);
+	if (dataspace < 0){
+	  cprintf(RED,"could not open dataspace\n");
+	  CloseResources ();
+	  return 0;
+	}
+
+
+	// Get Dimensions
+	int rank = H5Sget_simple_extent_dims (dataspace, file_dims, NULL);
+	cprintf (BLUE, "Number of Images: %llu\n", file_dims[Z]);
+
+	// validate file dimensions
+	if (validate) {
+
+		// validate rank
+		if(rank != RANK) {
+			cprintf(RED,"rank found %d. Expected %d\n", rank, RANK);
+			CloseResources ();
+			return 0;
+		}
+
+		// validate file dimensions of x and y
+		if (file_dims[X] != DEFAULT_X_DIMS) {
+			cprintf(RED,"file dimensions of x found %llu. Expected %d\n", file_dims[X], DEFAULT_X_DIMS);
+			CloseResources ();
+			return 0;
+		}
+		if (file_dims[Y] != DEFAULT_Y_DIMS) {
+			cprintf(RED,"file dimensions of y found %llu. Expected %d\n", file_dims[Y], DEFAULT_Y_DIMS);
+			CloseResources ();
+			return 0;
+		}
+		cprintf(GREEN, "File rank & dimensions validated. "
+				"Rank: %d, Dimensions: %llu x %llu x %llu\n",
+					rank, file_dims[Z], file_dims[Y], file_dims[X]);
+	}
+
+
+
+	// Get layout
+	hid_t cparms = H5Dget_create_plist(dataset);
+
+	// validate chunk layout
+	if (validate) {
+		if (H5D_CHUNKED != H5Pget_layout (cparms))  {
+			cprintf(RED,"not chunked data file\n");
+			H5Pclose(cparms);
+			CloseResources ();
+			return 0;
+		}
+		cprintf(GREEN, "Chunk layout validated\n");
+	}
+
+	// Get Chunk Dimensions
+	int rank_chunk = H5Pget_chunk (cparms, RANK, chunk_dims);
+
+	// validate dimensions
+	if (validate) {
+
+		// validate rank
+		if(rank_chunk != RANK) {
+			cprintf(RED,"chunk rank found %d. Expected %d\n", rank, RANK);
+			H5Pclose(cparms);
+			CloseResources ();
+			return 0;
+		}
+
+		// validate file dimensions of x, y and z
+		if (chunk_dims[X] != DEFAULT_CHUNK_X_DIMS) {
+			cprintf(RED,"chunk dimensions of x found %llu. Expected %d\n", chunk_dims[X], DEFAULT_CHUNK_X_DIMS);
+			H5Pclose(cparms);
+			CloseResources ();
+			return 0;
+		}
+		if (chunk_dims[Y] != DEFAULT_CHUNK_Y_DIMS) {
+			cprintf(RED,"chunk dimensions of y found %llu. Expected %d\n", chunk_dims[Y], DEFAULT_CHUNK_Y_DIMS);
+			H5Pclose(cparms);
+			CloseResources ();
+			return 0;
+		}
+		/*if (chunk_dims[Z] != DEFAULT_CHUNK_Z_DIMS) {
+			cprintf(RED,"chunk dimensions of z found %llu. Expected %d\n", chunk_dims[Z], DEFAULT_CHUNK_Z_DIMS);
+			H5Pclose(cparms);
+			CloseResources ();
+			return 0;
+		}*/
+
+		cprintf(GREEN, "Chunk rank & dimensions validated. "
+					"Rank: %d, Dimensions: %llu x %llu x %llu\n",
+					rank_chunk, chunk_dims[Z], chunk_dims[Y], chunk_dims[X]);
+
+	}
+
+	H5Pclose (cparms);
+
+	// allocate chunk memory
+	current_image = new uint16_t[chunk_dims[Z]*DEFAULT_CHUNK_Y_DIMS*DEFAULT_CHUNK_X_DIMS];
+	//current_image = new uint16_t[DEFAULT_X_DIMS*DEFAULT_Y_DIMS];
+
+	// Define memory space
+	//memspace = H5Screate_simple (RANK, chunk_dims, NULL);
+
+
+	// Get all the frame numbers
+	// Open frame index dataset
+	hid_t fi_dataset = H5Dopen2 (file, INDEX_DATASETNAME, H5P_DEFAULT);
+	if (fi_dataset < 0){
+		cprintf (RED,"could not open frame index dataset %s\n", INDEX_DATASETNAME);
+		CloseResources ();
+		return 0;
+	}
+
+	// validate size of frame index dataset
+	if (validate) {
+		hsize_t fi_dims[2];
+		hid_t fi_dataspace = H5Dget_space (fi_dataset);
+		int fi_rank = H5Sget_simple_extent_dims (fi_dataspace, fi_dims, NULL);
+
+		// validate rank
+		if(fi_rank != 2) {
+			cprintf(RED,"Frame index dataset rank found %d. Expected %d\n", fi_rank, 2);
+			H5Sclose (fi_dataspace);
+			H5Dclose (fi_dataset);
+			CloseResources ();
+			return 0;
+		}
+
+		// validate size
+		if (fi_dims[Z] != file_dims[Z]) {
+			cprintf (RED,"Frame index dimensions of z found %llu. Expected %llu\n", fi_dims[Z], file_dims[Z]);
+			H5Sclose (fi_dataspace);
+			H5Dclose (fi_dataset);
+			CloseResources ();
+			return 0;
+		}
+		H5Sclose (fi_dataspace);
+	}
+
+
+	// allocate frame index memory
+	frame_index_list = new unsigned int[file_dims[Z]];
+
+	//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_U32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, frame_index_list) < 0) {
+		cprintf (RED,"Could not read frame index dataset %s\n", INDEX_DATASETNAME);
+		H5Dclose (fi_dataset);
+		CloseResources ();
+	}
+	H5Dclose(fi_dataset);
+
+	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;
+	}
+	//if (memspace >= 0)      H5Sclose(memspace); // VH: I am getting memory leaks after recompilation
+}
+
+
+int HDF5File::ReadImage (uint16_t** image, int& iFrame) {
+
+	// no images in this frame
+	if (frame_index_list[frame_offset[Z]] == 0) {
+		cprintf (RED,"No images in this frame offset %llu\n", frame_offset[Z]);
+		CloseResources ();
+		return -99;
+	}
+
+	if (frame_offset[Z] == file_dims[Z]-1) {
+		printf("end of file\n");
+		return -1;
+	}
+
+	hsize_t frame_size[RANK] = {1, file_dims[X], file_dims[Y]};
+
+	// Define memory space
+	hid_t memspace = H5Screate_simple (RANK, frame_size, NULL);
+
+	// create hyperslab
+	// If I understand correctly, this puts dataspace such that we read the correct frame
+	if (H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, frame_offset, NULL, frame_size, NULL) < 0 ) {
+		cprintf (RED,"Could not create hyperslab for frame count %llu\n", frame_offset[Z]);
+		CloseResources ();
+		return -99;
+	}
+
+	// read dataset into current_image
+	if (H5Dread(dataset, H5T_STD_U16LE, memspace, dataspace, H5P_DEFAULT, current_image) < 0 ) {
+		cprintf (RED,"Could not read dataset for frame count %llu\n", frame_offset[Z]);
+		CloseResources ();
+		return -99;
+	}
+	*image = current_image;
+
+	// return frame number and then increment frame count number
+	unsigned int retval = frame_index_list[frame_offset[Z]];
+	iFrame = (int)frame_offset[Z];
+	++frame_offset[Z];
+	return retval;
+}
+
+
+
+void HDF5File::PrintCurrentImage () {
+	printf("\n");
+	printf("Frame %llu, Image: %d\n", frame_offset[Z]-1, frame_index_list[frame_offset[Z]-1]);
+
+	unsigned long long int size = file_dims[Y] * file_dims[X];
+	for (unsigned int i = 0; i < size; ++i){
+		printf("%u ", current_image[i]);
+		if (!((i+1) % file_dims[X] ))
+			printf("\n\n");
+	}
+	printf("\n\n\n\n");
+}
+

slsDetectorCalibration/dataStructures/HDF5File.h

@@ -0,0 +1,119 @@
+#pragma once
+/************************************************
+ * @file HDF5Fle.h
+ * @short functions to open/close/read HDF5 File
+ ***********************************************/
+/**
+ *@short functions to open/close/read HDF5 File
+ */
+
+
+#include "hdf5.h"
+#include "hdf5_hl.h"
+
+
+#define MAX_STR_LENGTH	1000
+
+#define RANK	3                       // Dimension of the image
+#define DEFAULT_Z_DIMS	10000
+#define DEFAULT_Y_DIMS	1024
+#define DEFAULT_X_DIMS	512
+//#define DEFAULT_S_DIMS	1       // Storage cells
+
+#define DEFAULT_CHUNK_Z_DIMS	1
+#define DEFAULT_CHUNK_Y_DIMS	1024
+#define DEFAULT_CHUNK_X_DIMS	512
+//#define DEFAULT_CHUNK_S_DIMS	1
+
+/** Assuming each chunk is one image 1024 x 512*/
+
+#define DATA_DATASETNAME	"/data/JF18T01V01/data" //Furka JF
+#define INDEX_DATASETNAME	"/data/JF18T01V01/frame_index"
+
+enum{Z,X,Y}; //S is the storage cell
+
+class HDF5File {
+
+public:
+
+
+	/**
+	 * Constructor
+	 */
+	HDF5File ();
+
+	/**
+	 * Destructor
+	 */
+	~HDF5File ();
+
+	/**
+	 * Initialize Parameters for each new file
+	 */
+	void InitializeParameters ();
+
+	/**
+	 * 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
+	 * @returns frame number read,
+	 */
+	int ReadImage (uint16_t** image, int& iFrame);
+
+	/**
+	 * Print current image in memory
+	 */
+	void PrintCurrentImage ();
+
+private:
+
+	/** file name */
+	std::string file_name{};
+
+	/** file handle */
+	hid_t file{};
+
+	/** dataspace handle */
+	hid_t dataspace{};
+
+	/** memory space handle */
+	//hid_t memspace;
+
+	/** dataset handle */
+	hid_t dataset{};
+
+
+	/** file dimensions */
+	hsize_t file_dims[RANK]{}; //static array (dimensions are known) //I think the {} initialization should work...
+
+	/** chunk dimensions */
+	hsize_t chunk_dims[RANK]{};
+
+	/** number of frames */
+	unsigned int number_of_frames{};
+
+	/** frame index list */
+	unsigned int* frame_index_list{NULL}; //dynamic array
+
+
+	/** current image */
+	uint16_t* current_image{NULL}; //dynamic array
+	//uint16_t current_chunk[DEFAULT_CHUNK_Z_DIMS][DEFAULT_CHUNK_Y_DIMS][DEFAULT_CHUNK_X_DIMS];
+
+	/** current frame offset */
+	hsize_t frame_offset[RANK]{}; //array (frame_offset[Z], 0, 0) I believe
+
+};

slsDetectorCalibration/dataStructures/jungfrauLGADStrixelsDataQuadH5.h

@@ -0,0 +1,426 @@
+// 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.h" //this includes hdf5.h and hdf5_hl.h
+#include "HDF5File.cpp"
+
+// #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 ( ipx>=xmin && ipx<=xmax && ipy>=ymin && ipy <=ymax )
+		    dataMap[iy][ix] =
+		      (globalROI.nc * (ipy - globalROI.ymin) + (ipx - globalROI.xmin)) * 2;
+		} else {
+		  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);
+    }
+
+    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 ) {
+        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); //mayb another value is safer
+            }
+        }
+
+	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;
+    };
+
+    /**
+
+    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(HDF5File& hfile, 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/jungfrauExecutables/CMakeLists.txt

@@ -9,6 +9,16 @@ set(JUNGFRAU_EXECUTABLES)
 find_package(fmt REQUIRED)
 find_package(nlohmann_json 3.11.2 REQUIRED)
 
+# 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)
@@ -29,6 +39,17 @@ 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)
+       target_link_libraries(jungfrauRawDataProcessStrxQuadH5 PUBLIC ${HDF5_LIBRARIES})
+       target_include_directories(jungfrauRawDataProcessStrxQuadH5 PUBLIC ${HDF5_INCLUDE_DIRS}) #${CMAKE_INSTALL_PREFIX}/include)
+       list(APPEND JUNGFRAU_EXECUTABLES jungfrauRawDataProcessStrxQuadH5)
+    endif ()
+endif (SLS_USE_HDF5)
 
 # jungfrauRawDataProcessStrxChip1
 add_executable(jungfrauRawDataProcessStrxChip1 jungfrauRawDataProcess.cpp)

slsDetectorCalibration/jungfrauExecutables/jungfrauRawDataProcess_filetxtH5.cpp

@@ -0,0 +1,558 @@
+// 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
+
+#if !defined JFSTRX && !defined JFSTRXQ && !defined JFSTRXQH5 && !defined JFSTRXOLD && !defined JFSTRXCHIP1 &&           \
+    !defined JFSTRXCHIP6 && !defined CHIP
+#ifndef MODULE
+#include "jungfrauHighZSingleChipData.h"
+#endif
+#ifdef MODULE
+#include "jungfrauModuleData.h"
+#endif
+#endif
+
+#ifdef CHIP
+#include "jungfrauSingleChipData.h"
+#endif
+
+#ifdef JFSTRX
+#include "jungfrauLGADStrixelsData_new.h"
+#endif
+#ifdef JFSTRXQ
+#include "jungfrauLGADStrixelsDataQuad.h"
+#endif
+#ifdef JFSTRXQH5
+#include "jungfrauLGADStrixelsDataQuadH5.h"
+#endif
+#if defined JFSTRXCHIP1 || defined JFSTRXCHIP6
+#include "jungfrauLGADStrixelsDataSingleChip.h"
+#endif
+#ifdef JFSTRXOLD
+#include "jungfrauStrixelsHalfModuleOldDesign.h"
+#endif
+
+#include "multiThreadedCountingDetector.h"
+#include "singlePhotonDetector.h"
+
+#include <fstream>
+#include <map>
+#include <stdio.h>
+#include <sys/stat.h>
+
+#include <ctime>
+#include <fmt/core.h>
+
+#include <nlohmann/json.hpp>
+using json = nlohmann::json;
+
+
+std::string getRootString( const std::string& 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( const std::string& dir, const std::string& fprefix="run", const std::string& fsuffix="", const std::string& fext="raw", int outfilecounter=-1 ) {
+  if (outfilecounter >= 0)
+    return fmt::format("{:s}/{:s}_{:s}_f{:05d}.{:s}", dir, fprefix, fsuffix, outfilecounter, fext);
+  else if (fsuffix.length()!=0)
+    return fmt::format("{:s}/{:s}_{:s}.{:s}", dir, fprefix, fsuffix, fext);
+  else
+    return fmt::format("{:s}/{:s}.{:s}", dir, fprefix, fext);
+}
+
+
+//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 < 11) {
+        std::cout
+            << "Usage is " << argv[0]
+            << " 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
+            << "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 fifosize = 1000;
+    int nthreads = 10;
+    int csize = 3; // 3
+    int nsigma = 5;
+    int nped = 10000;
+
+    int cf = 0;
+
+    double *gainmap = NULL;
+    // float *gm;
+
+    int ff, np;
+    // cout << " data size is " << dsize;
+
+    const std::string txtfilename(argv[1]);
+    const std::string outdir(argv[2]);
+    const std::string jsonmastername(argv[3]);
+    const std::string pedfilename(argv[4]);
+
+    double thr = 0;
+    double thr1 = 1;
+    thr = atof(argv[9]);
+
+    int nframes = 0;
+    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{};
+    //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);
+	  filenames.emplace_back(line);
+	  std::cout << line << " 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;
+
+    // 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 JFSTRXQ && !defined JFSTRXQH5 && !defined JFSTRXOLD && !defined JFSTRXCHIP1 &&           \
+    !defined JFSTRXCHIP6 && !defined CHIP
+#ifndef MODULE
+    jungfrauHighZSingleChipData *decoder = new jungfrauHighZSingleChipData();
+    int nx = 256, ny = 256;
+#endif
+#ifdef MODULE
+    jungfrauModuleData *decoder = new jungfrauModuleData();
+    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
+    std::cout << "Jungfrau strixel full module readout" << std::endl;
+
+#ifndef ALDO
+    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(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 JFSTRXQH5
+    std::cout << "Jungfrau strixel quad h5" << std::endl;
+    jungfrauLGADStrixelsDataQuadH5 *decoder =
+        new jungfrauLGADStrixelsDataQuadH5(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 =
+        new jungfrauLGADStrixelsDataSingleChip(1);
+    int nx = 256 / 3, ny = 256 * 5;
+#endif
+#ifdef JFSTRXCHIP6
+    std::cout << "Jungfrau strixel LGAD single chip 6" << std::endl;
+    jungfrauLGADStrixelsDataSingleChip *decoder =
+        new jungfrauLGADStrixelsDataSingleChip(6);
+    int nx = 256 / 3, ny = 256 * 5;
+#endif
+#ifdef JFSTRXOLD
+    std::cout << "Jungfrau strixels old design" << std::endl;
+    jungfrauStrixelsHalfModuleOldDesign *decoder =
+        new jungfrauStrixelsHalfModuleOldDesign();
+    int nx = 1024 * 3, ny = 512 / 3;
+#endif
+
+    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;
+    }
+    if ( ymin == ymax ) {
+      ymin = 0;
+      ymax = ny;
+    }
+    std::cout << xmin << " " << xmax << " " << ymin << " " << ymax << " "
+              << std::endl;
+    */
+
+    /*
+    char *gainfname = NULL;
+    if (argc > 14) {
+        gainfname = argv[14];
+        std::cout << "Gain map file name is: " << gainfname << 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, gainmap, NULL);
+
+    /*
+    if (gainfname) {
+
+        if (filter->readGainMap(gainfname))
+            std::cout << "using gain map " << gainfname << std::endl;
+        else
+            std::cout << "Could not open gain map " << gainfname << std::endl;
+    } else
+    */
+    thr = 0.15 * thr;
+    filter->newDataSet();
+    // 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);
+    mt->setClusterSize(csize, csize);
+
+#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);
+
+    //  cout << "mt " << endl;
+
+    int ifr = 0;
+
+    if (pedfilename.length()>1) {
+
+      std::string froot = getRootString(pedfilename);
+
+      std::cout << "PEDESTAL " << std::endl;
+
+        if (pedfilename.find(".tif") == std::string::npos) {
+	  const std::string fname(pedfilename);
+	  std::cout << fname << std::endl;
+	  std::time(&end_time);
+	  std::cout << "aaa " << std::ctime(&end_time) << std::endl;
+
+	  mt->setFrameMode(ePedestal);
+
+	  std::ifstream pedefile(fname, ios::in | ios::binary);
+	  //      //open file
+	  if (pedefile.is_open()) {
+	    std::cout << "bbbb " << std::ctime(&end_time) << std::endl;
+	    
+	    ff = -1;
+	    while (decoder->readNextFrame(pedefile, ff, np, buff)) {
+	      // if (np == 40) {
+	      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 << " " << ff << " " << np
+			  << std::endl;
+	      } // else
+	      //std::cout << ifr << " " << ff << " " << np << std::endl;
+	      if (ifr >= 1000)
+		break;
+	      ff = -1;
+	    }
+	    pedefile.close();
+	    while (mt->isBusy()) {
+	      ;
+	    }
+
+	    std::cout << "froot " << froot << std::endl;
+	    auto imgfname = createFileName( outdir, froot, "ped", "tiff");
+	    mt->writePedestal(imgfname.c_str());
+	    imgfname = createFileName( outdir, froot, "rms", "tiff");
+	    mt->writePedestalRMS(imgfname.c_str());
+
+	  } else
+	    std::cout << "Could not open pedestal file " << fname
+		      << " for reading " << std::endl;
+        } else {
+	  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 = NULL;
+    std::ifstream filebin{};
+    
+    //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 ";
+      std::string fsuffix{};
+      const std::string fprefix( getRootString(filenames[ifile]) );
+      std::string imgfname( createFileName( outdir, fprefix, fsuffix, "tiff" ) );
+      const std::string cfname( createFileName( outdir, fprefix, fsuffix, "clust" ) );
+      std::cout << filenames[ifile] << " ";
+      std::cout << imgfname << std::endl;
+      std::time(&end_time);
+      std::cout << std::ctime(&end_time) << std::endl;
+
+      std::ifstream filebin(filenames[ifile], ios::in | ios::binary);
+      //      //open file
+      ioutfile = 0;
+      if (filebin.is_open()) {
+	if (thr <= 0 && cf != 0) { // cluster finder
+	  if (of == NULL) {
+	    of = fopen(cfname.c_str(), "w");
+	    if (of) {
+	      mt->setFilePointer(of);
+	      std::cout << "file pointer set " << std::endl;
+	    } else {
+	      std::cout << "Could not open " << cfname
+			<< " for writing " << std::endl;
+	      mt->setFilePointer(NULL);
+	      return 1;
+	    }
+	  }
+	}
+	//     //while read frame
+	ff = -1;
+	ifr = 0;
+	while (decoder->readNextFrame(filebin, ff, np, buff)) {
+	  //  if (np == 40) {
+	  //         //push
+	  
+	  if ((ifr + 1) % 100 == 0) {
+	    std::cout << " ****"
+		      << decoder->getValue(buff, 20, 20); // << std::endl;
+	  }
+	  mt->pushData(buff);
+	  // 	//         //pop
+	  mt->nextThread();
+	  mt->popFree(buff);
+	  
+	  ifr++;
+	  if (ifr % 100 == 0)
+	    std::cout << " " << ifr << " " << ff << std::endl;
+	  if (nframes > 0) {
+	    if (ifr % nframes == 0) {
+	      imgfname = createFileName( outdir, fprefix, fsuffix, "tiff", ioutfile );
+	      mt->writeImage(imgfname.c_str(), thr1);
+	      mt->clearImage();
+	      ++ioutfile;
+	    }
+	  }
+	  // } else
+	  //std::cout << ifr << " " << ff << " " << np << std::endl;
+	  ff = -1;
+	}
+	std::cout << "aa --" << std::endl;
+	filebin.close();
+	std::cout << "bb --" << std::endl;
+	while (mt->isBusy()) {
+	  ;
+	}
+	std::cout << "cc --" << std::endl;
+	if (nframes >= 0) {
+	  if (nframes > 0)
+	    imgfname = createFileName( outdir, fprefix, fsuffix, "tiff", ioutfile );
+	  std::cout << "Writing tiff to " << imgfname << " " << thr1
+		    << std::endl;
+	  mt->writeImage(imgfname.c_str(), thr1);
+	  mt->clearImage();
+	  if (of) {
+	    fclose(of);
+	    of = NULL;
+	    mt->setFilePointer(NULL);
+	  }
+	}
+	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 fname(argv[10]);
+      std::string fprefix( getRootString(filenames[0]) ); //This might by a non-ideal name choice for that file
+      std::string imgfname( createFileName( outdir, fprefix, "sum", "tiff" ) );
+      std::cout << "Writing tiff to " << imgfname << " " << thr1 << std::endl;
+      mt->writeImage(imgfname.c_str(), thr1);
+    }
+
+    return 0;
+}