From a6a4ea7f218dbc1a65ec4bd01423db4f9ba412b6 Mon Sep 17 00:00:00 2001 From: vhinger182 Date: Wed, 7 Aug 2024 11:05:23 +0200 Subject: [PATCH] first draft --- .../dataStructures/HDF5File.cpp | 291 +++++++++ .../dataStructures/HDF5File.h | 119 ++++ .../jungfrauLGADStrixelsDataQuadH5.h | 426 +++++++++++++ .../jungfrauExecutables/CMakeLists.txt | 21 + .../jungfrauRawDataProcess_filetxtH5.cpp | 558 ++++++++++++++++++ 5 files changed, 1415 insertions(+) create mode 100644 slsDetectorCalibration/dataStructures/HDF5File.cpp create mode 100644 slsDetectorCalibration/dataStructures/HDF5File.h create mode 100644 slsDetectorCalibration/dataStructures/jungfrauLGADStrixelsDataQuadH5.h create mode 100644 slsDetectorCalibration/jungfrauExecutables/jungfrauRawDataProcess_filetxtH5.cpp diff --git a/slsDetectorCalibration/dataStructures/HDF5File.cpp b/slsDetectorCalibration/dataStructures/HDF5File.cpp new file mode 100644 index 000000000..62e4c577b --- /dev/null +++ b/slsDetectorCalibration/dataStructures/HDF5File.cpp @@ -0,0 +1,291 @@ +#include "HDF5File.h" + +#include "ansi.h" + +#include + +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"); +} + diff --git a/slsDetectorCalibration/dataStructures/HDF5File.h b/slsDetectorCalibration/dataStructures/HDF5File.h new file mode 100644 index 000000000..803dc0293 --- /dev/null +++ b/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 + +}; diff --git a/slsDetectorCalibration/dataStructures/jungfrauLGADStrixelsDataQuadH5.h b/slsDetectorCalibration/dataStructures/jungfrauLGADStrixelsDataQuadH5.h new file mode 100644 index 000000000..9e1117695 --- /dev/null +++ b/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 +#include +#include + +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 { + + 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 mods{ 0, 1, 2 }; + + void reverseVector( std::vector& 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 <=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( + 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 diff --git a/slsDetectorCalibration/jungfrauExecutables/CMakeLists.txt b/slsDetectorCalibration/jungfrauExecutables/CMakeLists.txt index 931ecfa2b..a97df3242 100644 --- a/slsDetectorCalibration/jungfrauExecutables/CMakeLists.txt +++ b/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) diff --git a/slsDetectorCalibration/jungfrauExecutables/jungfrauRawDataProcess_filetxtH5.cpp b/slsDetectorCalibration/jungfrauExecutables/jungfrauRawDataProcess_filetxtH5.cpp new file mode 100644 index 000000000..c7b6fd111 --- /dev/null +++ b/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 +#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 +#include +#include +#include + +#include +#include + +#include +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 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 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; +}