diff --git a/slsDetectorCalibration/dataStructures/jungfrauLGADStrixelsData.h b/slsDetectorCalibration/dataStructures/jungfrauLGADStrixelsData.h index e5700f985..863ea3f71 100644 --- a/slsDetectorCalibration/dataStructures/jungfrauLGADStrixelsData.h +++ b/slsDetectorCalibration/dataStructures/jungfrauLGADStrixelsData.h @@ -2,7 +2,11 @@ // Copyright (C) 2021 Contributors to the SLS Detector Package #ifndef JUNGFRAULGADSTRIXELDATA_H #define JUNGFRAULGADSTRIXELDATA_H +#ifdef CINT +#include "sls/sls_detector_defs_CINT.h" +#else #include "sls/sls_detector_defs.h" +#endif #include "slsDetectorData.h" /* @@ -32,7 +36,7 @@ typedef struct { uint64_t bunchNumber; /**< is the frame number */ uint64_t pre; /**< something */ -} jf_header; +} jf_header; //Aldo's header using namespace std; @@ -41,7 +45,11 @@ class jungfrauLGADStrixelsData : public slsDetectorData { private: int iframe; +#ifdef ALDO //VH + using header = jf_header; //VH +#else //VH using header = sls::defs::sls_receiver_header; +#endif //VH public: /** Implements the slsReceiverData structure for the moench02 prototype read @@ -53,6 +61,9 @@ class jungfrauLGADStrixelsData : public slsDetectorData { : slsDetectorData(1024/5, 512*5, 512 * 1024 * 2 + sizeof(header)) { cout << "aaa" << endl; +#ifdef ALDO //VH + cout<< "using reduced jf_header" << endl; //VH +#endif //VH for (int ix = 0; ix < 1024/5; ix++) { for (int iy = 0; iy < 512*5; iy++) { dataMap[iy][ix] = sizeof(header);//+ ( 1024 * 5 + 300) * 2; //somewhere on the guardring of the LGAD @@ -61,44 +72,51 @@ class jungfrauLGADStrixelsData : public slsDetectorData { #endif } } - int x0=256+10, x1=256+246; - int y0=10, y1=256-10; + + //chip1 + /* + * TL;DR comments: y0 is too high by 1, group 1 and group 2 are by one row too short + * group34 by 2 rows, x is by one column too short + * NOTE: If x0, x1, y0 are changed, likely also ox (and oy and ooy) will be affected! + */ + + cout << "G0" << endl; //chip coordinates of chip1 group1: x=255+10 to x=255+246, y=10 to y=64 + //9 pixels guard ring, bonding shift by one pixel in y, one square pixel in x on the left + + int x0=256+10, x1=256+246; //excludes first column (in chip coordinates) + int y0=10, y1=256-10; //y1 does nothing int ix,iy; int ox=0, oy=0, ooy=0; ox=0; - cout << "G0" << endl; - - //chip1 - for (int ipx=x0; ipx { } //chip 6 + /* + * TL;DR comments: y0 is too high by 3, group34 and group 1 are by one row too short + * x is by two columns too short + * NOTE: If x0, x1, y0 are changed, likely also ox (and oy and ooy) will be affected! + */ + cout << "G0" << endl; //chip coordinates of chip6 group34: x=255+256+9 to x=255+256+244, y=255+8 to y=255+126 + //9 pixels guard ring, bonding shift by one pixel in -y, 2 square pixels in x on the right x0=256*2+10; y0=256+10; x1=256*2+246; ooy=256*5; oy=0; ox=1; - cout << "G0" << endl; for (int ipx=x0; ipx { int getFrameNumber(char *buff) { - return ((header *)buff)->detHeader.frameNumber; +#ifdef ALDO //VH + return ((header *)buff)->bunchNumber; //VH +#else //VH + return ((header *)buff)->detHeader.frameNumber; +#endif //VH }; /** @@ -219,7 +248,11 @@ class jungfrauLGADStrixelsData : public slsDetectorData { */ int getPacketNumber(char *buff) { - return ((header *)buff)->detHeader.packetNumber; +#ifdef ALDO //VH + return -1; //VH //TODO: Keep in mind in case of bugs! +#else //VH + return ((header *)buff)->detHeader.packetNumber; +#endif //VH }; @@ -265,8 +298,6 @@ class jungfrauLGADStrixelsData : public slsDetectorData { 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 */ diff --git a/slsDetectorCalibration/dataStructures/jungfrauLGADStrixelsDataSingleChip.h b/slsDetectorCalibration/dataStructures/jungfrauLGADStrixelsDataSingleChip.h new file mode 100644 index 000000000..ad67ca51f --- /dev/null +++ b/slsDetectorCalibration/dataStructures/jungfrauLGADStrixelsDataSingleChip.h @@ -0,0 +1,305 @@ +// SPDX-License-Identifier: LGPL-3.0-or-other +// Copyright (C) 2021 Contributors to the SLS Detector Package +#ifndef JUNGFRAULGADSTRIXELSDATASINGLECHIP_H +#define JUNGFRAULGADSTRIXELSDATASINGLECHIP_H +#ifdef CINT +#include "sls/sls_detector_defs_CINT.h" +#else +#include "sls/sls_detector_defs.h" +#endif +#include "slsDetectorData.h" + +/* +/afs/psi.ch/project/mythen/Anna/slsDetectorPackageDeveloperMpc2011/slsDetectorCalibration/jungfrauExecutables +make -f Makefile.rawdataprocess jungfrauRawDataProcessStrx + ../dataStructures/jungfrauLGADStrixelsData.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 +*/ + +namespace strixelSingleChip { + constexpr int nc_chip = 256; + constexpr int nr_chip = 256; + constexpr int gr = 9; + + //Group 1: 25um pitch, groups of 3, 1 column of square pixels + constexpr int g1_ncols{ (nc_chip-(2*gr)-1)/3 }; //79 + constexpr int g1_nrows{ ( (nr_chip/4)-gr )*3 }; //165 + + //Group 2: 15um pitch, groups of 5, 3 columns of square pixels + constexpr int g2_ncols{ (nc_chip-(2*gr)-3)/5 }; //47 + constexpr int g2_nrows{ (nr_chip/4)*5 }; //320 + + //Group 3: 18.75um pitch, groups of 4, 2 columns of square pixels (double the size of the other groups) + constexpr int g3_ncols{ (nc_chip-(2*gr)-2)/4 }; //59 + constexpr int g3_nrows{ ( ((nr_chip/4)*2)-gr )*4 }; //476 + + constexpr int nc_strixel = 2*gr + 1 + g1_ncols; //group 1 is the "longest" group in x and has one extra square pixel + constexpr int nr_strixel = 2*gr + g1_nrows + g2_nrows + g3_nrows; +} + +typedef struct { + uint64_t bunchNumber; /**< is the frame number */ + uint64_t pre; /**< something */ + +} jf_header; //Aldo's header + + +using namespace strixelSingleChip; +class jungfrauLGADStrixelsDataSingleChip : public slsDetectorData { + + private: + int iframe; + int mchip; + + void remapGroup( const int group ) { + int ix, iy; + int x0, y0, x1, y1, shifty; + int multiplicator; + + switch (group) { + default: + case 1: + multiplicator = 3; + break; + case 2: + multiplicator = 5; + break; + case 3: + multiplicator = 4; + break; + } + + if ( mchip == 1 ) { + switch (group) { + default: + case 1: + x0 = 10; + x1 = 247; + y0 = 10; + y1 = 65; + shifty = 0; + break; + case 2: + x0 = 12; + x1 = 247; + y0 = 65; + y1 = 129; + shifty = ( (nr_chip/4)-gr )*3; + break; + case 3: + x0 = 11; + x1 = 247; + y0 = 129; + y1 = 248; + shifty = ( (nr_chip/4)-gr )*3 + (nr_chip/4)*5; + break; + } + } + + if ( mchip == 6 ) { + switch (group) { + default: + case 1: + x0 = 9; + x1 = 246; + y0 = 191; + y1 = 246; + shifty = ( (nr_chip/4)-gr+(nr_chip/4) )*4 + (nr_chip/4)*5; + break; + case 2: + x0 = 9; + x1 = 244; + y0 = 127; + y1 = 191; + shifty = ( (nr_chip/4)-gr+(nr_chip/4) )*4; + break; + case 3: + x0 = 9; + x1 = 245; + y0 = 8; + y1 = 127; + shifty = 0; + break; + } + } + + //remapping loop + for ( int ipx=x0; ipx!=x1; ++ipx ) { + for ( int ipy=y0; ipy!=y1; ++ipy) { + ix = (ipx-x0)/multiplicator; + for ( int m=0; m!=multiplicator; ++m ) { + if ( (ipx-x0)%multiplicator==m ) iy=(ipy-y0)*multiplicator + m + shifty; + } + dataMap[iy][ix] = sizeof(header) + (nc_chip * ipy + ipx) * 2; + } + } + + } + + public: + +#ifdef ALDO //VH + using header = jf_header; //VH +#else //VH + using header = sls::defs::sls_receiver_header; +#endif //VH + + jungfrauLGADStrixelsDataSingleChip( const int chip ) + : slsDetectorData( /*nc_strixel*/nc_chip/3, /*nr_strixel*/ nr_chip*5, + nc_chip * nr_chip * 2 + sizeof(header) ) { + std::cout << "Jungfrau strixels single chip" << std::endl; +#ifdef ALDO //VH + std::cout<< "using reduced jf_header" << std::endl; //VH +#endif //VH + + mchip = chip; + //Fill all strixels with dummy values + for (int ix = 0; ix != nc_strixel; ++ix) { + for (int iy = 0; iy != nr_strixel; ++iy) { + dataMap[iy][ix] = sizeof(header); +#ifdef HIGHZ + dataMask[iy][ix] = 0x3fff; +#endif + } + } + + remapGroup(1); + remapGroup(2); + remapGroup(3); + + 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 ((header *)buff)->bunchNumber; //VH +#else //VH + return ((header *)buff)->detHeader.frameNumber; +#endif //VH + }; + + /** + + 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; +#else //VH + return ((header *)buff)->detHeader.packetNumber; +#endif //VH + }; + + + + 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; + + // 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 diff --git a/slsDetectorCalibration/dataStructures/jungfrauStrixelsHalfModuleOldDesign.h b/slsDetectorCalibration/dataStructures/jungfrauStrixelsHalfModuleOldDesign.h new file mode 100644 index 000000000..25de42f80 --- /dev/null +++ b/slsDetectorCalibration/dataStructures/jungfrauStrixelsHalfModuleOldDesign.h @@ -0,0 +1,271 @@ +// SPDX-License-Identifier: LGPL-3.0-or-other +// Copyright (C) 2021 Contributors to the SLS Detector Package +#ifndef JUNGFRAUSTRIXELSHALFMODULEOLD_H +#define JUNGFRAUSTRIXELSHALFMODULEOLD_H +#ifdef CINT +#include "sls/sls_detector_defs_CINT.h" +#else +#include "sls/sls_detector_defs.h" +#endif +#include "slsDetectorData.h" + +/* +/afs/psi.ch/project/mythen/Anna/slsDetectorPackageDeveloperMpc2011/slsDetectorCalibration/jungfrauExecutables +make -f Makefile.rawdataprocess jungfrauRawDataProcessStrx + ../dataStructures/jungfrauLGADStrixelsData.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 +*/ + +namespace strixelsOldDesign { + constexpr int NC_STRIXEL = (1024*3); + constexpr int NR_TOTAL = (512/3); + constexpr int NR_STRIXEL = ( (256-4)/3 ); +} + +typedef struct { + uint64_t bunchNumber; /**< is the frame number */ + uint64_t pre; /**< something */ + +} jf_header; //Aldo's header + + +using namespace strixelsOldDesign; +class jungfrauStrixelsHalfModuleOldDesign : public slsDetectorData { + + private: + int iframe; + +#ifdef ALDO //VH + using header = jf_header; //VH +#else //VH + using header = sls::defs::sls_receiver_header; +#endif //VH + 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 + + */ + jungfrauStrixelsHalfModuleOldDesign() + : slsDetectorData( 1024*3, 512/3, + 512 * 1024 * 2 + sizeof(header) ) { + std::cout << "Jungfrau strixels old design" << std::endl; +#ifdef ALDO //VH + std::cout<< "using reduced jf_header" << std::endl; //VH +#endif //VH + for (int ix = 0; ix != 1024*3; ++ix) { + for (int iy = 0; iy != 512/3; ++iy) { + dataMap[iy][ix] = sizeof(header);//+ ( 1024 * 5 + 300) * 2; //somewhere on the guardring of the LGAD +#ifdef HIGHZ + dataMask[iy][ix] = 0x3fff; +#endif + } + } + + //remap + int ix, iy; + for (int ipx=0; ipx!=1024; ++ipx) { + for (int ipy=0; ipy!=256-4; ++ipy) { + iy=ipy/3; + /* //1 + if (ipy%3==0) ix=ipx*3+2; + if (ipy%3==1) ix=ipx*3+1; + if (ipy%3==2) ix=ipx*3; + */ + /* //2 + if (ipy%3==2) ix=ipx*3+2; + if (ipy%3==0) ix=ipx*3+1; + if (ipy%3==1) ix=ipx*3; + */ + /* //3 + if (ipy%3==1) ix=ipx*3+2; + if (ipy%3==2) ix=ipx*3+1; + if (ipy%3==0) ix=ipx*3; + */ + //4 //This seems to be correct //corresponds to looking from the backside of the sensor + if (ipy%3==0) ix=ipx*3; + if (ipy%3==1) ix=ipx*3+1; + if (ipy%3==2) ix=ipx*3+2; + + /* //5 + if (ipy%3==2) ix=ipx*3; + if (ipy%3==0) ix=ipx*3+1; + if (ipy%3==1) ix=ipx*3+2; + */ + /* //6 + if (ipy%3==1) ix=ipx*3; + if (ipy%3==2) ix=ipx*3+1; + if (ipy%3==0) ix=ipx*3+2; + */ + + if ( ipx!=255 && ipx!=256 && ipx!=511 && ipx!=512 && ipx!=767 && ipx!=768 ) //avoid double pixels + // ( !( ipx%256==0 || ipx%256==255 ) || ipx==0 || ipx==1023 ) + dataMap[iy][ix] = sizeof(header) + (1024 * ipy + ipx) * 2; + // cout << ipx << " " << ipy << " " << ix << " " << iy << endl; + } + } + + 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; + }; + + /* virtual void calcGhost(char *data, int ix, int iy) { */ + /* double val=0; */ + /* ghost[iy][ix]=0; */ + + /* } */ + + /* virtual void calcGhost(char *data) { */ + /* for (int ix=0; ix<25; ix++){ */ + /* for (int iy=0; iy<200; iy++) { */ + /* calcGhost(data, ix,iy); */ + /* } */ + /* } */ + /* // cout << "*" << endl; */ + /* } */ + + /* double getGhost(int ix, int iy) { */ + /* return 0; */ + /* }; */ + + /** + + 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) { +#ifdef ALDO //VH + return ((header *)buff)->bunchNumber; //VH +#else //VH + return ((header *)buff)->detHeader.frameNumber; +#endif //VH + }; + + /** + + 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; +#else //VH + return ((header *)buff)->detHeader.packetNumber; +#endif //VH + }; + + + + 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; + + // 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 diff --git a/slsDetectorCalibration/dataStructures/slsDetectorData.h b/slsDetectorCalibration/dataStructures/slsDetectorData.h index f7b3f1a5b..adf26ae4a 100644 --- a/slsDetectorCalibration/dataStructures/slsDetectorData.h +++ b/slsDetectorCalibration/dataStructures/slsDetectorData.h @@ -217,6 +217,7 @@ void slsDetectorData::setDataMap(int **dMap) { } } } + /* //commented this part because it causes out-of-bound issues if nx or ny are larger than dataMap bounds (single-chip readout of strixel with groups of different pitches) VH 2023-02-24 for (iy = 0; iy < ny; iy++) { for (ix = 0; ix < nx; ix++) { ip = dataMap[iy][ix] / sizeof(dataType); @@ -224,7 +225,7 @@ void slsDetectorData::setDataMap(int **dMap) { ymap[ip] = iy; } } - + */ // cout << "nx:" < #include -using namespace std; +//using namespace std; //#ifdef MYROOT1 //: public TObject @@ -128,7 +128,7 @@ class slsInterpolation { nSubPixelsY * nPixelsY); delete[] gm; } else - cout << "Could not allocate float image " << endl; + std::cout << "Could not allocate float image " << std::endl; return NULL; } @@ -323,7 +323,7 @@ class slsInterpolation { if (ix > 1) sumR += cl[ix + iy * 3]; if (iy < 1) - sumB = cl[ix + iy * 3]; + sumB = cl[ix + iy * 3]; //???? not "+="? VH if (iy > 1) sumT += cl[ix + iy * 3]; } @@ -472,13 +472,13 @@ class slsInterpolation { val = cl[ix + 3 * iy]; sum += val; if (iy == 0) - l += val; - if (iy == 2) - r += val; - if (ix == 0) b += val; - if (ix == 2) + if (iy == 2) t += val; + if (ix == 0) + l += val; + if (ix == 2) + r += val; } } if (sum > 0) { @@ -526,6 +526,185 @@ class slsInterpolation { return calcEta3X(cli, etax, etay, sum); } + /************************************************/ + /* Additional strixel eta functions by Viktoria */ + /************************************************/ + //Etax: only central row, etay: only central column + static int calcEta1x3( double* cl, double& etax, double& etay, double& toth, double& totv ) { + double l, r, t, b; + //sum = cl[0] + cl[1] + cl[2] + cl[3] + cl[4] + cl[5] + cl[6] + cl[7] + cl[8]; + toth = cl[3] + cl[4] + cl[5]; + if (toth > 0) { + l = cl[3]; + r = cl[5]; + } + etax = (-l + r) / toth; + totv = cl[1] + cl[4] + cl[7]; + if (toth > 0) { + b = cl[1]; + t = cl[7]; + } + etay = (-b + t) / totv; + return -1; + } + + static int calcEta1x3( int* cl, double& etax, double& etay, double& toth, double& totv ) { + double cli[9]; + for ( int ix = 0; ix != 9; ++ix ) + cli[ix] = cl[ix]; + return calcEta1x3( cli, etax, etay, toth , totv ); + } + + //Eta 1x2 essentially the same as etaL, but we also return toth and totv + static int calcEta1x2(double totquad, int corner, double sDum[2][2], + double &etax, double &etay, double& toth, double& totv) { + double t, r; + if (totquad > 0) { + switch (corner) { + case TOP_LEFT: + t = sDum[1][1]; + r = sDum[0][1]; + toth = sDum[0][1] + sDum[0][0]; + totv = sDum[0][1] + sDum[1][1]; + break; + case TOP_RIGHT: + t = sDum[1][0]; + r = sDum[0][1]; + toth = sDum[0][1] + sDum[0][0]; + totv = sDum[1][0] + sDum[0][0]; + break; + case BOTTOM_LEFT: + r = sDum[1][1]; + t = sDum[1][1]; + toth = sDum[1][0] + sDum[1][1]; + totv = sDum[0][1] + sDum[1][1]; + break; + case BOTTOM_RIGHT: + t = sDum[1][0]; + r = sDum[1][1]; + toth = sDum[1][0] + sDum[1][1]; + totv = sDum[1][0] + sDum[0][0]; + break; + default: + etax = -1000; + etay = -1000; + return 0; + } + // etax=r/totquad; + // etay=t/totquad; + etax = r / toth; + etay = t / totv; + } + return 0; + } + + static int calcEta1x2( double *cl, double &etax, double &etay, double &sum, + double &totquad, double sDum[2][2], double& toth, double& totv ) { + int corner = calcQuad( cl, sum, totquad, sDum ); + calcEta1x2( totquad, corner, sDum, etax, etay, toth, totv ); + return corner; + } + + static int calcEta1x2( int *cl, double &etax, double &etay, double &sum, + double &totquad, double sDum[2][2], double& toth, double& totv ) { + int corner = calcQuad( cl, sum, totquad, sDum ); + calcEta1x2( totquad, corner, sDum, etax, etay, toth , totv ); + return corner; + } + + //Two functions to calculate 2x3 or 3x2 eta + static int calcEta2x3( double* cl, double totquad, int corner, double& etax, double& etay, double& tot6 ) { + double t, b, r; + if (totquad > 0) { + switch (corner) { + case TOP_LEFT: + case BOTTOM_LEFT: + t = cl[6] + cl[7]; + b = cl[0] + cl[1]; + r = cl[1] + cl[4] + cl[7]; + tot6 = cl[0] + cl[1] + cl[3] + cl[4] + cl[6] + cl[7]; + break; + case TOP_RIGHT: + case BOTTOM_RIGHT: + t = cl[7] + cl[8]; + b = cl[1] + cl[2]; + r = cl[2] + cl[5] + cl[8]; + tot6 = cl[1] + cl[2] + cl[4] + cl[5] + cl[7] + cl[8]; + break; + default: + etax = -1000; + etay = -1000; + return -1; + } + etax = r / tot6; + etay = (-b + t) / tot6; + } + return -1; + } + + static int calcEta3x2( double* cl, double totquad, int corner, double& etax, double& etay, double& tot6 ) { + double t, l, r; + if (totquad > 0) { + switch (corner) { + case TOP_LEFT: + case TOP_RIGHT: + l = cl[3] + cl[6]; + r = cl[5] + cl[8]; + t = cl[6] + cl[7] + cl[8]; + tot6 = cl[3] + cl[4] + cl[5] + cl[6] + cl[7] + cl[8]; + break; + case BOTTOM_LEFT: + case BOTTOM_RIGHT: + l = cl[0] + cl[3]; + r = cl[2] + cl[5]; + t = cl[3] + cl[4] + cl[5]; + tot6 = cl[0] + cl[1] + cl[2] + cl[3] + cl[4] + cl[5]; + break; + default: + etax = -1000; + etay = -1000; + return -1; + } + etax = (-l + r) / tot6; + etay = t / tot6; + } + return -1; + } + + //overload including both eta2x3 and eta3x2 + //ornt (orientation of long side (3) of eta) decides which eta is chosen + enum orientation { + HORIZONTAL_ORIENTATION = 0, + VERTICAL_ORIENTATION = 1, + UNDEFINED_ORIENTATION = -1 + }; + static int calcEta2x3( int ornt, double* cl, double& etax, double& etay, double& tot6 ) { + double sum{}; + double totquad{}; + double sDum[2][2]{}; + int corner = calcQuad( cl, sum, totquad, sDum ); + switch (ornt) { + case HORIZONTAL_ORIENTATION: + calcEta3x2( cl, totquad, corner, etax, etay, tot6 ); + break; + case VERTICAL_ORIENTATION: + calcEta2x3( cl, totquad, corner, etax, etay, tot6 ); + break; + default: + etax = -1000; + etay = -1000; + return -1; + } + return corner; + } + + static int calcEta2x3( int strxo, int* cl, double& etax, double& etay, double& tot6 ) { + double cli[9]{}; + for ( int ix = 0; ix != 9; ++ix ) + cli[ix] = cl[ix]; + return calcEta2x3( strxo, cli, etax, etay, tot6 ); + } + /* static int calcMyEta(double totquad, int quad, double *cl, double &etax, * double &etay) { */ /* double l,r,t,b, sum; */ diff --git a/slsDetectorCalibration/jungfrauExecutables/Makefile.rawdataprocess b/slsDetectorCalibration/jungfrauExecutables/Makefile.rawdataprocess new file mode 100644 index 000000000..8f16632d1 --- /dev/null +++ b/slsDetectorCalibration/jungfrauExecutables/Makefile.rawdataprocess @@ -0,0 +1,76 @@ +# SPDX-License-Identifier: LGPL-3.0-or-other +# Copyright (C) 2021 Contributors to the SLS Detector Package +#module add CBFlib/0.9.5 +INCDIR=-I. -I../ -I../interpolations -I../interpolations/etaVEL -I../dataStructures -I../../slsSupportLib/include/ -I../../slsReceiverSoftware/include/ -I../tiffio/include + +LDFLAG= ../tiffio/src/tiffIO.cpp -L/usr/lib64/ -lpthread -lm -lstdc++ -pthread -lrt -ltiff -O3 -std=c++11 + +MAIN=jungfrauClusterFinder.cpp + + +all: jungfrauRawDataProcess + +jungfrauRawDataProcess: jungfrauRawDataProcess.cpp $(INCS) clean + g++ -o jungfrauRawDataProcess jungfrauRawDataProcess.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL -DMODULE + +jungfrauRawDataProcessStrx: jungfrauRawDataProcess.cpp $(INCS) clean + g++ -o jungfrauRawDataProcessStrx jungfrauRawDataProcess.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL -DJFSTRX + +jungfrauRawDataProcessStrxChip1: jungfrauRawDataProcess.cpp $(INCS) clean + g++ -o jungfrauRawDataProcessStrxChip1 jungfrauRawDataProcess.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL -DJFSTRXCHIP1 + +jungfrauRawDataProcessStrxChip6: jungfrauRawDataProcess.cpp $(INCS) clean + g++ -o jungfrauRawDataProcessStrxChip6 jungfrauRawDataProcess.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL -DJFSTRXCHIP6 + +jungfrauRawDataProcessStrxChip1Aldo: jungfrauRawDataProcess.cpp $(INCS) clean + g++ -o jungfrauRawDataProcessStrxChip1Aldo jungfrauRawDataProcess.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL -DJFSTRXCHIP1 -DALDO + +jungfrauRawDataProcessStrxChip6Aldo: jungfrauRawDataProcess.cpp $(INCS) clean + g++ -o jungfrauRawDataProcessStrxChip6Aldo jungfrauRawDataProcess.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL -DJFSTRXCHIP6 -DALDO + +jungfrauRawDataProcessStrxAldo: jungfrauRawDataProcess.cpp $(INCS) clean + g++ -o jungfrauRawDataProcessStrxAldo jungfrauRawDataProcess.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL -DJFSTRX -DALDO + +jungfrauRawDataProcessStrxOld: jungfrauRawDataProcess.cpp $(INCS) clean + g++ -o jungfrauRawDataProcessStrxOld jungfrauRawDataProcess.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL -DJFSTRXOLD + +jungfrauRawDataProcessStrxOldAldo: jungfrauRawDataProcess.cpp $(INCS) clean + g++ -o jungfrauRawDataProcessStrxOldAldo jungfrauRawDataProcess.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL -DJFSTRXOLD -DALDO + + + +jungfrauClusterFinder: jungfrauClusterFinder.cpp $(INCS) clean + g++ -o jungfrauClusterFinder jungfrauClusterFinder.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL + + +jungfrauClusterFinderHighZ: jungfrauClusterFinder.cpp $(INCS) clean + g++ -o jungfrauClusterFinderHighZ jungfrauClusterFinder.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL -DHIGHZ + + + + +jungfrauMakeEta: jungfrauInterpolation.cpp $(INCS) clean + g++ -o jungfrauMakeEta jungfrauInterpolation.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DFF + +jungfrauInterpolation: jungfrauInterpolation.cpp $(INCS) clean + g++ -o jungfrauInterpolation jungfrauInterpolation.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) + +jungfrauNoInterpolation: jungfrauNoInterpolation.cpp $(INCS) clean + g++ -o jungfrauNoInterpolation jungfrauNoInterpolation.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) + +jungfrauPhotonCounter: jungfrauPhotonCounter.cpp $(INCS) clean + g++ -o jungfrauPhotonCounter jungfrauPhotonCounter.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DNEWRECEIVER + +jungfrauAnalog: jungfrauPhotonCounter.cpp $(INCS) clean + g++ -o jungfrauAnalog jungfrauPhotonCounter.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DNEWRECEIVER -DANALOG + +jungfrauPhotonCounterHighZ: jungfrauPhotonCounter.cpp $(INCS) clean + g++ -o jungfrauPhotonCounterHighZ jungfrauPhotonCounter.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DNEWRECEIVER -DHIGHZ + +jungfrauAnalogHighZ: jungfrauPhotonCounter.cpp $(INCS) clean + g++ -o jungfrauAnalogHighZ jungfrauPhotonCounter.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DNEWRECEIVER -DANALOG -DHIGHZ + +clean: + rm -f jungfrauClusterFinder jungfrauMakeEta jungfrauInterpolation jungfrauNoInterpolation jungfrauPhotonCounter jungfrauAnalog + + diff --git a/slsDetectorCalibration/jungfrauExecutables/jungfrauRawDataProcess.cpp b/slsDetectorCalibration/jungfrauExecutables/jungfrauRawDataProcess.cpp new file mode 100644 index 000000000..a0bf7c5cc --- /dev/null +++ b/slsDetectorCalibration/jungfrauExecutables/jungfrauRawDataProcess.cpp @@ -0,0 +1,429 @@ +// 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 JFSTRXOLD && !defined JFSTRXCHIP1 && !defined JFSTRXCHIP6 +#ifndef MODULE +#include "jungfrauHighZSingleChipData.h" +#endif +#ifdef MODULE +#include "jungfrauModuleData.h" +#endif +#endif + +#ifdef JFSTRX +#include "jungfrauLGADStrixelsData.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 +using namespace std; + +int main(int argc, char *argv[]) { + + if (argc < 5) { + cout << "Usage is " << argv[0] + << "indir outdir fname(no extension) fextension [runmin] [runmax] [pedfile (raw or tiff)] [threshold] " + "[nframes] [xmin xmax ymin ymax] [gainmap]" + << endl; + cout << "threshold <0 means analog; threshold=0 means cluster finder; " + "threshold>0 means photon counting" + << endl; + cout << "nframes <0 means sum everything; nframes=0 means one file per " + "run; nframes>0 means one file every nframes" + << endl; + return 1; + } + + int fifosize = 1000; + int nthreads = 10; + int csize = 3; //3 + int nsigma = 5; + int nped = 10000; + + int cf = 0; + + +#if !defined JFSTRX && !defined JFSTRXOLD && !defined JFSTRXCHIP1 && !defined JFSTRXCHIP6 +#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 JFSTRX + cout << "Jungfrau strixel full module readout" << endl; + jungfrauLGADStrixelsData *decoder = new jungfrauLGADStrixelsData(); + int nx = 1024/5, ny = 512*5; +#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); + + + + + + cout << "Detector size is " << nx << " " << ny << endl; + + + + + + + + + double *gainmap = NULL; + //float *gm; + + int ff, np; + // cout << " data size is " << dsize; + + ifstream filebin; + char *indir = argv[1]; + char *outdir = argv[2]; + char *fformat = argv[3]; + char *fext = argv[4]; + int runmin = 0; + + // cout << "argc is " << argc << endl; + if (argc >= 6) { + runmin = atoi(argv[5]); + } + + int runmax = runmin; + + if (argc >= 7) { + runmax = atoi(argv[6]); + } + + char *pedfile = NULL; + if (argc >= 8) { + pedfile = argv[7]; + } + double thr = 0; + double thr1 = 1; + + if (argc >= 9) { + thr = atof(argv[8]); + } + + int nframes = 0; + + if (argc >= 10) { + nframes = atoi(argv[9]); + } + + 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]); + } + + char *gainfname = NULL; + if (argc > 14) { + gainfname = argv[14]; + cout << "Gain map file name is: " << gainfname << endl; + } + + char ffname[10000]; + char fname[10000]; + char imgfname[10000]; + char cfname[10000]; + + std::time_t end_time; + + FILE *of = NULL; + cout << "input directory is " << indir << endl; + cout << "output directory is " << outdir << endl; + cout << "input file is " << fformat << endl; + cout << "runmin is " << runmin << endl; + cout << "runmax is " << runmax << endl; + if (pedfile) + cout << "pedestal file is " << pedfile << endl; + if (thr > 0) + cout << "threshold is " << thr << endl; + cout << "Nframes is " << nframes << 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)) + cout << "using gain map " << gainfname << endl; + else + cout << "Could not open gain map " << gainfname << endl; + } else + thr = 0.15 * thr; + filter->newDataSet(); + //int dsize = decoder->getDataSize(); + + if (thr > 0) { + cout << "threshold is " << thr << endl; + filter->setThreshold(thr); + cf = 0; + + } else + cf = 1; + + filter->setROI(xmin, xmax, ymin, ymax); + std::time(&end_time); + cout << std::ctime(&end_time) << endl; + + char *buff; + + // multiThreadedAnalogDetector *mt=new + // multiThreadedAnalogDetector(filter,nthreads,fifosize); + multiThreadedCountingDetector *mt = + new multiThreadedCountingDetector(filter, nthreads, fifosize); + mt->setClusterSize(csize,csize); +#ifndef ANALOG + mt->setDetectorMode(ePhotonCounting); + cout << "Counting!" << endl; + if (thr > 0) { + cf = 0; + } +#endif +//{ +#ifdef ANALOG + mt->setDetectorMode(eAnalog); + cout << "Analog!" << endl; + cf = 0; + // thr1=thr; +#endif + // } + + mt->StartThreads(); + mt->popFree(buff); + + // cout << "mt " << endl; + + int ifr = 0; + + char froot[1000]; + double *ped=new double[nx * ny];//, *ped1; + + int pos,pos1; + + if (pedfile) { + + if (string(pedfile).find(".dat") != std::string::npos) { + pos1=string(pedfile).rfind("/"); + strcpy(froot,pedfile+pos1); + pos=string(froot).find(".dat"); + froot[pos]='\0'; + } + + cout << "PEDESTAL " << endl; + // sprintf(imgfname, "%s/pedestals.tiff", outdir); + + if (string(pedfile).find(".tif") == std::string::npos) { + sprintf(fname, "%s", pedfile); + cout << fname << endl; + std::time(&end_time); + //cout << "aaa" << std::ctime(&end_time) << endl; + + mt->setFrameMode(ePedestal); + // sprintf(fn,fformat,irun); + filebin.open((const char *)(fname), ios::in | ios::binary); + // //open file + if (filebin.is_open()) { + ff = -1; + while (decoder->readNextFrame(filebin, ff, np, buff)) { + // if (np == 40) { + if ((ifr+1) % 100 == 0) { + cout << " ****" << decoder->getValue(buff,20,20);// << endl; + } + mt->pushData(buff); + mt->nextThread(); + mt->popFree(buff); + ifr++; + if (ifr % 100 == 0) { + cout << " ****" << ifr << " " << ff << " " << np << endl; + } //else + //cout << ifr << " " << ff << " " << np << endl; + if (ifr>=1000) + break; + ff = -1; + } + filebin.close(); + while (mt->isBusy()) { + ; + } + + sprintf(imgfname, "%s/%s_ped.tiff", outdir, froot); + mt->writePedestal(imgfname); + sprintf(imgfname, "%s/%s_rms.tiff", outdir, froot); + mt->writePedestalRMS(imgfname); + + } else + cout << "Could not open pedestal file " << fname + << " for reading " << endl; + } else { + float *pp = ReadFromTiff(pedfile, 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); + cout << "Pedestal set from tiff file " << pedfile << endl; + } else { + cout << "Could not open pedestal tiff file " << pedfile + << " for reading " << endl; + } + } + std::time(&end_time); + cout << std::ctime(&end_time) << endl; + } + + ifr = 0; + int ifile = 0; + + mt->setFrameMode(eFrame); + + for (int irun = runmin; irun <= runmax; irun++) { + cout << "DATA "; + // sprintf(fn,fformat,irun); + sprintf(ffname, "%s/%s.%s", indir, fformat, fext); + sprintf(fname, (const char*)ffname, irun); + sprintf(ffname, "%s/%s.tiff", outdir, fformat); + sprintf(imgfname, (const char*)ffname, irun); + sprintf(ffname, "%s/%s.clust", outdir, fformat); + sprintf(cfname, (const char*)ffname, irun); + cout << fname << " "; + cout << imgfname << endl; + std::time(&end_time); + cout << std::ctime(&end_time) << endl; + // cout << fname << " " << outfname << " " << imgfname << endl; + filebin.open((const char *)(fname), ios::in | ios::binary); + // //open file + ifile = 0; + if (filebin.is_open()) { + if (thr <= 0 && cf != 0) { // cluster finder + if (of == NULL) { + of = fopen(cfname, "w"); + if (of) { + mt->setFilePointer(of); + cout << "file pointer set " << endl; + } else { + cout << "Could not open " << cfname << " for writing " + << 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) { + cout << " ****" << decoder->getValue(buff,20,20);// << endl; + } + mt->pushData(buff); + // // //pop + mt->nextThread(); + mt->popFree(buff); + + ifr++; + if (ifr % 100 == 0) + cout << " " << ifr << " " << ff << endl; + if (nframes > 0) { + if (ifr % nframes == 0) { + sprintf(ffname, "%s/%s_f%05d.tiff", outdir, fformat, + ifile); + sprintf(imgfname, (const char*)ffname, irun); + mt->writeImage(imgfname, thr1); + mt->clearImage(); + ifile++; + } + } + // } else + // cout << ifr << " " << ff << " " << np << endl; + ff = -1; + } + cout << "--" << endl; + filebin.close(); + while (mt->isBusy()) { + ; + } + if (nframes >= 0) { + if (nframes > 0) { + sprintf(ffname, "%s/%s_f%05d.tiff", outdir, fformat, ifile); + sprintf(imgfname, (const char*)ffname, irun); + } else { + sprintf(ffname, "%s/%s.tiff", outdir, fformat); + sprintf(imgfname, (const char*)ffname, irun); + } + cout << "Writing tiff to " << imgfname << " " << thr1 << endl; + mt->writeImage(imgfname, thr1); + mt->clearImage(); + if (of) { + fclose(of); + of = NULL; + mt->setFilePointer(NULL); + } + } + std::time(&end_time); + cout << std::ctime(&end_time) << endl; + } else + cout << "Could not open " << fname << " for reading " << endl; + } + if (nframes < 0) { + sprintf(ffname, "%s/%s.tiff", outdir, fformat); + strcpy(imgfname, ffname); + cout << "Writing tiff to " << imgfname << " " << thr1 << endl; + mt->writeImage(imgfname, thr1); + } + + return 0; +}