diff --git a/libs/pybind11 b/libs/pybind11 index 8de7772cc..4f72ef846 160000 --- a/libs/pybind11 +++ b/libs/pybind11 @@ -1 +1 @@ -Subproject commit 8de7772cc72daca8e947b79b83fea46214931604 +Subproject commit 4f72ef846fe8453596230ac285eeaa0ce3278bb4 diff --git a/slsDetectorCalibration/dataStructures/jungfrauHighZSingleChipData.h b/slsDetectorCalibration/dataStructures/jungfrauHighZSingleChipData.h new file mode 100644 index 000000000..d3bc61f62 --- /dev/null +++ b/slsDetectorCalibration/dataStructures/jungfrauHighZSingleChipData.h @@ -0,0 +1,238 @@ +#ifndef JUNGFRAUHIGHZSINGLECHIPDATA_H +#define JUNGFRAUHIGHZSINGLECHIPDATA_H +#include "slsDetectorData.h" + +//#define VERSION_V2 + /** + @short structure for a Detector Packet or Image Header + @li frameNumber is the frame number + @li expLength is the subframe number (32 bit eiger) or real time exposure time in 100ns (others) + @li packetNumber is the packet number + @li bunchId is the bunch id from beamline + @li timestamp is the time stamp with 10 MHz clock + @li modId is the unique module id (unique even for left, right, top, bottom) + @li xCoord is the x coordinate in the complete detector system + @li yCoord is the y coordinate in the complete detector system + @li zCoord is the z coordinate in the complete detector system + @li debug is for debugging purposes + @li roundRNumber is the round robin set number + @li detType is the detector type see :: detectorType + @li version is the version number of this structure format + */ + typedef struct { + uint64_t bunchNumber; /**< is the frame number */ + uint64_t pre; /**< something */ + + } jf_header; + + + + +class jungfrauHighZSingleChipData : public slsDetectorData { + + private: + + int iframe; + public: + /** + Implements the slsReceiverData structure for the moench02 prototype read out by a module i.e. using the slsReceiver + (160x160 pixels, 40 packets 1286 large etc.) + \param c crosstalk parameter for the output buffer + + */ + jungfrauHighZSingleChipData(): slsDetectorData(256, 256, 256*256*2+sizeof(jf_header)) { + + + for (int ix=0; ix<256; ix++) { + for (int iy=0; iy<256; iy++) { + dataMap[iy][ix]=sizeof(jf_header)+(256*iy+ix)*2; +#ifdef HIGHZ + dataMask[iy][ix]=0x3fff; +#endif + } + } + + iframe=0; + // cout << "data struct created" << endl; + }; + + + + /** + Returns the value of the selected channel for the given dataset as double. + \param data pointer to the dataset (including headers etc) + \param ix pixel number in the x direction + \param iy pixel number in the y direction + \returns data for the selected channel, with inversion if required as double + + */ + virtual double getValue(char *data, int ix, int iy=0) { + + uint16_t val=getChannel(data, ix, iy)&0x3fff; + 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){return ((jf_header*)buff)->bunchNumber;};//*((int*)(buff+5))&0xffffff;}; + + /** + + Returns the packet number for the given dataset. purely virtual func + \param buff pointer to the dataset + \returns packet number number + + + + */ + int getPacketNumber(char *buff){return 0;}//((*(((int*)(buff+4))))&0xff)+1;}; + +/* /\** */ + +/* 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){ndata=dsize; setDataSize(dsize); return data;}; */ + + +/* /\** */ + +/* Loops over a file stream until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). Can be overloaded for different kind of detectors! */ +/* \param filebin input file stream (binary) */ +/* \returns pointer to the begin of the last good frame, NULL if no frame is found or last frame is incomplete */ + +/* *\/ */ +/* virtual char *readNextFrame(ifstream &filebin){ */ +/* // int afifo_length=0; */ +/* uint16_t *afifo_cont; */ +/* int ib=0; */ +/* if (filebin.is_open()) { */ +/* afifo_cont=new uint16_t[dataSize/2]; */ +/* while (filebin.read(((char*)afifo_cont)+ib,2)) { */ +/* ib+=2; */ +/* if (ib==dataSize) break; */ +/* } */ +/* if (ib>0) { */ +/* iframe++; */ +/* // cout << ib << "-" << endl; */ +/* return (char*)afifo_cont; */ +/* } else { */ +/* delete [] afifo_cont; */ +/* return NULL; */ +/* } */ +/* } */ +/* return NULL; */ +/* }; */ + + + virtual char *readNextFrame(ifstream &filebin) { + int ff=-1, np=-1; + return readNextFrame(filebin, ff, np); + }; + + virtual char *readNextFrame(ifstream &filebin, int &ff) { + int np=-1; + return readNextFrame(filebin, ff, np); + }; + + virtual char *readNextFrame(ifstream &filebin, int& ff, int &np) { + char *data=new char[dataSize]; + char *d=readNextFrame(filebin, ff, np, data); + if (d==NULL) {delete [] data; data=NULL;} + return data; + } + + + + + virtual char *readNextFrame(ifstream &filebin, int& ff, int &np, char *data) { + char *retval=0; + int nd; + int fnum = -1; + np=0; + int pn; + + // cout << dataSize << endl; + if (ff>=0) + fnum=ff; + + if (filebin.is_open()) { + if (filebin.read(data, dataSize) ){ + ff=getFrameNumber(data); + np=getPacketNumber(data); + return data; + } + } + return NULL; + }; + + + + /** + + Loops over a memory slot until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). purely virtual func + \param data pointer to the memory to be analyzed + \param ndata reference to the amount of data found for the frame, in case the frame is incomplete at the end of the memory slot + \param dsize size of the memory slot to be analyzed + \returns pointer to the beginning of the last good frame (might be incomplete if ndata smaller than dataSize), or NULL if no frame is found + + */ + virtual char *findNextFrame(char *data, int &ndata, int dsize){ + if (dsize +#include +#include +#include +#endif + +#include + +#define max(a,b) ((a) > (b) ? (a) : (b)) +#define min(a,b) ((a) < (b) ? (a) : (b)) +#define ELEM_SWAP(a,b) { register int t=(a);(a)=(b);(b)=t; } + + +using namespace std; + +#ifdef MYROOT + +Double_t energyCalibrationFunctions::pedestal(Double_t *x, Double_t *par) { + return par[0]-par[1]*sign*x[0]; +} + + +Double_t energyCalibrationFunctions::gaussChargeSharing(Double_t *x, Double_t *par) { + Double_t f, arg=0; + // Gaussian exponent + if (par[3]!=0) { + arg=sign*(x[0]-par[2])/par[3]; + } + // the Gaussian + f=TMath::Exp(-1*arg*arg/2.); + // Gaussian + error function + f=f+par[5]/2.*(TMath::Erfc(arg/(TMath::Sqrt(2.)))); + // Gaussian + error function + pedestal + return par[4]*f+pedestal(x,par); +} + +Double_t energyCalibrationFunctions::gaussChargeSharingKb(Double_t *x, Double_t *par) { + Double_t f, arg=0,argb=0; + // Gaussian exponent + if (par[3]!=0) { + arg=sign*(x[0]-par[2])/par[3]; + argb=sign*(x[0]-(par[6]*par[2]))/par[3]; // using absolute kb mean might seem better but like this the ratio can be fixed + } + // the Gaussian + f=TMath::Exp(-1*arg*arg/2.); + f=f+par[7]*(TMath::Exp(-1*argb*argb/2.)); + // Gaussian + error function + f=f+par[5]/2.*(TMath::Erfc(arg/(TMath::Sqrt(2.)))); + f=f+par[7]*par[5]/2.*(TMath::Erfc(argb/(TMath::Sqrt(2.)))); + // Gaussian + error function + pedestal + return par[4]*f+pedestal(x,par); +} + +Double_t energyCalibrationFunctions::gaussChargeSharingKaDoublet(Double_t *x, Double_t *par) { + Double_t f, f2, arg=0, arg2=0; + // Gaussian exponent + if (par[3]!=0) { + arg=sign*(x[0]-par[2])/par[3]; + arg2=sign*(x[0]-par[6])/par[3]; + } + // the Gaussian + f=TMath::Exp(-1*arg*arg/2.); + f2=TMath::Exp(-1*arg2*arg2/2.); + // Gaussian + error function + f=f+par[5]/2.*(TMath::Erfc(arg/(TMath::Sqrt(2.)))); + f2=f2+par[5]/2.*(TMath::Erfc(arg/(TMath::Sqrt(2.)))); // shouldn't this be arg2? + // Gaussian + error function + pedestal + return par[4]*f+par[7]*f2+pedestal(x,par); +} + +Double_t energyCalibrationFunctions::gaussChargeSharingPixel(Double_t *x, Double_t *par) { + Double_t f; + if (par[3]<=0 || par[2]*(*x)<=0 || par[5]<0 || par[4]<=0) return 0; + + Double_t pp[3]; + + pp[0]=0; + pp[1]=par[2]; + pp[2]=par[3]; + + + f=(par[5]-par[6]*(TMath::Log(*x/par[2])))*erfBox(x,pp); + f+=par[4]*TMath::Gaus(*x, par[2], par[3], kTRUE); + return f+pedestal(x,par); +} + +Double_t energyCalibrationFunctions::erfBox(Double_t *z, Double_t *par) { + + + + Double_t m=par[0]; + Double_t M=par[1]; + + if (par[0]>par[1]) { + m=par[1]; + M=par[0]; + } + + if (m==M) + return 0; + + + if (par[2]<=0) { + if (*z>=m && *z<=M) + return 1./(M-m); + else + return 0; + + } + + return (TMath::Erfc((z[0]-M)/par[2])-TMath::Erfc((z[0]-m)/par[2]))*0.5/(M-m); + +} + + +// basic erf function +Double_t energyCalibrationFunctions::erfFunction(Double_t *x, Double_t *par) { + double arg=0; + if (par[1]!=0) arg=(par[0]-x[0])/par[1]; + return ((par[2]/2.*(1+TMath::Erf(sign*arg/(TMath::Sqrt(2)))))); +}; + + +Double_t energyCalibrationFunctions::erfFunctionChargeSharing(Double_t *x, Double_t *par) { + Double_t f; + + f=erfFunction(x, par+2)*(1+par[5]*(par[2]-x[0]))+par[0]-par[1]*x[0]*sign; + return f; +}; + + +Double_t energyCalibrationFunctions::erfFuncFluo(Double_t *x, Double_t *par) { + Double_t f; + f=erfFunctionChargeSharing(x, par)+erfFunction(x, par+6)*(1+par[9]*(par[6]-x[0])); + return f; +}; +#endif + +double energyCalibrationFunctions::median(double *x, int n){ + // sorts x into xmed array and returns median + // n is number of values already in the xmed array + double xmed[n]; + int k,i,j; + + for (i=0; i*(x+j)) + k++; + if (*(x+i)==*(x+j)) { + if (i>j) + k++; + } + } + xmed[k]=*(x+i); + } + k=n/2; + return xmed[k]; +} + + +int energyCalibrationFunctions::quick_select(int arr[], int n){ + int low, high ; + int median; + int middle, ll, hh; + + low = 0 ; high = n-1 ; median = (low + high) / 2; + for (;;) { + if (high <= low) /* One element only */ + return arr[median] ; + + if (high == low + 1) { /* Two elements only */ + if (arr[low] > arr[high]) + ELEM_SWAP(arr[low], arr[high]) ; + return arr[median] ; + } + + /* Find median of low, middle and high items; swap into position low */ + middle = (low + high) / 2; + if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]) ; + if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ; + if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]) ; + + /* Swap low item (now in position middle) into position (low+1) */ + ELEM_SWAP(arr[middle], arr[low+1]) ; + + /* Nibble from each end towards middle, swapping items when stuck */ + ll = low + 1; + hh = high; + for (;;) { + do ll++; while (arr[low] > arr[ll]) ; + do hh--; while (arr[hh] > arr[low]) ; + + if (hh < ll) + break; + + ELEM_SWAP(arr[ll], arr[hh]) ; + } + + /* Swap middle item (in position low) back into correct position */ + ELEM_SWAP(arr[low], arr[hh]) ; + + /* Re-set active partition */ + if (hh <= median) + low = ll; + if (hh >= median) + high = hh - 1; + } +} + +int energyCalibrationFunctions::kth_smallest(int *a, int n, int k){ + register int i,j,l,m ; + register double x ; + + l=0 ; m=n-1 ; + while (lSetParNames("Background Offset","Background Slope","Inflection Point","Noise RMS", "Number of Photons","Charge Sharing Slope"); + + fspectrum=new TF1("fspectrum",funcs,&energyCalibrationFunctions::spectrum,0,1000,6,"energyCalibrationFunctions","spectrum"); + fspectrum->SetParNames("Background Pedestal","Background slope", "Peak position","Noise RMS", "Number of Photons","Charge Sharing Pedestal"); + fspectrumkb=new TF1("fspectrumkb",funcs,&energyCalibrationFunctions::spectrumkb,0,1000,8,"energyCalibrationFunctions","spectrumkb"); + fspectrumkb->SetParNames("Background Pedestal","Background slope", "Peak position","Noise RMS", "Number of Photons","Charge Sharing Pedestal","kb mean","kb frac"); + + fspectrumkadoublet=new TF1("fspectrumkadoublet",funcs,&energyCalibrationFunctions::spectrumkadoublet,0,1000,8,"energyCalibrationFunctions","spectrumkadoublet"); + fspectrumkadoublet->SetParNames("Background Pedestal","Background slope", "Peak position","Noise RMS", "Number of Photons","Charge Sharing Pedestal","ka2 mean","n2"); + + fspixel=new TF1("fspixel",funcs,&energyCalibrationFunctions::spectrumPixel,0,1000,7,"energyCalibrationFunctions","spectrumPixel"); + fspixel->SetParNames("Background Pedestal","Background slope", "Peak position","Noise RMS", "Number of Photons","Charge Sharing Pedestal","Corner"); + +#endif + + +} + + + +void energyCalibration::fixParameter(int ip, Double_t val){ + + fscurve->FixParameter(ip, val); + fspectrum->FixParameter(ip, val); + fspectrumkb->FixParameter(ip, val); + fspectrumkadoublet->FixParameter(ip, val); +} + + +void energyCalibration::releaseParameter(int ip){ + + fscurve->ReleaseParameter(ip); + fspectrum->ReleaseParameter(ip); + fspectrumkb->ReleaseParameter(ip); + fspectrumkadoublet->ReleaseParameter(ip); +} + + + + + + + +energyCalibration::~energyCalibration(){ +#ifdef MYROOT + delete fscurve; + delete fspectrum; + delete fspectrumkb; + delete fspectrumkadoublet; +#endif + +} + +#ifdef MYROOT + + + + + +TH1F* energyCalibration::createMedianHistogram(TH2F* h2, int ch0, int nch, int direction) { + + if (h2==NULL || nch==0) + return NULL; + + double *x=new double[nch]; + TH1F *h1=NULL; + + double val=-1; + + if (direction==0) { + h1=new TH1F("median","Median",h2->GetYaxis()->GetNbins(),h2->GetYaxis()->GetXmin(),h2->GetYaxis()->GetXmax()); + for (int ib=0; ibGetXaxis()->GetNbins(); ib++) { + for (int ich=0; ichGetBinContent(ch0+ich+1,ib+1); + } + val=energyCalibrationFunctions::median(x, nch); + h1->SetBinContent(ib+1,val); + } + } else if (direction==1) { + h1=new TH1F("median","Median",h2->GetXaxis()->GetNbins(),h2->GetXaxis()->GetXmin(),h2->GetXaxis()->GetXmax()); + for (int ib=0; ibGetYaxis()->GetNbins(); ib++) { + for (int ich=0; ichGetBinContent(ib+1,ch0+ich+1); + } + val=energyCalibrationFunctions::median(x, nch); + h1->SetBinContent(ib+1,val); + } + } + delete [] x; + + return h1; + +} + + + + + + + + + + + + + + +void energyCalibration::setStartParameters(Double_t *par){ + bg_offset=par[0]; + bg_slope=par[1]; + flex=par[2]; + noise=par[3]; + ampl=par[4]; + cs_slope=par[5]; +} + +void energyCalibration::setStartParametersKb(Double_t *par){ + bg_offset=par[0]; + bg_slope=par[1]; + flex=par[2]; + noise=par[3]; + ampl=par[4]; + cs_slope=par[5]; + kb_mean=par[6]; + kb_frac=par[7]; + //fit_min = 400; // used for soleil flat field + //fit_max = 800; +} + +void energyCalibration::setStartParametersKaDoublet(Double_t *par){ + bg_offset=par[0]; + bg_slope=par[1]; + flex=par[2]; + noise=par[3]; + ampl=par[4]; + cs_slope=par[5]; + mean2=par[6]; + ampl2=par[7]; + //fit_min = 400; // used for soleil flat field + //fit_max = 800; +} + + +void energyCalibration::getStartParameters(Double_t *par){ + par[0]=bg_offset; + par[1]=bg_slope; + par[2]=flex; + par[3]=noise; + par[4]=ampl; + par[5]=cs_slope; +} + +#endif +int energyCalibration::setChargeSharing(int p) { + if (p>=0) { + cs_flag=p; +#ifdef MYROOT + if (p) { + fscurve->ReleaseParameter(5); + fspectrum->ReleaseParameter(1); + fspectrumkb->ReleaseParameter(1); + fspectrumkadoublet->ReleaseParameter(1); + } else { + fscurve->FixParameter(5,0); + fspectrum->FixParameter(1,0); + fspectrumkb->FixParameter(1,0); + fspectrumkadoublet->FixParameter(1,0); + } +#endif + } + + return cs_flag; +} + + +#ifdef MYROOT +void energyCalibration::initFitFunction(TF1 *fun, TH1 *h1) { + + Double_t min=fit_min, max=fit_max; + + Double_t mypar[6]; + + if (max==-1) + max=h1->GetXaxis()->GetXmax(); + + if (min==-1) + min=h1->GetXaxis()->GetXmin(); + + + if (bg_offset==-1) + mypar[0]=0; + else + mypar[0]=bg_offset; + + + if (bg_slope==-1) + mypar[1]=0; + else + mypar[1]=bg_slope; + + + if (flex==-1) + mypar[2]=(min+max)/2.; + else + mypar[2]=flex; + + + if (noise==-1) + mypar[3]=0.1; + else + mypar[3]=noise; + + if (ampl==-1) + mypar[4]=h1->GetBinContent(h1->GetXaxis()->FindBin(0.5*(max+min))); + else + mypar[4]=ampl; + + if (cs_slope==-1) + mypar[5]=0; + else + mypar[5]=cs_slope; + + fun->SetParameters(mypar); + + fun->SetRange(min,max); + +} + +void energyCalibration::initFitFunctionKb(TF1 *fun, TH1 *h1) { + + Double_t min=fit_min, max=fit_max; + + Double_t mypar[8]; + + if (max==-1) + max=h1->GetXaxis()->GetXmax(); + + if (min==-1) + min=h1->GetXaxis()->GetXmin(); + + + if (bg_offset==-1) + mypar[0]=0; + else + mypar[0]=bg_offset; + + + if (bg_slope==-1) + mypar[1]=0; + else + mypar[1]=bg_slope; + + + if (flex==-1) + mypar[2]=(min+max)/2.; + else + mypar[2]=flex; + + + if (noise==-1) + mypar[3]=0.1; + else + mypar[3]=noise; + + if (ampl==-1) + mypar[4]=h1->GetBinContent(h1->GetXaxis()->FindBin(0.5*(max+min))); + else + mypar[4]=ampl; + + if (cs_slope==-1) + mypar[5]=0; + else + mypar[5]=cs_slope; + + if (kb_mean==-1) + mypar[6]=0; + else + mypar[6]=kb_mean; + + if (kb_frac==-1) + mypar[7]=0; + else + mypar[7]=kb_frac; + + fun->SetParameters(mypar); + + fun->SetRange(min,max); + +} + +void energyCalibration::initFitFunctionKaDoublet(TF1 *fun, TH1 *h1) { + + Double_t min=fit_min, max=fit_max; + + Double_t mypar[8]; + + if (max==-1) + max=h1->GetXaxis()->GetXmax(); + + if (min==-1) + min=h1->GetXaxis()->GetXmin(); + + + if (bg_offset==-1) + mypar[0]=0; + else + mypar[0]=bg_offset; + + + if (bg_slope==-1) + mypar[1]=0; + else + mypar[1]=bg_slope; + + + if (flex==-1) + mypar[2]=(min+max)/2.; + else + mypar[2]=flex; + + + if (noise==-1) + mypar[3]=0.1; + else + mypar[3]=noise; + + if (ampl==-1) + mypar[4]=h1->GetBinContent(h1->GetXaxis()->FindBin(0.5*(max+min))); + else + mypar[4]=ampl; + + if (cs_slope==-1) + mypar[5]=0; + else + mypar[5]=cs_slope; + + if (mean2==-1) + mypar[6]=0; + else + mypar[6]=mean2; + + if (ampl2==-1) + mypar[7]=0; + else + mypar[7]=ampl2; + + fun->SetParameters(mypar); + + fun->SetRange(min,max); + +} + +TF1* energyCalibration::fitFunction(TF1 *fun, TH1 *h1, Double_t *mypar, Double_t *emypar) { + + + TF1* fitfun; + + char fname[100]; + + strcpy(fname, fun->GetName()); + + if (plot_flag) { + h1->Fit(fname,"R0Q"); + } else + h1->Fit(fname,"R0Q"); + + fitfun= h1->GetFunction(fname); + fitfun->GetParameters(mypar); + for (int ip=0; ip<6; ip++) { + emypar[ip]=fitfun->GetParError(ip); + } + return fitfun; +} + +TF1* energyCalibration::fitFunctionKb(TF1 *fun, TH1 *h1, Double_t *mypar, Double_t *emypar) { + + + TF1* fitfun; + + char fname[100]; + + strcpy(fname, fun->GetName()); + + if (plot_flag) { + h1->Fit(fname,"R0Q"); + } else + h1->Fit(fname,"R0Q"); + + fitfun= h1->GetFunction(fname); + fitfun->GetParameters(mypar); + for (int ip=0; ip<8; ip++) { + emypar[ip]=fitfun->GetParError(ip); + } + return fitfun; +} + +TF1* energyCalibration::fitFunctionKaDoublet(TF1 *fun, TH1 *h1, Double_t *mypar, Double_t *emypar) { + + + TF1* fitfun; + + char fname[100]; + + strcpy(fname, fun->GetName()); + + if (plot_flag) { + h1->Fit(fname,"R0Q"); + } else + h1->Fit(fname,"R0Q"); + + + fitfun= h1->GetFunction(fname); + fitfun->GetParameters(mypar); + for (int ip=0; ip<8; ip++) { + emypar[ip]=fitfun->GetParError(ip); + } + return fitfun; +} + +TF1* energyCalibration::fitSCurve(TH1 *h1, Double_t *mypar, Double_t *emypar) { + initFitFunction(fscurve,h1); + return fitFunction(fscurve, h1, mypar, emypar); +} + + + + + +TF1* energyCalibration::fitSpectrum(TH1 *h1, Double_t *mypar, Double_t *emypar) { + initFitFunction(fspectrum,h1); + return fitFunction(fspectrum, h1, mypar, emypar); +} + +TF1* energyCalibration::fitSpectrumKb(TH1 *h1, Double_t *mypar, Double_t *emypar) { + initFitFunctionKb(fspectrumkb,h1); + return fitFunctionKb(fspectrumkb, h1, mypar, emypar); +} + +TF1* energyCalibration::fitSpectrumKaDoublet(TH1 *h1, Double_t *mypar, Double_t *emypar) { + initFitFunctionKaDoublet(fspectrumkadoublet,h1); + return fitFunctionKaDoublet(fspectrumkadoublet, h1, mypar, emypar); +} + + +TGraphErrors* energyCalibration::linearCalibration(int nscan, Double_t *en, Double_t *een, Double_t *fl, Double_t *efl, Double_t &gain, Double_t &off, Double_t &egain, Double_t &eoff) { + + TGraphErrors *gr; + + Double_t mypar[2]; + + gr = new TGraphErrors(nscan,en,fl,een,efl); + + if (plot_flag) { + gr->Fit("pol1"); + gr->SetMarkerStyle(20); + } else + gr->Fit("pol1","0Q"); + + TF1 *fitfun= gr->GetFunction("pol1"); + fitfun->GetParameters(mypar); + + egain=fitfun->GetParError(1); + eoff=fitfun->GetParError(0); + + gain=funcs->setScanSign()*mypar[1]; + + off=mypar[0]; + + return gr; +} + + +TGraphErrors* energyCalibration::calibrate(int nscan, Double_t *en, Double_t *een, TH1F **h1, Double_t &gain, Double_t &off, Double_t &egain, Double_t &eoff, int integral) { + + TH1F *h; + + Double_t mypar[6], emypar[6]; + Double_t fl[nscan], efl[nscan]; + + + for (int ien=0; ien + + +//#include "moench03T1ZmqData.h" +#include "jungfrauHighZSingleChipData.h" + + +#include "multiThreadedAnalogDetector.h" +#include "singlePhotonDetector.h" + +#include +#include +#include +#include + +#include +using namespace std; + + +int main(int argc, char *argv[]) { + + + if (argc<6) { + cout << "Usage is " << argv[0] << "indir outdir fname runmin runmax " << endl; + return 1; + } + int p=10000; + int fifosize=1000; + int nthreads=1; + int nsubpix=25; + int etabins=nsubpix*10; + double etamin=-1, etamax=2; + int csize=3; + int nx=400, ny=400; + int save=1; + int nsigma=5; + int nped=1000; + int ndark=100; + int ok; + int iprog=0; + + + + + jungfrauHighZSingleChipData *decoder=new jungfrauHighZSingleChipData(); + + decoder->getDetectorSize(nx,ny); + cout << "nx " << nx << " ny " << ny << endl; + + //moench03T1ZmqData *decoder=new moench03T1ZmqData(); + singlePhotonDetector *filter=new singlePhotonDetector(decoder,csize, nsigma, 1, 0, nped, 200); + // char tit[10000]; + cout << "filter " << endl; + + + + + int* image; + filter->newDataSet(); + + + int ff, np; + int dsize=decoder->getDataSize(); + cout << " data size is " << dsize; + + + char data[dsize]; + + ifstream filebin; + char *indir=argv[1]; + char *outdir=argv[2]; + char *fformat=argv[3]; + int runmin=atoi(argv[4]); + int runmax=atoi(argv[5]); + + char fname[10000]; + char outfname[10000]; + char imgfname[10000]; + char pedfname[10000]; + char fn[10000]; + + std::time_t end_time; + + FILE *of=NULL; + cout << "input directory is " << indir << endl; + cout << "output directory is " << outdir << endl; + cout << "fileformat is " << fformat << endl; + + + std::time(&end_time); + cout << std::ctime(&end_time) << endl; + + char* buff; + multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize); + + + mt->setDetectorMode(ePhotonCounting); + mt->setFrameMode(eFrame); + mt->StartThreads(); + mt->popFree(buff); + + + cout << "mt " << endl; + + int ifr=0; + + + for (int irun=runmin; irunsetFilePointer(of); + // cout << "file pointer set " << endl; + } else { + cout << "Could not open "<< outfname << " for writing " << endl; + mt->setFilePointer(NULL); + return 1; + } + // //while read frame + ff=-1; + while (decoder->readNextFrame(filebin, ff, np,buff)) { + + mt->pushData(buff); + mt->nextThread(); + mt->popFree(buff); + ifr++; + if (ifr%10000==0) cout << ifr << " " << ff << endl; + ff=-1; + } + cout << "--" << endl; + filebin.close(); + while (mt->isBusy()) {;}//wait until all data are processed from the queues + if (of) + fclose(of); + + mt->writeImage(imgfname); + mt->clearImage(); + + std::time(&end_time); + cout << std::ctime(&end_time) << endl; + + } else + cout << "Could not open "<< fname << " for reading " << endl; + + + } + + + return 0; +} + diff --git a/slsDetectorCalibration/jungfrauExecutables/jungfrauClusterFinder.cpp~ b/slsDetectorCalibration/jungfrauExecutables/jungfrauClusterFinder.cpp~ new file mode 100644 index 000000000..36a6675af --- /dev/null +++ b/slsDetectorCalibration/jungfrauExecutables/jungfrauClusterFinder.cpp~ @@ -0,0 +1,237 @@ +//#include "sls/ansi.h" +#include + + +//#include "moench03T1ZmqData.h" +#ifdef NEWRECEIVER +#ifndef RECT +#include "moench03T1ReceiverDataNew.h" +#endif + +#ifdef RECT +#include "moench03T1ReceiverDataNewRect.h" +#endif + +#endif + + +#ifdef CSAXS_FP +#include "moench03T1ReceiverData.h" +#endif +#ifdef OLDDATA +#include "moench03Ctb10GbT1Data.h" +#endif + +#ifdef REORDERED +#include "moench03T1ReorderedData.h" +#endif + +// #include "interpolatingDetector.h" +//#include "etaInterpolationPosXY.h" +// #include "linearInterpolation.h" +// #include "noInterpolation.h" +#include "multiThreadedAnalogDetector.h" +#include "singlePhotonDetector.h" +//#include "interpolatingDetector.h" + +#include +#include +#include +#include + +#include +using namespace std; + + +int main(int argc, char *argv[]) { + + + if (argc<6) { + cout << "Usage is " << argv[0] << "indir outdir fname runmin runmax " << endl; + return 1; + } + int p=10000; + int fifosize=1000; + int nthreads=1; + int nsubpix=25; + int etabins=nsubpix*10; + double etamin=-1, etamax=2; + int csize=3; + int nx=400, ny=400; + int save=1; + int nsigma=5; + int nped=1000; + int ndark=100; + int ok; + int iprog=0; + + + + + +#ifdef NEWRECEIVER +#ifdef RECT + cout << "Should be rectangular!" <getDetectorSize(nx,ny); + cout << "nx " << nx << " ny " << ny << endl; + + //moench03T1ZmqData *decoder=new moench03T1ZmqData(); + singlePhotonDetector *filter=new singlePhotonDetector(decoder,csize, nsigma, 1, 0, nped, 200); + // char tit[10000]; + cout << "filter " << endl; + + + + // filter->readPedestals("/scratch/ped_100.tiff"); + // interp->readFlatField("/scratch/eta_100.tiff",etamin,etamax); + // cout << "filter "<< endl; + + + int size = 327680;////atoi(argv[3]); + + int* image; + //int* image =new int[327680/sizeof(int)]; + filter->newDataSet(); + + + int ff, np; + int dsize=decoder->getDataSize(); + cout << " data size is " << dsize; + + + char data[dsize]; + + ifstream filebin; + char *indir=argv[1]; + char *outdir=argv[2]; + char *fformat=argv[3]; + int runmin=atoi(argv[4]); + int runmax=atoi(argv[5]); + + char fname[10000]; + char outfname[10000]; + char imgfname[10000]; + char pedfname[10000]; + // strcpy(pedfname,argv[6]); + char fn[10000]; + + std::time_t end_time; + + FILE *of=NULL; + cout << "input directory is " << indir << endl; + cout << "output directory is " << outdir << endl; + cout << "fileformat is " << fformat << endl; + + + std::time(&end_time); + cout << std::ctime(&end_time) << endl; + + + + + + + + + char* buff; + multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize); + + + mt->setDetectorMode(ePhotonCounting); + mt->setFrameMode(eFrame); + mt->StartThreads(); + mt->popFree(buff); + + + cout << "mt " << endl; + + int ifr=0; + + + for (int irun=runmin; irunsetFilePointer(of); + // cout << "file pointer set " << endl; + } else { + cout << "Could not open "<< outfname << " for writing " << endl; + mt->setFilePointer(NULL); + return 1; + } + // //while read frame + ff=-1; + while (decoder->readNextFrame(filebin, ff, np,buff)) { + // cout << "*"<getChannel(buff, ix, iy)<3000 || decoder->getChannel(buff, ix, iy)>8000) { + // cout << ifr << " " << ff << " " << ix << " " << iy << " " << decoder->getChannel(buff, ix, iy) << endl ; + // } + // } + + mt->pushData(buff); + // // //pop + mt->nextThread(); + // // // cout << " " << (void*)buff; + mt->popFree(buff); + ifr++; + if (ifr%10000==0) cout << ifr << " " << ff << endl; + ff=-1; + } + cout << "--" << endl; + filebin.close(); + // //close file + // //join threads + while (mt->isBusy()) {;}//wait until all data are processed from the queues + if (of) + fclose(of); + + mt->writeImage(imgfname); + mt->clearImage(); + + std::time(&end_time); + cout << std::ctime(&end_time) << endl; + + } else + cout << "Could not open "<< fname << " for reading " << endl; + + + } + + + return 0; +} + diff --git a/slsDetectorCalibration/jungfrauExecutables/jungfrauInterpolation.cpp b/slsDetectorCalibration/jungfrauExecutables/jungfrauInterpolation.cpp new file mode 100644 index 000000000..eeeed90c6 --- /dev/null +++ b/slsDetectorCalibration/jungfrauExecutables/jungfrauInterpolation.cpp @@ -0,0 +1,250 @@ + +#include "sls/ansi.h" +#include + +//#include "moench03T1ZmqData.h" +//#define DOUBLE_SPH +//#define MANYFILES + +#ifdef DOUBLE_SPH +#include "single_photon_hit_double.h" +#endif + +#ifndef DOUBLE_SPH +#include "single_photon_hit.h" +#endif + +//#include "etaInterpolationPosXY.h" +#include "noInterpolation.h" +#include "etaInterpolationPosXY.h" +//#include "etaInterpolationCleverAdaptiveBins.h" +//#include "etaInterpolationRandomBins.h" +using namespace std; +#define NC 400 +#define NR 400 +#define MAX_ITERATIONS (nSubPixels*100) + +#define XTALK + +int main(int argc, char *argv[]) { + +#ifndef FF + if (argc<9) { + cout << "Wrong usage! Should be: "<< argv[0] << " infile etafile outfile runmin runmax ns cmin cmax" << endl; + return 1; + } +#endif + +#ifdef FF + if (argc<7) { + cout << "Wrong usage! Should be: "<< argv[0] << " infile etafile runmin runmax cmin cmax" << endl; + return 1; + } +#endif + int iarg=4; + char infname[10000]; + char fname[10000]; + char outfname[10000]; +#ifndef FF + iarg=4; +#endif + +#ifdef FF + iarg=3; +#endif + int runmin=atoi(argv[iarg++]); + int runmax=atoi(argv[iarg++]); + cout << "Run min: " << runmin << endl; + cout << "Run max: " << runmax << endl; + + int nsubpix=4; +#ifndef FF + nsubpix=atoi(argv[iarg++]); + cout << "Subpix: " << nsubpix << endl; +#endif + float cmin=atof(argv[iarg++]); + float cmax=atof(argv[iarg++]); + cout << "Energy min: " << cmin << endl; + cout << "Energy max: " << cmax << endl; + //int etabins=500; + int etabins=1000;//nsubpix*2*100; + double etamin=-1, etamax=2; + //double etamin=-0.1, etamax=1.1; + double eta3min=-2, eta3max=2; + int quad; + double sum, totquad; + double sDum[2][2]; + double etax, etay, int_x, int_y; + double eta3x, eta3y, int3_x, int3_y, noint_x, noint_y; + int ok; + int f0=-1; + int ix, iy, isx, isy; + int nframes=0, lastframe=-1; + double d_x, d_y, res=5, xx, yy; + int nph=0, badph=0, totph=0; + FILE *f=NULL; + +#ifdef DOUBLE_SPH + single_photon_hit_double cl(3,3); +#endif + +#ifndef DOUBLE_SPH + single_photon_hit cl(3,3); +#endif + + int nSubPixels=nsubpix; +#ifndef NOINTERPOLATION + eta2InterpolationPosXY *interp=new eta2InterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax); + //eta2InterpolationCleverAdaptiveBins *interp=new eta2InterpolationCleverAdaptiveBins(NC, NR, nsubpix, etabins, etamin, etamax); +#endif +#ifdef NOINTERPOLATION + noInterpolation *interp=new noInterpolation(NC, NR, nsubpix); +#endif + + + +#ifndef FF +#ifndef NOINTERPOLATION + cout << "read ff " << argv[2] << endl; + sprintf(fname,"%s",argv[2]); + interp->readFlatField(fname); + interp->prepareInterpolation(ok);//, MAX_ITERATIONS); +#endif + // return 0; +#endif +#ifdef FF + cout << "Will write eta file " << argv[2] << endl; +#endif + + int *img; + float *totimg=new float[NC*NR*nsubpix*nsubpix]; + for (ix=0; ixcalcQuad(cl.get_cluster(), sum, totquad, sDum); + quad=interp->calcEta(cl.get_cluster(), etax, etay, sum, totquad, sDum); + if (sum>cmin && totquad/sum>0.8 && totquad/sum<1.2 && sum200 && sum<580) { + // interp->getInterpolatedPosition(cl.x,cl.y, totquad,quad,cl.get_cluster(),int_x, int_y); +// #ifdef SOLEIL +// if (cl.x>210 && cl.x<240 && cl.y>210 && cl.y<240) { +// #endif +#ifndef FF + // interp->getInterpolatedPosition(cl.x,cl.y, cl.get_cluster(),int_x, int_y); + interp->getInterpolatedPosition(cl.x,cl.y, etax, etay, quad,int_x, int_y); + // cout <<"**************"<< endl; + // cout << cl.x << " " << cl.y << " " << sum << endl; + // cl.print(); + // cout << int_x << " " << int_y << endl; + // cout <<"**************"<< endl; + // if (etax!=0 && etay!=0 && etax!=1 && etay!=1) + interp->addToImage(int_x, int_y); + if (int_x<0 || int_y<0 || int_x>400 || int_y>400) { + cout <<"**************"<< endl; + cout << cl.x << " " << cl.y << " " << sum << endl; + cl.print(); + cout << int_x << " " << int_y << endl; + cout <<"**************"<< endl; + } +#endif +#ifdef FF + // interp->addToFlatField(cl.get_cluster(), etax, etay); +// #ifdef UCL +// if (cl.x>50) +// #endif +// if (etax!=0 && etay!=0 && etax!=1 && etay!=1) + interp->addToFlatField(etax, etay); + // if (etax==0 || etay==0) cout << cl.x << " " << cl.y << endl; + +#endif +// #ifdef SOLEIL +// } +// #endif + + if (nph%1000000==0) cout << nph << endl; + if (nph%10000000==0) { +#ifndef FF + interp->writeInterpolatedImage(outfname); +#endif +#ifdef FF + interp->writeFlatField(outfname); +#endif + + } + } + + } + + + fclose(f); +#ifdef FF + interp->writeFlatField(outfname); +#endif + +#ifndef FF + interp->writeInterpolatedImage(outfname); + + img=interp->getInterpolatedImage(); + for (ix=0; ixwriteFlatField(outfname); +#endif + + cout << "Filled " << nph << " (/"<< totph <<") " << endl; + return 0; +} + diff --git a/slsDetectorCalibration/jungfrauExecutables/jungfrauPhotonCounter.cpp b/slsDetectorCalibration/jungfrauExecutables/jungfrauPhotonCounter.cpp new file mode 100644 index 000000000..7b765bd5e --- /dev/null +++ b/slsDetectorCalibration/jungfrauExecutables/jungfrauPhotonCounter.cpp @@ -0,0 +1,453 @@ +//#include "sls/ansi.h" +#include +#define CORR + +#define C_GHOST 0.0004 + +#define CM_ROWS 50 + +//#define VERSION_V1 + +//#include "moench03T1ZmqData.h" +#ifdef NEWRECEIVER +#ifndef RECT +#include "moench03T1ReceiverDataNew.h" +#endif + +#ifdef RECT +#include "moench03T1ReceiverDataNewRect.h" +#endif + +#endif + + +#ifdef CSAXS_FP +#include "moench03T1ReceiverData.h" +#endif +#ifdef OLDDATA +#include "moench03Ctb10GbT1Data.h" +#endif + +// #include "interpolatingDetector.h" +//#include "etaInterpolationPosXY.h" +// #include "linearInterpolation.h" +// #include "noInterpolation.h" +#include "multiThreadedCountingDetector.h" +//#include "multiThreadedAnalogDetector.h" +#include "singlePhotonDetector.h" +#include "moench03GhostSummation.h" +#include "moench03CommonMode.h" +//#include "interpolatingDetector.h" + +#include +#include +#include +#include + +#include +using namespace std; + + +int main(int argc, char *argv[]) { + + + if (argc<4) { + cout << "Usage is " << argv[0] << "indir outdir fname [runmin] [runmax] [pedfile] [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 p=10000; + int fifosize=1000; + int nthreads=10; + int nsubpix=25; + int etabins=nsubpix*10; + double etamin=-1, etamax=2; + int csize=3; + int save=1; + int nsigma=5; + int nped=10000; + int ndark=100; + int ok; + int iprog=0; + + int cf=0; + +#ifdef NEWRECEIVER +#ifdef RECT + cout << "Should be rectangular!" <getDetectorSize(nx,ny); + + int ncol_cm=CM_ROWS; + double xt_ghost=C_GHOST; + moench03CommonMode *cm=NULL; + moench03GhostSummation *gs; + double *gainmap=NULL; + float *gm; + + + + int size = 327680;////atoi(argv[3]); + + int* image; + //int* image =new int[327680/sizeof(int)]; + + int ff, np; + //cout << " data size is " << dsize; + + + + ifstream filebin; + char *indir=argv[1]; + char *outdir=argv[2]; + char *fformat=argv[3]; + int runmin=0; + + // cout << "argc is " << argc << endl; + if (argc>=5) { + runmin=atoi(argv[4]); + } + + int runmax=runmin; + + if (argc>=6) { + runmax=atoi(argv[5]); + } + + char *pedfile=NULL; + if (argc>=7) { + pedfile=argv[6]; + } + double thr=0; + double thr1=1; + + if (argc>=8) { + thr=atof(argv[7]); + } + + + int nframes=0; + + if (argc>=9) { + nframes=atoi(argv[8]); + } + + int xmin=0, xmax=nx, ymin=0, ymax=ny; + if (argc>=13) { + xmin=atoi(argv[9]); + xmax=atoi(argv[10]); + ymin=atoi(argv[11]); + ymax=atoi(argv[12]); + } + + + char *gainfname=NULL; + if (argc>13) { + gainfname=argv[13]; + cout << "Gain map file name is: " << gainfname << endl; + } + + + + + char ffname[10000]; + char fname[10000]; + char imgfname[10000]; + char cfname[10000]; + char fn[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 nnx, nny; + double *gmap; + + // if (gainfname) { + // gm=ReadFromTiff(gainfname, nny, nnx); + // if (gm && nnx==nx && nny==ny) { + // gmap=new double[nx*ny]; + // for (int i=0; ireadGainMap(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(); + + + char data[dsize]; + + + + + + //#ifndef ANALOG + if (thr>0) { + cout << "threshold is " << thr << endl; + //#ifndef ANALOG + filter->setThreshold(thr); + //#endif + cf=0; + + } else + cf=1; + //#endif + + + 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); +#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; + + double ped[nx*ny], *ped1; + + + if (pedfile) { + + cout << "PEDESTAL " << endl; + sprintf(imgfname,"%s/pedestals.tiff",outdir); + + if (string(pedfile).find(".tif")==std::string::npos){ + sprintf(fname,"%s.raw",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) { + mt->pushData(buff); + mt->nextThread(); + mt->popFree(buff); + ifr++; + if (ifr%100==0) + cout << ifr << " " << ff << " " << np << endl; + } else + cout << ifr << " " << ff << " " << np << endl; + ff=-1; + } + filebin.close(); + while (mt->isBusy()) {;} + + } else + cout << "Could not open pedestal file "<< fname << " for reading " << endl; + } else { + float *pp=ReadFromTiff(pedfile, nny, nnx); + if (pp && nnx==nx && nny==ny) { + for (int i=0; isetPedestal(ped); + // ped1=mt->getPedestal(); + + // for (int i=0; iwritePedestal(imgfname); + 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.raw",indir,fformat); + sprintf(fname,ffname,irun); + sprintf(ffname,"%s/%s.tiff",outdir,fformat); + sprintf(imgfname,ffname,irun); + sprintf(ffname,"%s/%s.clust",outdir,fformat); + sprintf(cfname,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) { + // cout << "*"<pushData(buff); + // // //pop + mt->nextThread(); + // // // cout << " " << (void*)buff; + mt->popFree(buff); + + + + + ifr++; + if (ifr%100==0) cout << ifr << " " << ff << endl; + if (nframes>0) { + if (ifr%nframes==0) { + //The name has an additional "_fXXXXX" at the end, where "XXXXX" is the initial frame number of the image (0,1000,2000...) + + sprintf(ffname,"%s/%s_f%05d.tiff",outdir,fformat,ifile); + sprintf(imgfname,ffname,irun); + //cout << "Writing tiff to " << imgfname << " " << thr1 << endl; + mt->writeImage(imgfname, thr1); + mt->clearImage(); + ifile++; + } + } + } else + cout << ifr << " " << ff << " " << np << endl; + ff=-1; + } + cout << "--" << endl; + filebin.close(); + // //close file + // //join threads + while (mt->isBusy()) {;} + if (nframes>=0) { + if (nframes>0) { + sprintf(ffname,"%s/%s_f%05d.tiff",outdir,fformat,ifile); + sprintf(imgfname,ffname,irun); + } else { + sprintf(ffname,"%s/%s.tiff",outdir,fformat); + sprintf(imgfname,ffname,irun); + } + cout << "Writing tiff to " << imgfname << " " << thr1 <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 <writeImage(imgfname, thr1); + } + + + + return 0; +} + diff --git a/slsDetectorCalibration/jungfrauExecutables/jungfrauZmqProcess.cpp b/slsDetectorCalibration/jungfrauExecutables/jungfrauZmqProcess.cpp new file mode 100644 index 000000000..47e424306 --- /dev/null +++ b/slsDetectorCalibration/jungfrauExecutables/jungfrauZmqProcess.cpp @@ -0,0 +1,1026 @@ +//#define WRITE_QUAD +#define DEVELOPER +#undef CORR +#define C_GHOST 0.0004 + +#define CM_ROWS 20 + +#include "sls/sls_detector_defs.h" +#include "sls/ZmqSocket.h" +#ifndef RECT +#ifndef MOENCH04 +#include "moench03T1ZmqDataNew.h" +#endif +#ifdef MOENCH04 +#include "moench04CtbZmq10GbData.h" +#endif + +#endif +#ifdef RECT +#include "moench03T1ZmqDataNewRect.h" +#endif +#include "moench03GhostSummation.h" +#include "moench03CommonMode.h" +#include +#include +#include +#include +#include +#include "tiffIO.h" + +#include //json header in zmq stream + +#include + +//#include "analogDetector.h" +//#include "multiThreadedAnalogDetector.h" +//#include "singlePhotonDetector.h" +//#include "interpolatingDetector.h" +//#include "multiThreadedCountingDetector.h" +#include "multiThreadedInterpolatingDetector.h" +#include "etaInterpolationPosXY.h" +#include "sls/ansi.h" +#include + +#include +#include // time_t +#include + +using namespace std; +using namespace std::chrono; + +//#define SLS_DETECTOR_JSON_HEADER_VERSION 0x2 + + // myDet->setNetworkParameter(ADDITIONAL_JSON_HEADER, " \"what\":\"nothing\" "); + +int main(int argc, char *argv[]) { +/** + * trial.o [socket ip] [starting port number] [send_socket ip] [send port number] + * + */ + FILE *of=NULL; + int fifosize=5000; + int etabins=1000, etabinsy=1000;//nsubpix*2*100; + double etamin=-1, etamax=2; + int nSubPixelsX=2; + // int emin, emax; + int nSubPixelsY=2; + // help + if (argc < 3 ) { + cprintf(RED, "Help: ./trial [receive socket ip] [receive starting port number] [send_socket ip] [send starting port number] [nthreads] [nsubpix] [gainmap] [etafile]\n"); + return EXIT_FAILURE; + } + + // receive parameters + bool send = false; + char* socketip=argv[1]; + uint32_t portnum = atoi(argv[2]); + // send parameters if any + char* socketip2 = 0; + uint32_t portnum2 = 0; + + zmqHeader zHeader, outHeader; + zHeader.jsonversion = SLS_DETECTOR_JSON_HEADER_VERSION; + outHeader.jsonversion = SLS_DETECTOR_JSON_HEADER_VERSION; + + uint32_t nSigma=5; + + int ok; + + high_resolution_clock::time_point t1; + high_resolution_clock::time_point t2 ; + std::chrono::steady_clock::time_point begin,end,finished; + //time_t begin,end,finished; + int rms=0; + + + + if (argc > 4) { + socketip2 = argv[3]; + portnum2 = atoi(argv[4]); + if (portnum2>0) + send = true; + } + cout << "\nrx socket ip : " << socketip << + "\nrx port num : " << portnum ; + if (send) { + cout << "\ntx socket ip : " << socketip2 << + "\ntx port num : " << portnum2; + } + int nthreads=5; + if (argc>5) + nthreads=atoi(argv[5]); + + cout << "Number of threads is: " << nthreads << endl; + if (argc>6) { + nSubPixelsX=atoi(argv[6]); + nSubPixelsY=nSubPixelsX; +#ifdef RECT + nSubPixelsX=2; +#endif + } + cout << "Number of subpixels is: " << nSubPixelsX << " " << nSubPixelsY << endl; + + char *gainfname=NULL; + if (argc>7) { + gainfname=argv[7]; + cout << "Gain map file name is: " << gainfname << endl; + } + + char *etafname=NULL; + if (argc>8) { + etafname=argv[8]; + cout << "Eta file name is: " << etafname << endl; + } + + //slsDetectorData *det=new moench03T1ZmqDataNew(); +#ifndef MOENCH04 + moench03T1ZmqDataNew *det=new moench03T1ZmqDataNew(); +#endif +#ifdef MOENCH04 + moench04CtbZmq10GbData *det=new moench04CtbZmq10GbData(); +#endif + cout << endl << " det" <getDetectorSize(npx, npy); + + int send_something=0; + + + int maxSize = npx*npy*2;//32*2*8192;//5000;//atoi(argv[3]); + int size= maxSize;//32*2*5000; + //int multisize=size; + //int dataSize=size; + + char dummybuff[size]; + + + moench03CommonMode *cm=NULL; + moench03GhostSummation *gs=NULL; +#ifdef CORR + + //int ncol_cm=CM_ROWS; + //double xt_ghost=C_GHOST; + + cm=new moench03CommonMode(CM_ROWS); + gs=new moench03GhostSummation(det, C_GHOST); +#endif + double *gainmap=NULL; + float *gm; + double *gmap=NULL; + + uint32_t nnnx, nnny; + if (gainfname) { + gm=ReadFromTiff(gainfname, nnny, nnnx); + if (gm && nnnx==(uint)npx && nnny==(uint)npy) { + gmap=new double[npx*npy]; + for (int i=0; i *filter=new analogDetector(det,1,NULL,1000); +#ifndef INTERP + singlePhotonDetector *filter=new singlePhotonDetector(det,3, nSigma, 1, cm, 1000, 100, -1, -1, gainmap, gs); + + multiThreadedCountingDetector *mt=new multiThreadedCountingDetector(filter,nthreads,fifosize); + + // multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize); +#endif +#ifdef INTERP + eta2InterpolationPosXY *interp=new eta2InterpolationPosXY(npx, npy, nSubPixelsX,nSubPixelsY, etabins, etabinsy, etamin, etamax); + + if (etafname) interp->readFlatField(etafname); + + interpolatingDetector *filter=new interpolatingDetector(det,interp, nSigma, 1, cm, 1000, 10, -1, -1, gainmap, gs); + multiThreadedInterpolatingDetector *mt=new multiThreadedInterpolatingDetector(filter,nthreads,fifosize); +#endif + + + + char* buff; + mt->setFrameMode(eFrame); + mt->StartThreads(); + mt->popFree(buff); + + ZmqSocket* zmqsocket=NULL; + +#ifdef NEWZMQ + // receive socket + try{ +#endif + + zmqsocket = new ZmqSocket(socketip,portnum); + + +#ifdef NEWZMQ + } catch (...) { + cprintf(RED, "Error: Could not create Zmq socket on port %d with ip %s\n", portnum, socketip); + delete zmqsocket; + return EXIT_FAILURE; + } +#endif + +#ifndef NEWZMQ + if (zmqsocket->IsError()) { + cprintf(RED, "Error: Could not create Zmq socket on port %d with ip %s\n", portnum, socketip); + delete zmqsocket; + return EXIT_FAILURE; + } +#endif + if (zmqsocket->Connect()) { + cprintf(RED, "Error: Could not connect to socket %s\n", + (zmqsocket->GetZmqServerAddress()).c_str()); + delete zmqsocket; + return EXIT_FAILURE; + } else + printf("Zmq Client at %s\n", zmqsocket->GetZmqServerAddress().c_str()); + + // send socket + ZmqSocket* zmqsocket2 = 0; + // cout << "zmq2 " << endl; + if (send) { +#ifdef NEWZMQ + // receive socket + try{ +#endif + zmqsocket2 = new ZmqSocket(portnum2, socketip2); + + + +#ifdef NEWZMQ + } catch (...) { + cprintf(RED, "Error: Could not create Zmq socket server on port %d and ip %s\n", portnum2, socketip2); + // delete zmqsocket2; + // zmqsocket2=NULL; + // delete zmqsocket; + // return EXIT_FAILURE; + send = false; + } +#endif + +#ifndef NEWZMQ + if (zmqsocket2->IsError()) { + cprintf(RED, "AAA Error: Could not create Zmq socket server on port %d and ip %s\n", portnum2, socketip2); + // delete zmqsocket2; + //delete zmqsocket; + // return EXIT_FAILURE; + send = false; + } +#endif + if (zmqsocket2->Connect()) { + cprintf(RED, "BBB Error: Could not connect to socket %s\n", + zmqsocket2->GetZmqServerAddress().c_str()); + // delete zmqsocket2; + send = false; + // return EXIT_FAILURE; + } else + printf("Zmq Client at %s\n", zmqsocket2->GetZmqServerAddress().c_str()); + + } + + + // header variables + uint64_t acqIndex = -1; + uint64_t frameIndex = -1; +#ifdef MOENCH_BRANCH + uint32_t subFrameIndex = -1; + int* flippedData = 0; +#endif + + uint64_t subframes=0; + //uint64_t isubframe=0; + uint64_t insubframe=0; + double subnorm=1; + uint64_t f0=-1, nsubframes=0, nnsubframe=0; + + uint64_t fileindex = -1; + string filename = ""; + // char* image = new char[size]; + //int* image = new int[(size/sizeof(int))](); + //uint32_t flippedDataX = -1; + //int *nph; + int iframe=0; + char ofname[10000]; + + string fname; + // int length; + int *detimage=NULL; + int nnx, nny,nnsx, nnsy; + //uint32_t imageSize = 0, nPixelsX = 0, nPixelsY = 0, + //uint32_t dynamicRange = 0; + // infinite loop + uint32_t packetNumber = 0; + uint64_t bunchId = 0; + uint64_t timestamp = 0; + int16_t modId = 0; + uint32_t expLength=0; + uint16_t xCoord = 0; + uint16_t yCoord = 0; + //uint16_t zCoord = 0; + uint32_t debug = 0; + //uint32_t dr = 16; + //int16_t *dout;//=new int16_t [nnx*nny]; + uint32_t dr = 32; + int32_t *dout=NULL;//=new int32_t [nnx*nny]; + float *doutf=NULL;//=new int32_t [nnx*nny]; + uint16_t roundRNumber = 0; + uint8_t detType = 0; + uint8_t version = 0; + string additionalJsonHeader="" ; + + int32_t threshold=0; + + int32_t xmin=0, xmax=400, ymin=0, ymax=400; + + string frameMode_s, detectorMode_s, intMode_s; + + // int resetFlat=0; + //int resetPed=0; + // int nsubPixels=1; + //int isPedestal=0; + //int isFlat=0; + int newFrame=1; + detectorMode dMode=eAnalog; + frameMode fMode=eFrame; + double *ped; + + filter->getImageSize(nnx, nny,nnsx, nnsy); + + + std::map addJsonHeader; + + + + + while(1) { + + + // cout << "+++++++++++++++++++++++++++++++LOOP" << endl; + // get header, (if dummy, fail is on parse error or end of acquisition) + + + + // rapidjson::Document doc; + if (!zmqsocket->ReceiveHeader(0, zHeader, SLS_DETECTOR_JSON_HEADER_VERSION)) { + /* zmqsocket->CloseHeaderMessage();*/ + + // if (!zmqsocket->ReceiveHeader(0, acqIndex, frameIndex, subframeIndex, filename, fileindex)) { + cprintf(RED, "Got Dummy\n"); + // t1=high_resolution_clock::now(); + //time(&end); + //cout << "Measurement lasted " << difftime(end,begin) << endl; + + end = std::chrono::steady_clock::now(); + cout << "Measurement lasted " << (end-begin).count()*0.000001 << " ms" << endl; + + while (mt->isBusy()) {;}//wait until all data are processed from the queues + + if (of) { + mt->setFilePointer(NULL); + fclose(of); + of=NULL; + } + if (newFrame>0) { + cprintf(RED,"DIDn't receive any data!\n"); + if (send) { + + //zHeader.data = false; + outHeader.data=false; + // zmqsocket2->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION); + zmqsocket2->SendHeader(0,outHeader); + cprintf(RED, "Sent Dummy\n"); + } + } else { + send_something=0; + if (fMode==ePedestal) { + sprintf(ofname,"%s_%ld_ped.tiff",fname.c_str(),fileindex); + mt->writePedestal(ofname); + cout << "Writing pedestal to " << ofname << endl; + if (rms){ + sprintf(ofname,"%s_%ld_var.tiff",fname.c_str(),fileindex); + mt->writePedestalRMS(ofname); + + } + send_something=1; + } +#ifdef INTERP + else if (fMode==eFlat) { + mt->prepareInterpolation(ok); + sprintf(ofname,"%s_%ld_eta.tiff",fname.c_str(),fileindex); + mt->writeFlatField(ofname); + cout << "Writing eta to " << ofname << endl; + send_something=1; + } +#endif + else { + if (subframes>0 ) { + if (insubframe>0) { + sprintf(ofname,"%s_sf%ld_%ld.tiff",fname.c_str(),nnsubframe,fileindex); + // mt->writeImage(ofname); + doutf= new float[nnx*nny]; + if (subframes>0 && insubframe!=subframes && insubframe>0) + subnorm=((double)subframes)/((double)insubframe); + else + subnorm=1.; + for (int ix=0; ixwriteImage(ofname); + send_something=1; + } + + cout << "Writing image to " << ofname << endl; + } + // cout << nns*nnx*nny*nns*dr/8 << " " << length << endl; + + if (send) { + + if (fMode==ePedestal) { + cprintf(MAGENTA,"Get pedestal!\n"); + nnsx=1; + nnsy=1; + + nnx=npx; + nny=npy; + //dout= new int16_t[nnx*nny*nns*nns]; + dout= new int32_t[nnx*nny*nnsx*nnsy]; + // cout << "get pedestal " << endl; + ped=mt->getPedestal(); + // cout << "got pedestal " << endl; + for (int ix=0; ixgetFlatField(nb, emi, ema); + nnx=nb; + nny=nb; + dout= new int32_t[nb*nb]; + for (int ix=0; ixgetImage(nnx,nny,nnsx, nnsy); + cprintf(MAGENTA,"Get image!\n"); + cout << nnx << " " << nny << " " << nnsx << " " << nnsy << endl; + // nns=1; + // nnx=npx; + // nny=npy; + // nnx=nnx*nns; + //nny=nny*nns; + dout= new int32_t[nnx*nny]; + if (subframes>0 && insubframe!=subframes && insubframe>0) + subnorm=((double)subframes)/((double)insubframe); + else + subnorm=1.; + for (int ix=0; ix0 && subframes>0) || (subframes<=0) ){ + + if(send_something) { + + // zmqsocket2->SendHeaderData (0, false,SLS_DETECTOR_JSON_HEADER_VERSION , dr, fileindex, 1,1,nnx,nny,nnx*nny*dr/8,acqIndex, frameIndex, fname,acqIndex,0 , packetNumber,bunchId, timestamp, modId,xCoord, yCoord, zCoord,debug, roundRNumber, detType, version, 0,0, 0,&additionalJsonHeader); + + outHeader.data=true; + outHeader.dynamicRange=dr; + outHeader.fileIndex=fileindex; + outHeader.ndetx=1; + outHeader.ndety=1; + outHeader.npixelsx=nnx; + outHeader.npixelsy=nny; + outHeader.imageSize=nnx*nny*dr/8; + outHeader.acqIndex=acqIndex; + outHeader.frameIndex=frameIndex; + outHeader.fname=fname; + outHeader.frameNumber=acqIndex; + outHeader.expLength=expLength; + outHeader.packetNumber=packetNumber; + outHeader.bunchId=bunchId; + outHeader.timestamp=timestamp; + outHeader.modId=modId; + outHeader.row=xCoord; + outHeader.column=yCoord; + outHeader.debug=debug; + outHeader.roundRNumber=roundRNumber; + outHeader.detType=detType; + outHeader.version=version; + + zmqsocket2->SendHeader(0,outHeader); + zmqsocket2->SendData((char*)dout,nnx*nny*dr/8); + cprintf(GREEN, "Sent Data\n"); + } + outHeader.data=false; + zmqsocket2->SendHeader(0,outHeader); + // zmqsocket2->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION); + cprintf(RED, "Sent Dummy\n"); + if (dout) + delete [] dout; + dout=NULL; + + } + } + + mt->clearImage(); + + newFrame=1; + + + //time(&finished); + //cout << "Processing lasted " << difftime(finished,begin) << endl; + + finished = std::chrono::steady_clock::now(); + cout << "Processing lasted " << (finished-begin).count()*0.000001 << " ms" << endl; +#ifdef OPTIMIZE + return 0; +#endif + continue; //continue to not get out + + + } + + //#ifdef NEWZMQ + if (newFrame) { + begin = std::chrono::steady_clock::now(); + + size = zHeader.imageSize;//doc["size"].GetUint(); + + // dynamicRange = zheader.dynamicRange; //doc["bitmode"].GetUint(); + // nPixelsX = zHeader.npixelsx; //doc["shape"][0].GetUint(); + // nPixelsY = zHeader.npixelsy;// doc["shape"][1].GetUint(); + filename = zHeader.fname;//doc["fname"].GetString(); + acqIndex = zHeader.acqIndex; //doc["acqIndex"].GetUint64(); + // frameIndex = zHeader.frameIndex;//doc["fIndex"].GetUint64(); + fileindex = zHeader.fileIndex;//doc["fileIndex"].GetUint64(); + expLength = zHeader.expLength;//doc["expLength"].GetUint(); + packetNumber=zHeader.packetNumber;//doc["packetNumber"].GetUint(); + bunchId=zHeader.bunchId;//doc["bunchId"].GetUint(); + timestamp=zHeader.timestamp;//doc["timestamp"].GetUint(); + modId=zHeader.modId;//doc["modId"].GetUint(); + debug=zHeader.debug;//doc["debug"].GetUint(); + // roundRNumber=r.roundRNumber;//doc["roundRNumber"].GetUint(); + detType=zHeader.detType;//doc["detType"].GetUint(); + version=zHeader.version;//doc["version"].GetUint(); + /*document["bitmode"].GetUint(); zHeader.dynamicRange + +document["fileIndex"].GetUint64(); zHeader.fileIndex + +document["detshape"][0].GetUint(); +zHeader.ndetx + +document["detshape"][1].GetUint(); +zHeader.ndety + +document["shape"][0].GetUint(); +zHeader.npixelsx + +document["shape"][1].GetUint(); +zHeader.npixelsy + +document["size"].GetUint(); zHeader.imageSize + +document["acqIndex"].GetUint64(); zHeader.acqIndex + +document["frameIndex"].GetUint64(); zHeader.frameIndex + +document["fname"].GetString(); zHeader.fname + +document["frameNumber"].GetUint64(); zHeader.frameNumber + +document["expLength"].GetUint(); zHeader.expLength + +document["packetNumber"].GetUint(); zHeader.packetNumber + +document["bunchId"].GetUint64(); zHeader.bunchId + +document["timestamp"].GetUint64(); zHeader.timestamp + +document["modId"].GetUint(); zHeader.modId + +document["row"].GetUint(); zHeader.row + +document["column"].GetUint(); zHeader.column + +document["reserved"].GetUint(); zHeader.reserved + +document["debug"].GetUint(); zHeader.debug + +document["roundRNumber"].GetUint(); zHeader.roundRNumber + +document["detType"].GetUint(); zHeader.detType + +document["version"].GetUint(); zHeader.version + +document["flippedDataX"].GetUint(); zHeader.flippedDataX + +document["quad"].GetUint(); zHeader.quad + +document["completeImage"].GetUint(); zHeader.completeImage + */ + //dataSize=size; + + //strcpy(fname,filename.c_str()); + fname=filename; + // cprintf(BLUE, "Header Info:\n" + // "size: %u\n" + // "multisize: %u\n" + // "dynamicRange: %u\n" + // "nPixelsX: %u\n" + // "nPixelsY: %u\n" + // "currentFileName: %s\n" + // "currentAcquisitionIndex: %lu\n" + // "currentFrameIndex: %lu\n" + // "currentFileIndex: %lu\n" + // "currentSubFrameIndex: %u\n" + // "xCoordX: %u\n" + // "yCoordY: %u\n" + // "zCoordZ: %u\n" + // "flippedDataX: %u\n" + // "packetNumber: %u\n" + // "bunchId: %u\n" + // "timestamp: %u\n" + // "modId: %u\n" + // "debug: %u\n" + // "roundRNumber: %u\n" + // "detType: %u\n" + // "version: %u\n", + // size, multisize, dynamicRange, nPixelsX, nPixelsY, + // filename.c_str(), acqIndex, + // frameIndex, fileindex, subFrameIndex, + // xCoord, yCoord,zCoord, + // flippedDataX, packetNumber, bunchId, timestamp, modId, debug, roundRNumber, detType, version); + + addJsonHeader=zHeader.addJsonHeader; + + /* Analog detector commands */ + //isPedestal=0; + //isFlat=0; + rms=0; + fMode=eFrame; + frameMode_s="frame"; + cprintf(MAGENTA, "Frame mode: "); + // if (doc.HasMember("frameMode")) { + if (addJsonHeader.find("frameMode")!= addJsonHeader.end()) { + // if (doc["frameMode"].IsString()) { + frameMode_s=addJsonHeader.at("frameMode");//doc["frameMode"].GetString(); + if (frameMode_s == "pedestal"){ + fMode=ePedestal; + //isPedestal=1; + } else if (frameMode_s == "newPedestal"){ + mt->newDataSet(); //resets pedestal + // cprintf(MAGENTA, "Resetting pedestal\n"); + fMode=ePedestal; + //isPedestal=1; + } else if (frameMode_s == "variance"){ + mt->newDataSet(); //resets pedestal + // cprintf(MAGENTA, "Resetting pedestal\n"); + fMode=ePedestal; + //isPedestal=1; + rms=1; + } +#ifdef INTERP + else if (frameMode_s == "flatfield") { + fMode=eFlat; + //isFlat=1; + } else if (frameMode_s == "newFlatfield") { + mt->resetFlatField(); + //isFlat=1; + cprintf(MAGENTA, "Resetting flatfield\n"); + fMode=eFlat; + } + //#endif + else { + fMode=eFrame; + //isPedestal=0; + //isFlat=0; + fMode=eFrame; + frameMode_s="frame"; + } + //} + } + cprintf(MAGENTA, "%s\n" , frameMode_s.c_str()); + mt->setFrameMode(fMode); + + // threshold=0; + cprintf(MAGENTA, "Threshold: "); + if (addJsonHeader.find("threshold")!= addJsonHeader.end()) { + istringstream(addJsonHeader.at("threshold")) >>threshold; + // threshold=atoi(addJsonHeader.at("threshold").c_str());//doc["frameMode"].GetString(); + } + //if (doc.HasMember("threshold")) { + //if (doc["threshold"].IsInt()) { + // threshold=doc["threshold"].GetInt(); + mt->setThreshold(threshold); + // } + // } + cprintf(MAGENTA, "%d\n", threshold); + + xmin=0; + xmax=npx; + ymin=0; + ymax=npy; + cprintf(MAGENTA, "ROI: "); + + if (addJsonHeader.find("roi")!= addJsonHeader.end()) { + istringstream(addJsonHeader.at("roi")) >> xmin >> xmax >> ymin >> ymax ; + // if (doc.HasMember("roi")) { + //if (doc["roi"].IsArray()) { + // if (doc["roi"].Size() > 0 ) + // if (doc["roi"][0].IsInt()) + // xmin=doc["roi"][0].GetInt(); + + // if (doc["roi"].Size() > 1 ) + // if (doc["roi"][1].IsInt()) + // xmax=doc["roi"][1].GetInt(); + + // if (doc["roi"].Size() > 2 ) + // if (doc["roi"][2].IsInt()) + // ymin=doc["roi"][2].GetInt(); + + // if (doc["roi"].Size() > 3 ) + // if (doc["roi"][3].IsInt()) + // ymax=doc["roi"][3].GetInt(); + // } + } + + cprintf(MAGENTA, "%d %d %d %d\n", xmin, xmax, ymin, ymax); + mt->setROI(xmin, xmax, ymin, ymax); + if (addJsonHeader.find("dynamicRange")!= addJsonHeader.end()) { + istringstream(addJsonHeader.at("dynamicRange")) >> dr ; + dr=32; + } + // if (doc.HasMember("dynamicRange")) { + // dr=doc["dynamicRange"].GetUint(); + // dr=32; + // } + + dMode=eAnalog; + detectorMode_s="analog"; + cprintf(MAGENTA, "Detector mode: "); + if (addJsonHeader.find("detectorMode")!= addJsonHeader.end()) {; + //if (doc.HasMember("detectorMode")) { + //if (doc["detectorMode"].IsString()) { + detectorMode_s=addJsonHeader.at("detectorMode");//=doc["detectorMode"].GetString(); +#ifdef INTERP + if (detectorMode_s == "interpolating"){ + dMode=eInterpolating; + mt->setInterpolation(interp); + } else +#endif + if (detectorMode_s == "counting"){ + dMode=ePhotonCounting; +#ifdef INTERP + mt->setInterpolation(NULL); +#endif + } else { + dMode=eAnalog; +#ifdef INTERP + mt->setInterpolation(NULL); +#endif + } + // } + + } + + mt->setDetectorMode(dMode); + cprintf(MAGENTA, "%s\n" , detectorMode_s.c_str()); + + // cout << "done " << endl; + + // /* Single Photon Detector commands */ + // nSigma=5; + // if (doc.HasMember("nSigma")) { + // if (doc["nSigma"].IsInt()) + // nSigma=doc["nSigma"].GetInt(); + // mt->setNSigma(nSigma); + // } + + // emin=-1; + // emax=-1; + // if (doc.HasMember("energyRange")) { + // if (doc["energyRange"].IsArray()) { + // if (doc["energyRange"].Size() > 0 ) + // if (doc["energyRange"][0].IsInt()) + // emin=doc["energyRange"][0].GetInt(); + + // if (doc["energyRange"].Size() > 1 ) + // if (doc["energyRange"][1].IsInt()) + // emax=doc["energyRange"][1].GetUint(); + // } + // } + // if (doc.HasMember("eMin")) { + // if (doc["eMin"][1].IsInt()) + // emin=doc["eMin"].GetInt(); + // } + // if (doc.HasMember("eMax")) { + // if (doc["eMax"][1].IsInt()) + // emin=doc["eMax"].GetInt(); + // } + // mt->setEnergyRange(emin,emax); + + // /* interpolating detector commands */ + + // if (doc.HasMember("nSubPixels")) { + // if (doc["nSubPixels"].IsUint()) + // nSubPixels=doc["nSubPixels"].GetUint(); + // mt->setNSubPixels(nSubPixels); + // } + + // threshold=0; + // cprintf(MAGENTA, "Subframes: "); + // subframes=0; + // //isubframe=0; + // insubframe=0; + // subnorm=1; + // f0=0; + // nnsubframe=0; + // if (doc.HasMember("subframes")) { + // if (doc["subframes"].IsInt()) { + // subframes=doc["subframes"].GetInt(); + // } + // } + // cprintf(MAGENTA, "%ld\n", subframes); + + + newFrame=0; + /* zmqsocket->CloseHeaderMessage();*/ + } +#endif + + // cout << "file" << endl; + // cout << "data " << endl; + if (of==NULL) { +#ifdef WRITE_QUAD + sprintf(ofname,"%s_%ld.clust2",filename.c_str(),fileindex); +#endif +#ifndef WRITE_QUAD + sprintf(ofname,"%s_%ld.clust",filename.c_str(),fileindex); +#endif + of=fopen(ofname,"w"); + if (of) { + mt->setFilePointer(of); + }else { + cout << "Could not open "<< ofname << " for writing " << endl; + mt->setFilePointer(NULL); + } + } + + + + // cout << "data" << endl; + // get data + // acqIndex = doc["acqIndex"].GetUint64(); + + frameIndex = zHeader.frameIndex;////doc["fIndex"].GetUint64(); + + // subFrameIndex = doc["expLength"].GetUint(); + + // bunchId=doc["bunchId"].GetUint(); + // timestamp=doc["timestamp"].GetUint(); + packetNumber=zHeader.packetNumber; //doc["packetNumber"].GetUint(); + // cout << acqIndex << " " << frameIndex << " " << subFrameIndex << " "<< bunchId << " " << timestamp << " " << packetNumber << endl; + //cprintf(GREEN, "frame\n"); + if (packetNumber>=40) { + //*((int*)buff)=frameIndex; + if (insubframe==0) f0=frameIndex; + memcpy(buff,&frameIndex,sizeof(int)); + //length = + zmqsocket->ReceiveData(0, buff+sizeof(int), size); + mt->pushData(buff); + mt->nextThread(); + mt->popFree(buff); + insubframe++; + nsubframes=frameIndex+1-f0; + } else { + cprintf(RED, "Incomplete frame: received only %d packet\n", packetNumber); + //length = + zmqsocket->ReceiveData(0, dummybuff, size); + } + + + + + if (subframes>0 && insubframe>=subframes && fMode==eFrame) { + while (mt->isBusy()) {;}//wait until all data are processed from the queues + detimage=mt->getImage(nnx,nny,nnsx, nnsy); + cprintf(MAGENTA,"Get image!\n"); + dout= new int32_t[nnx*nny]; + doutf= new float[nnx*nny]; + if (subframes>0 && insubframe!=subframes && insubframe>0) + subnorm=((double)subframes)/((double)insubframe); + else + subnorm=1.; + for (int ix=0; ixSendHeaderData (0, false,SLS_DETECTOR_JSON_HEADER_VERSION , dr, fileindex, 1,1,nnx,nny,nnx*nny*dr/8,acqIndex, frameIndex, fname,acqIndex,0 , packetNumber,bunchId, timestamp, modId,xCoord, yCoord, zCoord,debug, roundRNumber, detType, version, 0,0, 0,&additionalJsonHeader); + zHeader.data = true; + zmqsocket2->SendHeader(0,zHeader); + zmqsocket2->SendData((char*)dout,nnx*nny*dr/8); + cprintf(GREEN, "Sent subdata\n"); + + + if (dout) + delete [] dout; + dout=NULL; + + if (doutf) + delete [] doutf; + doutf=NULL; + + mt->clearImage(); + + } + + + + + + + + + + + + + + + iframe++; + + } // exiting infinite loop + + + + delete zmqsocket; + if (send) + delete zmqsocket2; + + + cout<<"Goodbye"<< endl; + return 0; +} + diff --git a/slsDetectorCalibration/multiThreadedAnalogDetector.h b/slsDetectorCalibration/multiThreadedAnalogDetector.h index dcdab30f5..4e0fb1b90 100644 --- a/slsDetectorCalibration/multiThreadedAnalogDetector.h +++ b/slsDetectorCalibration/multiThreadedAnalogDetector.h @@ -17,7 +17,7 @@ #include #include "analogDetector.h" -#include "sls/CircularFifo.h" +#include "circularFifo.h" #include #include #include