From e65e7ac42f49e54cb98f70a2d854b6d5943f756e Mon Sep 17 00:00:00 2001 From: Anna Bergamaschi Date: Mon, 29 Mar 2021 15:35:27 +0200 Subject: [PATCH] energyCalibration updated with Sophie's version --- libs/pybind11 | 2 +- .../jungfrauHighZSingleChipData.h | 238 ++++ slsDetectorCalibration/energyCalibration.cpp | 798 +++++++++++++ slsDetectorCalibration/energyCalibration.h | 31 +- .../Makefile.cluster_finder | 47 + .../Makefile.cluster_finder~ | 47 + .../jungfrauExecutables/Makefile.zmq | 23 + .../jungfrauExecutables/jungfrauClusterFinder | Bin 0 -> 69248 bytes .../jungfrauClusterFinder.cpp | 161 +++ .../jungfrauClusterFinder.cpp~ | 237 ++++ .../jungfrauInterpolation.cpp | 250 ++++ .../jungfrauPhotonCounter.cpp | 453 ++++++++ .../jungfrauZmqProcess.cpp | 1026 +++++++++++++++++ .../multiThreadedAnalogDetector.h | 2 +- 14 files changed, 3301 insertions(+), 14 deletions(-) create mode 100644 slsDetectorCalibration/dataStructures/jungfrauHighZSingleChipData.h create mode 100644 slsDetectorCalibration/energyCalibration.cpp create mode 100644 slsDetectorCalibration/jungfrauExecutables/Makefile.cluster_finder create mode 100644 slsDetectorCalibration/jungfrauExecutables/Makefile.cluster_finder~ create mode 100644 slsDetectorCalibration/jungfrauExecutables/Makefile.zmq create mode 100755 slsDetectorCalibration/jungfrauExecutables/jungfrauClusterFinder create mode 100644 slsDetectorCalibration/jungfrauExecutables/jungfrauClusterFinder.cpp create mode 100644 slsDetectorCalibration/jungfrauExecutables/jungfrauClusterFinder.cpp~ create mode 100644 slsDetectorCalibration/jungfrauExecutables/jungfrauInterpolation.cpp create mode 100644 slsDetectorCalibration/jungfrauExecutables/jungfrauPhotonCounter.cpp create mode 100644 slsDetectorCalibration/jungfrauExecutables/jungfrauZmqProcess.cpp 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; ienCy3FcW%`!Hfb1jS+Bwz)ApNNYL;y z12Ju9BI{$9yXz`4-6D4y~55={I9_fnlegbH)I&5 z8SI9ch8cz_po+tPjn1yVX>eFgra?1k>50XE6aHt!P>7(NKorhkAp=U_k5)QCDsS|$ ztFIhxet`!yXfRkAg3=MMoiDN=U?F?apof^jV<#f^SlM4%RL>AfNDYLYdB;Am&0aIT8Mkpq&Vx3Wg`b-#AJBHdJgP`gN0(+d4_P?*nhe|MAb& zKv>5k5jIUyZuTVQo;3-5-z50!CdudgN#w9{5`5|;_}NgTiR4x}34Zz{aoQtn?Sk*8~t`VN8qQo}_CLvJBUG8ppNUxVRu2uk8j8Udm3};a3u4bToNd^8fqe}90Fy(4VdtY`@gJyi zvm}khUkN^^gP%f`e+%g7Ig0<(uLsA`Pgmu3K=@r)&8=YlA>Z7G*b-)?%eEx z;*y$@mF2bmlA7G?%!;bYlH8)ZDoP9mi&q!iqDEv^6xG(2)N=X>scTDC6|5}ruP&;n zD=|=@x~kT1Pzx(qRp&3MFQ};U-Ob{POKSZ!Rrf}tec;Gn5)J;Ss=SiXl~FZD*U}^Nsj+Ipa0(K61UHlc5Oj+!J6{olF9-q)3uria;mN?t*EOlW1$>>+Vuqm zkX_Z9g1d@}3w%YjexfQa@-kxNmNu^dQMKj1g3|H|%Juu@wFOnxC6&4iAQaSf6-ucs zsqwpgNR{R*E2=5**Px2ENaVT&xNCVnDz!VeAUCzZm2-W8cS+8*1ylzxX9LS>Hr^$x zDP2`AYK$m*GxOZ7Ks*sOkiDzs^@tQdHwY%siGC^7E5suv%$S z`A6&^gR0sRR^QxPZ%Rug1*k8&t9*4@+V}BJYFt2V4#|~Pxkz2=>kFz&YHF)0iz>?f z_ZF;1r?@@GpYE%w^Rr^BDr#Njc&l7PT9f0isC7@k_ympS%SDr2hbHr}eEvD%t8E1a zN-B#hCX~?xMMd?4b9`-eo!>>JjuXK=h@-ZA0x6(c6vwr@qNH-AzpS97rlzVU%KTW3 zQdOk50=kb5SmUS8mqxGt^3qZYTxZC2FIq&Zz+Q4n{ENylA^cSp247JH%nF4d2Gx8(0qbFYR+X%( zE%6(C{&Hkf%Jj@oS6fk1Qf(+*L#z_Ts;brCs(?ySjnJwh7_)+s;v#<$BJP6zG2;YO zWw9a4<6f9qFmLWW2rMfzv*23TWS~fP1?hOt%+Gf$bQh$}y_S>a=Yu4Du5C;xElnRJ zYafGEjtTz>#<7(Tt$Sjzrku*KF)Rjy1dESDh}K6jprv$_Dh8B17oHErP&^Nr809qJ z&SDsfGvl8z`j1K@41WlJ0An7l&#d@gcTxE?1JeE&<0d^6|FJIKd*!`n8?f5{g)Y88 z7vFyD^f|O5|BWu5hQzV?n{@e`46mv2UclG3+>MlG!|(Bz(n!s-AcCoe7u0ak@J|t* zV%V#O{j(lMRbmZqAxsbP8>_}qnAayJTHf&XLPE1OJSsZ+sN*+fLu*Ry(Z+ScXKF}5 zZC*n71v-3EsE~!G;;%XtU}38cPi?9_DLQ?+1g{*;q~i&QHLj< zw8yK%Yq~&j`8vE_K7~3w6d?MP>G1JUtie#N!$TpWPrVMG5XC~CI=rC6Z_wdMAGD`k zho?22_H^j*`gK*O4u76T47gc`Cm)3NY}MfxYDmB?9sU9xKBU86sKfW^@JTv+zYhN$ z9sZyWf3XffsKd|D;T0YJyE^<49sUv>eprW}slyw33R#GB#LC|IAF~dBsSa<^;V;wS zlXQ5e4sX@rFW2EybogW)K2?XmLWj@L;jh%;?K=GTba+vRPtoDMI{YjhK3|8Qt-}}U z@K@>ZWjg#E9llzJC*O0HX}^kDpJkD$0%&3a3KpH zqVO~dds+B>3ZF$`I}5)}VKTXqR2JS(VKTK5D+}+XFqznhg@t=5Or|wrVBu{PCX*UD z@+}ZoJWF9RrIA4veu~0mLL>bw{7VXx>5PO}_-7O*lNs5{!VgoJOl72#g&(Bw`4nzv z;dK-y(->)C;rl2|CNWaY!Zj2oQy3{^;U7_$Okl*z!go;^4gu7kg>R=YnZiga3*SOv zGJz2*3ooWHnZAgHg%?toOkTvm!Z%QuOkL#2H&p+56t+@$kcDScm`q!wpM@`@@MRPZ zvG7F{CX*G}%EISTm`qiqlZDTsFqx=GI}689m`qcog@sT3JHlj=BGoK>jKV1tE@a_D z6rM$4FAKjQJ73aWGf3lOkpw!kxmwVkiujNBJC`^ zj>2RDA}uU@ABCywN2*!4hQie4BZVydBMMVjk9ZOGit|+G%MqpH;{Bg0*cSH#5)#J5l)dBb4RGcsnLo z;wl438+ImJfn+aiNC$#-J;z@Z*a?Kp;9w+dO6~)dD=qXFvM$@mgm|Qnk+ECat>pd) z?EX9P^EUcAl{OD$xaBDx>Ewtvh4{5>-$L=*cs$sx{u+D*Ur268PRa+wq?`%kwGk6~j zMk#Kfj1K+)WG#4OBHTs-^GIMbEJ`~RrHaar?Ix@!y*J{uAVigTZ-m|$%1$rA=aO4A zn!@cTh`9!}l-E%!r|mN#cs5M6Bj_5DDw1=lDXAwKW+`c(iGh%jgjFa8Ta!CkHqPW! z3Ybwa#pDY!4RUUBZftIH3RS0n4e{kqwx@-ZbgDNZMKQRMM}VLc!$|ijumOVeK;V%Z zAVeL(R`mOB(@%)NBb`zXlW26&A?P#k-4;fn^s{JHraP#aGyz}v6#{J`w4~7bJ}O9Z zT?c9pDRl34TEr;lG7jU`{@Mzr?TSagXnkV}WYW@@VGvq9Bp%^Uu1ul( z2v2u|st?Cl>5sBX$!CdyUQ_V?k#E2K_T<2&eOI*L@kgW{mZT?Q`E5&!!5}z$#Ngr> zQE(0zLQJZVhZuM}R-DxbvSOGuEHyOs zdQom=P-wjbGZ3{ki7r;H)Sl2Q z3q(d{>jw~zQ##?14l6l`hhY_Errp`{GgFR$F{a*oN$S)9{!2CC!^$c}BsYluJIExG zh`Q7Kkh0o9S;;O#HxPtzQ4Td3~1C^saRDQ~ePu{0icuO)4KE}3;k8W}P) zdh{dJQ2vba)dJrfXN53Bc^=9k7H<1;tf7n`j4Ivo z)0fD_u$r1aWlRwN3v^X2K52q@8%m_4|3TAS-=YKdkO&w1H|Ip0haA25nbwm(XP>rAreG)DhT`4D_| zS3@mGq6bLIfCWg=1t}XM&(I)Gy+#=#siVo5l?Khzl<+Dw;S*}YE0GW_$;WtRpLjMj zLsq0?l=}{?4y-*e2T+JVr2>)4O8Rpb{*<;W%Ng}9BIsnVLqiDcAOq0O!tW_bR2rc^ zQ5scL98khcCs@V$a{F653JuD`h}xz+YBzLKa64@T(mFHp6!4=HR7U?c(Np=Ho`R$F z_W_^wnOhpTQkSAV!1V9W7Cy3HV zYhqUH=rZ7;1@9Ltwo}x|pxQsXD5TN%3}bTi5oQ2l&!HcR#(r^NxDJW8@x=R)7{f46 z|Gt{O<;HEe{P*Xng0fwp|?k+TiVXo z`%__(uNe6uM1U4JCd>|jw^phZ@|dDqE-?BOEGaM<2R|bi&DczxMlei|-yiW-BZm}S z3gn+ON-s)(=Ojs*{ga!K?uhLx;ij*o8+IW7@YL#mu^Xc0mkP^5;Vc@v4abD(Sr+<} zq|0BsGcatFvf}-XFl0+Xm*00xxGe;jm4rlup_FLi1*l+usgQzrO99Y~ z=%r)Aoo^o#?f_)8mZeCQ-vjJkU;!yH=BxC=(ts)Ii9BLdqS1IB(GLFFLQT_ABo#+V z@9d#`s9pryPBg_!Ys}J-9Wa6DA99viSG!V5{}01$C+ZVLB$BGjq%{^~qw1ZUWx=Z* z?U!W+EA~P)|55r1X-%S@K2fD_#u8{MX0$2P2asM3x6tBiWH1`GCykA_TgQf5NS}4| zEyLr&w35-scZ?0&tHxb69BrM32`Tr)Z|K{=a=d~a8Cq(t5C5h$7ybIEo zj$F5N8a@WZ*)kzn^N4tisYo{LCpM#qTFbJDM{<1g z$9(SOl#V#=)U}x-7hP)~&phxciG>Pa0!J6`8S$EVQM@6Vt{yFo_ctmJZD2}{mU!XF zEoNTi2WsPx)G%>l{bS<#c#|fn5c<#f`hb>I=%X#hX07(91G9nqC z;7W@_dcz^@b4lT}cUU{Rt; zDnIl~Z5VE*E)#=Q7Rs3HMx#M55-Vp~VS*UVp?ysHD<)0PI*zX1kR|UgW77YF^1c5| z`Trn&{$%O%+|oD9{RPG6Qb^~uaHW07?4M2!7$&P7({?&Ysvnc)JEcJ`*d3rA`i$s( z546MgBFYV6QbgF%!jVCpyv~d-8xt?H-hE8H$D4sAo8RNWB7}74deVP1anxQ>JE!gc3k~LQZm2Epk*C69kS*d; zdmD94AuB|ymmOFOCT~U^hn_MqBrAgUx{85)2|6>=3 zzUQ33)Y5n4@;>IbJ!7MiR(3%TB1$Yv*_+g>+!~3Vr+FU9JCAX|Ge-P0qCr z5XVZ&&;aDsb^*MBETbg`wwWn=A-o2KlrLxy+JJvWw@_M9dFIb41FD5pn$zran zEmS|nIp3bxbe>bTLy!9!r#dBjt|aCk-_ee@{Gp-IN=sty&?_iM@_L815UlnL{d_dK zJx5-ei1}kUYYv;gWi>!%hPFLTW%4phTBx+WZGSkko^VZ`w2>iGyPVk|$~j45hv*v- zdp?O1V;A%?y|VBayRaqG5}R2sc7#pqD4m>HFXvc-^+ww%X>EP;c6mu6D521{J}I-I z;Yo6fN|a*4>=uP1%#Hy7Eo{k3+@XutP@|P4 zJBFlTwUVKx9IE6V);_3b7O!WzD6g?&I5ac7ytF~sGJvT)>a(8JC#xQ*i2q(GtG@X? zWC)r-kI~j6Ep2G{v%c{mZ|aY4z@CWGmo!gj7JON*ye8jvOj_0;1UDj2%>CYRNH5Dv zz0le_l2I`D%=D0e>KlvEU+dzB?$hRRISuj>^qSfv(Kk#DAZ*zS-96AV;_K3bc z5acn3V99-jI*IcYWS@dL>vRnMBVIc-#}Y5G--Fznx5~Boa!x~_*9d-c8X7i?jt#t@ zA(A}cqVE{W6Sm}7s7&PR!Bir40R|aZ6@kPGYF+KeJjeT#=V%qq4HGp8-|fv z;ghcqdW~|X7a7?4n|Dy?FVMne* zdYxAa)?NCVTqJ?DR)c?OLD=8jaU$Zq(gF!H)?OGo%_)e4ReQzslFch zY^m2J4T|Jhe$O%T&;@I8B$q>Hb++={b$Xd#hM-VTpn%r zSYDTB9oMpdAJ?)yTFYYoK`qK=-GNW7ZTp1{ACmwv&`4iXjICD)z79;=9`rb2 z!ydfOUz#WcLwFJ8WvL>Pv>re}>Tw8LUN%P1p`c%hdod983-WIm*;)+rG)r$A5n7*T z7@A?`CS!>>dr68@`a42c0857;Rn}K5PK>u0=VG1G`_kuFV)j9ZPKlP0pd0b#b1M90 zncUGw>r9aIcqC%ko4GXFx1vZ!iiHF^BCR2NX>kf6SXrw=PldKTv?k)uk12zHNJN0y zGr}m?wToQOFz{z?*~v18Aw)S~%c!u<91mA-Eb$5sVK_L4R3PROSOGH*Bex-w(auBd zhV`T$(gt0egko&F8}G!BcKiHgiH&(={6u*vWW2{;AZ+P19^obRK(^;B0rh=T>P)l{ zhtWrm-!J;U<<^E+fVJsoI+%jIDU3-hEuP73!si2VNytK*M}X)Y1Hwap7|S0`?WEVh zi3M!k!t~-BTW=#+fw4j*`!j>M@I+}J5?S?-?s~X0AQ@?E5I>x-%p@YpCRJMS!k+TK;iPbMU zk6`OD^3Mo=4xunJ7}T4nap^Id^7PEkK`+ogbLOwHG=2}&znog-k|;+^*s7ROw3n@~ zyU@ZTIv%MzPhuBTi=ZF9slH?CW7hc5>Vv&z=@pvve*w!wOJ!~>*i2QVy|Wz<;us71 zYeHmVuxP={sNGl`X+#j^#kwSg?E}9z(ykCiCPr{|P9g@KZRk3)AvN;l%;}dba!x&R zBio8GkVYQZ*8R-3K4)8IZrTFd`ZU?re~a4IQ=_)^$6BH>wpG|z3ZYTsh!R^MVUjG$ zciSWBxHhho0oDeq9Wgk>n ztu1rTYwx1VcG4-W#fT#iwKVKXWBzc5?K7uzcYR~JOFBgT)umSULMPSML6ggPC)OGb zO-@^mrFp&sHFVnEa6&N6uQ@Q3xNN(fn7B3niR$N+hWR2u-NLXNC{@#5%4U~yer96h zcJK#7=alwqZ9>~dPDtHh3_H;_x&nlUDj}eue<8kLi&NabKDQ5JTH;>Wv<^ELZRLm_ zL;nLHA|&j#&zy3VN08oy`y4Wo@AeAPTkQRWQ(l`ZNN*tE#HbN@nXqXiksFZ;{Q=sa zD*8^KBT6xN&!88hPTw80G7tw2u_+u|fv3<^k&S&#K3Z!@uVeWXb$E!_*V4S7%tnlO zVaiO&5XhYQE_37eL>~+&I^QMp^Wq!7Bg)x!RB)iwA#{7bYbcfC7pm(8EWI4Icbm6_ zT_PF~e7eZWq^N7I==@>)d}Msskt?>H&<1s(HNZMVVmOn}^FvLm`yPB%Gu|HE+)avQ*tZ&hbuNjQH&s(bPL9CPCZG9q9c6Xq|8Z4})! zqW(5j{8->~7o|5b1<=u;xuipJb>-r zwP;S33buUHL+IS_VT!@KEwD;GN(?5J@jZhvSTK8h3e7SDp%{11M})_!+vEEL3!oUO zA$f+o=ffBgSW84g37Y@RZr+W<2guYB_;!l0z65V1+_n?`tI|SgCyonM@iEW?K$Edl z__pV`&}@oMYoTnKwa{pqO~DyW7a@2lf|pS+sp)Juvyv^c$!0X& z;0T17W7cU~8!v3)bKyo9E0}~77knF1Pt&`i{0PZN4B9QCtyc^#wAh75c8j(i;gJxH zXKge&Wqkh>6#?I83+sb;qbV~Q@=BF&lacYQvTFmoYa)zYp>;JfV5~B}CDUwSlTSq@ zG2SK9Wx^&eP>gkrX;eWn-N0CX$i*UAvy7Tqo!-J1!t_HJI@*k5rQoef~P~%z=ka*z_qR> z{|sG7XnlaRX8|%pd(&L(Jw%Ywl75M@%y_A^C+a+wA zNsS9t!h9^OY$HL*Oi*ZPzF-<{A*vpE8UJzh`l8kIuMUiyzxp0jYhX;Z_N&!0QMITp zRHHk4$JN7OJG$m=)I(^^28+?U(7Fb-nr<41Tn9Yq7s)u1!U8Q|4?bG1ZzQpPMo;9B z4Tj$4Ly>QpTB2RY)6qiF1VN1we-s&D#I2OAf%Y@pFS1vRMVO9;ccM^o3$V!#S6VoI znBKBk`b65MxP`3KO3Qvy$sLmu6flC_vJ2yurjd0P(n?p_XS{xH`3Gc8_;wYX9@1{1 zwVu>69%_lTuZ#4O*6VG1giQscmoPJOBlK|x^v`4a7)AyY{sZP-p|yn+@vEuA`ZZBS z>`Mzd1{4w!DPlY=t4&iNL{i9DYPKn)h@M#dy$u!sBMJ7A(TEgzA<$$wownYZu+z53 z|J`i(3P>6EEe@KnaZDS@Y-NsYr^B$3M()O5(l^4!DUyq#y6QjyN2{aNgw=s5A@T3z{5x#@Yu*I^LhEhd zDq07Xf3AoxY4c@>HZl)n8h@HyB*E>DCV_hGARB4xY5q783vx6jug^sEOec4cPHKt% z9@(ebRT}@lb>u-GXyq<1!%T5%>d)=ypUp|wXee4Nn-1R;O#X~Cmb$o8l=I;i{UA;> z?htJU{g;b@AuF5b(>YdA`YVM;EN*#G5|Ua)n^O0o7_1*welB&e0kvQU$CQfM$>#b>w2_hxHb09;LklZoM7wiiDc(;PIHp!w0e; ziDx?(832P`rD0}p45l3#<}8k}k79=VX&JN=1Md#lx03kgr?6|I=J+n+5VG?8_wi5) z4~2M$RU%G#mWSYtq0#2oJj8qtXH%Ub4+B%$POrM@6Z~uB?<8mTxJ_|%bU!?>^=iqv zNG;cz;RF@gaTtRkG}9^E^pQF-lJ2%RC10Yj`TQDESuUgwtY;^D}RaV+DHuo!oy!=PyCb4qAGbJESRCsW`$eD$2mNxfmhzj-N+~ zJogH>%V=e|R`#HQuY{R_C1yuRa9RW%5-LsIpH7V2roNNk>}K+Jh3>;D6C2WAS=}C@fa%kIAcb1&9fOfKk$7rIBJ^X0@3TI`Ie9$^L ztR9(*dXJ86G>)+R!&8s^i`_7WREcFGTFX?0JA2LWf-vpCgpX~_z|E9dYR@TiCDU*t5NK0php+N)gN4FL@&YuGlrc>r*hi%XzS8> z+0y%-*c|)P;2l@c>3`m@j%Nz9P3v4i(TKO{rhpsP$7zG#VM%@S2U2YVMq#({v~73Y zpI9p&R<1?cVrzxe)=UQhe}x8xEMbsUBh2+;Z2h7u5Vkq|kL$PYC00Lhh-J!{Dy|_v zT(mxSV4ZzVwEf5>!dKg;&@=k7eIZ(BFOFRhG7ew`?g>@|(YUnY#*h`FRq7H;sc?HA zt+@w|WBrtEpp|zQB+I=kG*Rs*k4?iywmv$n*kGox3}Q9MN=W%F8e(kCM^%RQYcTNR z{7m4}d~L#NCB=k>f(x9s&&3@Ggi~Q&m3#iOmac8Y}Q ztJDi^wgYRm^#Mam%sx1THB}lEwDr>&{(w1q$tK`oT`jJ?7mJ$Caxp)~}8z%VAFf5wYL3JJ^hJFwzP+mm);VQ>LcP9dcinj)NQ zIvZpE*=FHFOr>w{6K*#p)1beY2K~jD3!LX%45?q}T%05ZdN2XlU9;OZ?Eiu2d-b?r zNtcD&Phd1l-y2+t5s+LydjjDY#LSUx!@}*KA#TpzCJTHiz$HH~-2Ndj$*07&cW6F< z^ODdU=10OZ^~6Yf3%Eu5y-k8Z>4zmZ-N9XRy1kUL}BZ!9v~Sc;ZP#x(wf}SJM^lUX)XYDAZ>+qAyy|bl zuD>I#fq9?%7HR@Rj+1`1OL0K&a+>GBcoQ0D1j)wVhWGq~ASs|zfshkW1JQGeG!eWF zlkV9|kuj$BLOuqH>C13#5&Sk`e@UAujy=oB?JLcD&=zd-gzceaU~DkM^GGLuUZy$E ze8}qi^S$QAd1$V=^TqhaS@7SdK<&tPDVzFaQ#+1#Cqjvm_;h7-ykPFx_f-c^8HTa9 zH6>|8Q-kFkbw4W?4k&n7K4)86d>15^Z5gq&s18=^x@{sF=-9T6y5|MG9$_0-$leM8 zfLi*THk{ZtN6IK`wCe!*;d<3*?zBn8wxzJGea9~7U@S7nMmaAJFW8VG;p1&+)2Q)c z@?#G5oIrCKE+2g@w324BGOO(TG`r?B!Y+5NoSBbpaoVQc1TwkSi}RM)m&DF^eS>YN zit=OrrjJ=|$dQfhaoQA*&8df|jn2n_s7E?=JAo4N?SxEo{Vto=+!R6!yc@G}MAX5$*npxHzW2PCGl8$FnkVB57H9Gs!Q!iLR2 zt14rmzH~SA1Ls{58}H(quy^Bp$4V#>T)H?rz&jLti&AXYr~Mf;+*|CE`5hMR|K>ZWUd zh65L3f}f@Cd*RFlRIlyD_@n#>7Sl4N@b4pyX9qdCJ;(k~DoEMfE!dcQ$ z!HwJKm3lCrp`8#>%CV>IKUtR|)_eN=z|`7Q2UR(l#EW?EKS-0yW2E_(ja?&PMi+_{PZ{*Cx~p6 zL$!Wq&d+^RBYDF4xi(#=`p$qDMsve@X`KE05Naqdt-uca`AlC@p>cuCIYwq)*g6^8 z&i}|ZHV%+ojP7^+kFk%@N@)Fw8_JP(GLq3q$YYxYi0&h}2!g(hkyND9ZyCVx8)tHv zsANHi?1C3A1ux%DnZDsPE1C9%m6=Q@GtK-` zoAEI_FG7^dtjY%=wG6CBl=YJrGJ>pmA+*lY6%wNsvJDx*{Ft}VHwC(<>k7$Gn)HQS zI(Z>KQw#Ce>k9Z+;vNpu+o*sbPU%T+rJdLo%E8r{pCGH{ZJ`yyqapPiG%X~yb`vM+ zy*?mdiX(?qKMuvWPY{322%m>ieDwtJ)reP4zDPtZVDi7{U)cvf$NE=Z;252M`Uj9T-rRPvDOUn%4vJ74$Q!b7Tm#mPfD|H?dICi1WRSNunx zbBVz?Yu=)FIYw8%y&R*fUl%Zw)$eget<`T652^KQ;UTqtH9Vx&ub78e{o<5c5yIBm zVe%tA_=4Ud{nz|s2`pE=f9y0=0WGSl-yx3C)$brMlhyA9My=KFKX^#3UndW#^?R6y z)cUpZkXpah2x;|u>iIGC`+w!HJPxr=?yr1<$EsBcBg7~ExR&R*Gjej!GP{QTlV0c`QATX4Jt<@@SIN^=&#Kc8(~$s+^-` zg`MuU{cc<{Bj>Lh8<0s&WyB9A&nIsR&iN9IGnTW)`G1oq^Z#Ny6=_q*0sH~u`j==& zmN8?-5nM#xdgA7B-k;a0@v2N~5nR3k)}G(sjr{5_BvZD3LSrtw-11+N=VdJaiR5WU z%zrM=_tB`DxF?n8iy&f=|M%s23*!e3HHkdusqsvnqdxVEf#$kO`}W%W1%vctTc-T` z6!mF&Q2yi*x{oU1$rCPMwLqFdSQL znDt^`LUShV!hRClt9}g>fstv#`n`Cetv~FpJ9;89i-IQFyrzYaDH93==*fh$=mKXZ z97h*;GT|A%9o~c6_XrpCz>ySZfFCJN-ALM3)918>{i53!Lb^GaM&JI~abbEz8$PLF zH`<09=ZYN;90)*U@(xfU{zJqUw;7re;VPufrx}#!o#WU&yrT`5>#t5xb3CFGzD7fng{H74DLH zxndlGV#GfN3{#9ORWY)dVq`JJ$Qr8{rr78PozwPa%^sJn&%em+gY(EThuWR=Votyj zbimCWYMeuQ(Zlp23`$imk^@cWsrz{|Ue(lMh76g~q%#d4F(PsQmc%`(2z>sm*8?h! z^dK=nYT&Sipa;)Sq6hJ;?c(YR&d`G!nI6m>qX$2rZU{X%fd+`wLl4?n`#*O^`@f(lda4}jYBAIxe*)6XGYQD$D0Ztb%#z# zbf}5ti$=$FbMZjs}6hDY(pNd#&km_?O{;vVE+-Twf*-GUB>-6vuQ;PjY@>)`StmW_N9tmGma-UQF7_!qe^2US$N5_Idsc8FZEpdO>=LC?N3}-lhYfr)*=3gV zQqZifLQNe;_aLN9v+S}+_}-Jxl_=$;$a!<5oH=l&q~chYl!@z8IjJ)I99T_a6Eepm z=i&Z|SI*0oa&qN79Gst?St#Y?%Xwu|P8n#bL0TO>=SUj149tpBW)hA|Vxt=i0CzWu z(_k2mCEN2j-ZcM%$G6s^8nbNM$u(ztBz_Hx&h=T4O?DgZ?O|7@+U&oH6we)8`mHD( z(An31F%ZE-tJf;d8WzK*DRjL)yl_y7#Ybuul6wF%z=etmX|0Ra{_5L091Y6Et=65! z#Zr0-Hpk7*`PjeRSV4_|J#9?X=>@)N>Q~r(a!5P1Q>Hj=L?4yKP8IenwS6yP6a6#> z_P)+z^FmUny>OQb!9!FcT{-Gd51BeKh{Y}-*Bs?sTQB-@8rpFBIfC@_>nuWWkUAv} zN*^Rd9xijdiWg!^b#AI&NBDc;u0;k8IyUGc?H}|5;XrL-s_wpQ_rE5zU^2=?wbXs6EkL zS>^Tw_X@_@lWEA4*%Jd^sM8O_Br1pH#(9|i@35^b_**KF#zd>(#!HY^8y|4Jmh4I2s6FY^SrI;~p`G?K z?8$3SQTNcbg_9O(mKhB;kW*sW5#IPJxd-5FUs6i&Q zCtXjovM!yVtQ1{Y51v`p->3}C7k+{q&Y9J){E5*K zTD^(z>~ATr03F7|chE+MYA(~ZQ1jVCKP+wc$bTRTWMI*X3eC(mevae1n)$}r$PuL& zxzT8wgh;odb0OG zhRZV+u?rMD_$5ITj&d#7Wk$Txuo*TaSczLB_3*UQ zJY7n4n!>%54gwJ<{_b zfN-D>f1&+!4!<9-^x*X%f5QkLAZm2o)No)>MTP?(Mx(K5cgASL#{qS`Nj}ZTn>YUb z%<*RSH&iYkZz50dfR(b2T!8?MQcq5hst|GLz>#9)jL8~C8UGF0vpy%GWe1DHUD>uz z{Ijxgpr$KZ`Unybaf`w;YY=PCDI1VC_eGr`soDe{uQybjT$)xmXt@d(0Mc#M)!9xul-0H^Vr_h{hx9+ z_KQ+oyC~&#h9A)6DYVWIRMq$8B+7qi0lyyEXY~ZsK0jAAgf=Tood>>Tg7h{Ait`5<7gl@tl;zEy?ri|=XA`W?AAJ<&}gzZ*y)4w3; zW44>(OWT_^iSp{`K?2lZ%QM<@9}W z`cH7#{CYa&9L>mxjBr4k@jb<9t28%noj4n4Ms^2ldy(5uyq|O1Rb=&U>$^!dO%aWI z9br=jD|)o6LrRv!1G)|m%^48ZlOw??FK=+lkC1Id-Rj-G*RXs3X4rHCOc^fR?1wpY z8u2BN*9T5H22MLpzU#2PEiD%G>m3Fa0@2e-6rLoK4 zH_^tk6h3(XS=(E2&MS(2yQOg}eUqSd4BAKkpu zy4%Jd11APc4#MH0ef(2*UV=^Da8&_Dn7!709`Tbw*5-Gp`CG!epKxja732kLdo@|;EAW2PYP2LyrWLli$oEj)NL zvVpat&-dMi=?%t@V=g*$h{UH!LwZ&V#M?n&?>I<=J1Dpp+2m8VX7Ru(K-|H4_ZwPO z@?X`(AEEe@K%{(2wgn30lKuuoiR7d0N9{kVHdHeke)z#W$AxwK@y64lSt;!wSlZ`-rp6PdySn&EO#LW>=Eu=f z%ko0}$?HOve_MSIUeAfU59|rUA|i2=*o8A83>}{$UV**-!=BbH@PeHVs$95H=sKBA&14s}pmOW%KWN>d0k4 zMfF#dev=O)OJga9e^dl|hZz9spnW+d+*S<&;?;iv;?>%DMUE#~$ndKz$tQOaMz3-+ z5Sd;&`a|QiNBSL^|8QHOnu%QNvF(;>F|5r;kHC4$`f0W`X5pcCOk@{g^%F<9#Yg9k zT=d9fgt~xC0|{i_R)(^ACn;MxMXsu_ zs|8^U1h3+?lU@O70!6Sxg{)z;^9|B4&Lq!t7>8ZbP6y>rgOkg)vpL%(edP+yxzKTP z0Q*cc@x=h<-?Y8aaItg#-HDCoJA$Tk+mVKI=4V=jjbFg-Q9QEa#>zO1%xpWG_VKzn zB=obek(bqnf?T$6)4FWw0JKG>evQxVko!#7(nj)xF~+ThWpRJK95HBb;Wj1-w7Vsd zMulGCft_d?B<6USfG!ZxO$S9-Pe0O$w?o1MPl5<8RT zbavQ|WJQVusqR}njkay(_6a`$V3TRNlSjR>dI})X8mU49#|aoH<}r+F1PWwgy3n};7L<-iH)Yrr{cv2kgDmLv>!JH zo>#zt^5)~zirP7xB^i!!q$G?eNFQhYzjDMh4SkCeAra&3imA?}&+tVh;z*QV zq$Y$%e=Gd;n|X7PC~30c5w;w51mn})#$9OjV|1V75cNE_Zx`A=EE@Z8ZJAow)N2=> zpt+gw1Rd(e2sE8#(0m@vukgh>OWg+Gj{sk{0pA8Z*^gGl-K5dz4lTL^GqdwtttVKj zRh>X9NVTeys5Nw5i$>jo#&mUYV-G330e69ft)U;`p0;RJ|h1XRHZ31BRUC+4dF z43a$tZIAH4Z&)Q`r2(|LbWCk=YW$IQvL$1bJ8}tIh0`dmndtw_K1V))l<2G0dw)u< zfGMiy89yeH8rgi58m#5PsvuIaCIqf{|kB}TZKT3{g5u7mSM8R z^TSI}&CsG|XrT-(l%a((v`~f?mZ1eoRCWR)#XP5n@kHY&ihGgqMm({+gz3LfdCMl= zNX@6OBHl`fw>*uvAAmRHVtI_`V$pK3sJU3wTr89e(;*AxVxe3tEEi(v5-!*2RRUYc1v zBme94QOsdF0t6>j1a}Uc$)~%?7>T}q#RQ_EXUR5>PAwuQp%SB0KVF{pJrDi8>F1FN z9{3#T3t0S0XgG_0G!mEa;yJKb@Jo_nILl6y;VhBBb}L-bitRC+6s0$4PO5ShLU_r5 z>B8MW{tW;ey9EFHF~&P;cJ$5|(ed!EOosF$>3*=U_@PJ0KRtI)O; zLNm21<1IfFrMJb{yd<#3j=^|rM~O8NA1^U0Rg@R@py3I`cC`zDXzidicx(!zVfDzG zkj{Dr-OLH;aLk3qP=}oZF4z*l`%`V&D8+N6gKz35UTY5{_^e6C*)8<#Zc9 z8TkmCwo-1mtDJ_P139)>Ucfei70<)uiTtMSBN;j|f)=jky%TMw2I>y6AA%F{d0MwEBm@>=2uurwO;bTK z0w)vqCkH|%wdFpFIM|w*%Tx1(x}=w*S3z)*pAJppV4{R$WVy-tFoE!Ni$3}(5PZxX zy*i1QlW?gH(PbVf%PK8Rk?2#TbEs5N!iNVDQ?2%Ugfsl-$cxMpzLVrz6c0*mGXXY4 z`$Cd>M8(#pbwm6b0ZtRp0fNsI(~l`zcBFu|uwfhD1P9okz=t$>!ZVH!474+kw03-8 zss8vt)E_KjfBL_7egJIH5rVXlHtevcr)bxC=koaqMd7>(?ag7|3)@}jm~4NJ-EgJP z$G(c)yoQ6=cOef8eGD9KMz~HPKAh1Sy=1O_C51M)v>lgz`c#r~>>(x*I~KazQy>eB zFyH+cuIO+k(;Nh^ew77tSzFakPqV%f&a9?r3e${A)NeG=;j}RGSP=_iPxs-kZdHm;x)7Re(&(0_${$z* z^ZgCltx}J=K|iF!Q58|Qu1;>k2$_R}Dy&4%sq8pllV)hBXO;(-O{TEVU&qyKs2 zr#Apk`&ik=k2G1rqO#qqe?x|D$%Tb?A{qlDuGywEzlOwkq;8R>6q;Ih8!dBGjk-9h0gdXYX5UNGJi;i*uGoQ7VdewN) z2%$qLfy873TdU#Pob64Ci9+z~rezp!S8O(K+Jp?GV`EPLj5f-QD1Tf+GTZ4x8&UCY z<7a46lS9f&h5ZL~rQ_5aJ60Q~yo{D|VE+r_iYtAzQ)mtO*Nv)yi=CnXL6QJ#}L^rUhKubAofMab!@ zysL$J5Mt-AVU5`j7}6P0|13H<%)V5C{aBoLqF-#m9_6KC&qwUnZw@K4jGT>k)Zgh2M)-kc?6dYX#@NrpBF#Q{(HTZFI>ln>GhtLAqMp!+SAuKUG#xIpG z9;6D-=i7Z*;L#&po^AY8l_M5OfrreJ# zw5_*D33z1wqGeshDBdlN;9MO0=w$RPWs0;cK5)vojo72lTs&?!?8YQjaP^LUWZo_P z*^86E(ju&6;-T72*NDL_`M{_@*nne@?EEd+XLb_FH7w=RseT9OZ3mkNsjn2Hr_=d* zZ2j_a(V0Ah_J(m-4dWeN%%hWPK2 zGb@>t48Kxg!8DxC4YC{dVY``2%WhF2t5Dc(rC`{eOaR%=kTUUmU>RUU#1%*!_&lMU zM|y|*k4|*@6G|LHZfQjdD_%K3asPW|J%>)Im4)fRZ%Ca`S{nBJjwwGy2pfASFtynk zc@#WG?TcGpLDw|BO3nrdny!7aYnrJr^UOo#k&@}^Xl)AJ`Yb>`PTQCMJp6)|c{`Z_ zGy~>6)fT0g3Z!35+*$VwNE+vDOGI?{-O+~p08OX<0G)F30UQj)cCVqk1@Y84yQpEQ zJ@BJyNk)8#WdUMwM!=Opr2sN{qk%-i(@(FN(JxYv0>7q_6)!xr6Xw7nKS{s6uQ!|niUlW5pepk-%NSCe)@weXXUN(I$B zurmdf&@76hAaVh6RP743={HcGV-BpcpQ=^u4t-i9xYn%9!&)^$KT|die-YvTi0Z-5 zq1{n;UU_ysm5L8QpuaE!dz27#6nsee7_URG3=QIYz5M&a z(H|}|D%+q4sF4{r1Yj2(w|aa>)lb#pYykTOF{;&=3Af8?nRtNM$8_yfdP`y5<=*Hrt2Mq`Ogv3 zkF;^63lohTmzru-?OOCQMKM#j**b^&#T}*HFgw z<9^w00O=$5Kzo!2T1o3#z%l&sw(7AcKPE2KJ2Wq{+F?p{XtcE*DBL5x z&srP315&*Pn9W)MiFD?p4Lfl`l-qvPAnNi5JfazCrt zd&(myP2X;pgJAs6ND7olj;o)0M!PL$?PeX*ZVr+k+U+akoED0#V*U0#zFZVrRJgYs5F$jB}{tQM^Zh(D>7zTagRTEwFsEwne`DqL)P;RR?(lD zm8+4KK7oq6d^|n-Y$^tY5ywa(ei#Tpvjp8fME!l#eC=o-9gmOifEtsHAdA2_82iUa zYjoY?Eb*84{8cs9%DPo|mDE@bzN%HLsw%Cks)|dj<+WDB%Cf3jYi-@CRYm@C;0fic zsH!ZlTxs>KioUL_H53&WTm4m5UHTsyR63 zRaa4LEvl%f@)h|@tW}j2_gW3s(wd@GCA9{(zpSRDsF;&Meyib%+H0<;om*41 zMt$>D)YbaccYk?lsmgC+nSnkKgq&y_0=3fTr&h0GOr{hidM9=qKb-L zMU^s?SAttBcqC9-iMm^B>MB>2S29q9`i^}py@d?16~>63wCidsYPCjmLKtqpi$$n^ z+IJMPw<9WjrwY|*!@TeRsIGD)gkC3>uPnPg2dr0=WR{gzQ(hp_%Cr8|&XT0rs!NJX zYW+nOIdymWYl?iN6f-Zr&CSG0{Am8wU&v< zTUO<-svKV%jSA^hIS*#p8Pa+?hlUIe+vjt*RfU!x@NmmD9C{aUSg1m~nr<+Og3Okqb_?A2ly()g;@?`98;6Y#bI|$6==`j~12Qas$U#-^}4wHGN2h?JDe4p+SB6RrwmO;^_)i z*ssFWIUL`qLA75Es$sj2$6K#f)2pykg-29qnakr-RcKhq!<|kwor<>?sdNP#ZuM|D zsK#4WyJXOw4mExShwXQ87%Jz`uC<2>|G#g4wf~QA|5R0Og(}=y&E;bl^Bnmh4<6ah z;h-AsS7D-k3^~@dTG+=XQPD{W&1Vh4T%QB zLBPU`F;Bz{|ho z0Q&$_00#jx0EYoZKnvUx`G6^aWq@|Tdcb_Z4S>~v9e^7EHv@J8Du5xt6u5H+0XG09 znNe>*avUju`G6^KE!6{N0JZ~m0QLi#;deGnLwSI9z=MDr0F&Te>;lXOQ~)~wlg>gs zxg7yh;CI~sC<1l@mH~DFb^!JRCSjw;ayHTfrT~@!+5tNN^8vd6s{uu9F?0g<0d@hF zVXyBXU_Iawzz#rjJm>+ffGOBw%>c~6wro9MA7CfoLBK9R1+X7*7*GK;W4qdL4&nh* z0aJd0@&P*lTL6===iUX_1$YoJ1E&eB3E&$r1F#;@3)lfz2G|AI0yqrV0hsbLqz4oM z`w9LW^%lU_;kL@%*XBgcEHVmTM6HV{s34H$9T$xh<_e(1k88=bbv`O zLf--V01Zh9Z%4U+Wg)}^*7rib1a~6-JIHr8;sJ;EA|6l~fLsBSK0rHP1pLRS7huPy z$QQ8dGw=y}a_}hR2RIC<04m3k5A4f^FHwKMVL%0t9*BWG#`|tD80#&@i{q!7+l^B# zgeQ_W?;jbN%}LnbiRlf7WR>0m_u5u+CM8%FB_u5t&R%0~G5qkN`LnN0zKk-U^iKTm z!a^HWMIgblG{Ma9cOV!79wspQ5WOG&iT+B$IulZhRQinw_5%-*5Z;=r;(v|cLExd% z(TDPR0sn~}rX%_g{&oEC10GEmeF*<4{uBL^Q7q@5E@tP$3-R4?cyhB90#9kRhu+hH zF9Uvjy*w)Ycs$(|p#0S)Vp&l>4JtiMRP-S{ofMsj{uz}XrXl(eJ)J})dUVm~L-==r z_X0nj{tMv63Gnf7W!NXkCmHyR3GfSmPn`gt3w+81_zK{yz>nvr75JnH@Q(v;nIIp! zJZT<>C)vIXJibvgGLNL-NU#KA3eHL}d(U=2k0}0A#8dg-r+Db853%mE7ae&%yuGNFPG_3PK=uV~m)Pv_2Np0i;ij z$xpC45|U7%g$ZUVZw@IooOjDq`iaYHny@@4(o`V573qCy`iF?$O|c6TtU=Slgp`fM z@A@gugp9ycXM#5-D=en^4_A9E5+*SqUjH zfsj8KM!g+6J$er5JCGhjS@a=2c^UZ53Gg2RzZv+eqtpfi@o^gX5b%wJM=>fNC^pK6 zH|BOQl6fX0RNey274iExBQ{E}K7`K&J_Y#k{k#HrEAW?U#DGNK3Va6em`X$+YPZLM zw*wy}Y(Dizv~wVKae}o0t+oUGYj=V<7up)TK8DpF^qruui{`?|ybpmlV=n#^!n6Jv z6GYcvM23OVnP{$xc!(qV5Pk;mEx^x=VmZI*Dm|t$(TC_gz_){*>9-m~axMm50iMPg zot#;Fqn`%ivJ&hu{|Y$+6K!zSa3mwiuaUkBYlc?Jli4LE2N;={I;OsAJzfT#6>Ej@ z_4p9@B;d(BY7g;y8h8uvFayzt@QE}RABU&*m_vpV_!Ac=EcQMgg1My9SVD&=zsy$f*JLuKt zyYOtMJ}Bz?r3d*~zvGr3j$)1yjKY3O{1~B6cx{C{t_}qy63bE#6@*tk@ zP`m#ccv=(wln{(hj8m$;WPH+S!H2nM_o4jWMY?{Z8?WDA0N*zbPx>7XRzkqPtd_Hp z?9KXEuHS9YZ?utSiiv&!=nJt1{b%Ssps%2+VojPDmDjlOZwbjuTwIaHb&x)n1 ze6%dWq8~N;J?EY?nR#>HwCKmb_I}#A?|tvN=l(kP+;i`{_Xcz2TZP}A&_ch0K!k4q z15$WKZEWCztwofZ5pN80<~_c zk2K;vh7 z5fq#C$ur>B3x24MvZMa-Tj1@$*Aa@g!3=K+>C5GqXpI49I~nkGFvVKO*_>s5T?oGs z;gx!h;*o23M|>mCJQs(Lu5_KSnSUBPneK+}}4=_jg$2p?|5ik0{;&#B;G0@?v$o1$&0#eY-l|g`jQ3+Q$;CnfxR5%gS-8 z5PlB?jM@^}>tng1a(03KMy#!LP&!;cqdb$KG7!IEgx`(uX$p^0{SBY^7!@*y-!#2P z@j!v5eXa-n%b@QzP#OO?@K&tZbo%fw0zZGd*ndbr0cO;m2Y|PO{@Mz9lBZ{Y-wnJ& z^v;LK?yj^qNjvu$aVjg&_06FV`XK7rRL@>f0@{$JKx@)kFhWKF&xh^yUc`+6}g0Q@V*%Zl+}Plfg`PnaqGg;B%cdRjT( z4tF*%1B}H;1NH)?#^Y-a@ekfLKR;MMej94EYW%ScNL0@c5x>u1o>d-Sd0gT4rSWwb zj9RgT0TD(vjHOEPeua2=VBf<@;TIBBjoizAg`z&yNwZ(HGuHGu+<2Ht^mIh@%zI3{9*WwM~;i% z!%bFx1I<*vlTO3WgP5P9-iw{%hLZ^9pL{7}=LMb|A3zqRQ^ET^t-4qvvI9#%Jp zZx5LVZr_cFe{?_AuF2rEm2F*lePa;DOOnb~Yd;yiK;;|8{A~Au`T5&1p=&(tA$xtX zjV(4=Jp9YX0G`C}Ippu)x3LacmOr)^S$?#fY)>u^0#)`NO+^YJvTHDgWp9GaiUk~-_`V389W>L^UcRt{&;?|X_=d|2_KI)%dwEZjC}4Qjk5c^kh7O_%)wtpgU4; zYUcfQ?#-h3$Zfrj@jcU}lL-M`$quJO0a@UrLB z^Ps^`9l+y9X_4=f%U9s>yHwu)uda2K%SCTF&-kwRj^J7^B(miE{&x*r#)lgXpmd(^ z{rkT6!@l=0nuXj+#&hml;qnuvD(t>~v&;TU^XlyT>Igq^AI6h5rLJ{6DB-UbdNI45 z7k@vWC+OCYOP(XiIx{DjcMK`)j)YwHp_17BR0w%807TQHp`+53*BoX?7^WeUTtVVf zKPCqi$_`zZoqLMQzT$H&=3VnW!l4d*CTFMVW8XZN=jh2UYT!Nc_X%za;*m zGbDZSSBc*ve!uuz#2*!ZTKpTt-z)xI;y)n%ZZe536x{O~G>AB$eiMzTkaGbW{uJKZG8~oCw z@yGNjJAs|c1*Z*rb~0_N*TgO=CAZ7l_KB&sZByA?`l4*wX`{EXWO2-COYfQhKjW7s zO~`a+vPfHRnM__hgiU60Ns1s(gCn0 zgXq;xSta;C!*6iTbbj#I2qS8DJc>74y59m$eD>YsGQH;bBt#@!?H@n?ryd(&dgV`R zX*T}!2QwO++ICql{BKhGin%$CC-|%nr*B-8AHN-TBlADgdD533no~4b-d^ECgL zGM`#G?+1Pk%B%aRdF>yOlE7QB3Cl0fIpA0ps&VfDem44h{1%tP)ZU~$Da2p*`|>vm zQ2cuUpWg)FEie#?kM1Y$mK^+r;JVL!K=AFrDL=Z;9!r|+XzvK&dOiU)%kD1VwsL({ z=tn;8GTkiyI?mI;DL*Y27+m?E4Ex;1=XJm-UEQ~?^lxJNrOvGANfe*m2H?bB_pd9@ zU4rX=bN)t-7_`^+RiGdB^^1EM7p4{7LEx01*4tf1*Fi4H{apC#d4|0Pc4)53`Pu9; z(kb|A!S`Kc@J)j60Iqo!h`notUiWjO$FjRWfc^>Pf0@g4{H+a7jWE zba5=KCjxJk=SuHw0eBv`HUGk&{>&x0AM&r`d|q%@@}uSY&j3Ep3%#Dt`JV7ux>(xf zR+kZ-=g_`lDzBa^xL)XofZOtUT>$<$;iKmzcx@NF_6x4(1yJqme$ISq$M@%jUe8H< zN$6K%;zjwKYjYX7MDWXj+w%Diq4(>bzY|=~-K-Nnvx4h+80|;TDjqZG>U@~?Do}p( zTneqjk?RC*%TK@1>p2oydnK0%pub+|^*o2l!@Yv*Ija?hjq_c_uXGvVwQ)Q>FSwqQ zdAHEN787a8f2;H_)q8Ql^}GYE>5`iiT+fxfXkZ5mMTU>a?QMeZ7kut2ml0lT$J1e@ zf4#x+sgT{Pn}yH$245~XZNH&>>bVoBM|Rf%x7G9Q0r)S3kDeRTesL-ql=RiyZ7%br zzhh0V9k`AEM!{QTzG$c&=NjOYPd$gDdq(IRdaDvu`1b|ye^L2#yG%{?2f1_6MNgr8 ziarXzk0X0src?Q+qsGijh zTZN|5TY`pev5fnt0jG4GkGV|VDv{m@+*Ymu;FKRdFLglZCxnll2Yrp;HwnF-6FemN z{esW^waW;9E5OrJO8-WKZ{bkKX+{>$hJCB{_7#F(3Y_wz=aNnne23tAPUdXE{}#Bd zTz3h*o>$^;LWuK};ChY^-G|*GG!Uh$=b8AMIYQnHoa$50A!<8r5qwtk|4qVY1~}*c z0+*4u3Vt_mTRtBTz(e>}!}*c^3o&4Kq2PKBl=s{qFIDNYz_ckPYN@%5X^pO#M*PrdcDgPsbE<+s>up2nl=k8k! zj#qYj@u1HFv3C!^Qy?cJkB;zri}c?|fD?Z`FSb`A9Tr^A<*9yt)u|?3Jx_;fXV(SX zmTmz!@tGeQB(cQN&RTQ zdn0hN8{@)1BMH0-IOS8%A3iVm)q=NsE+hOc6;GcMJTCdya(#>W<2#E_4?idLdhT*U z_@9oBN%@(*$?)fIN(gHMPWiuYufg%kZU7G~4<9r*f2%`CNpZ=Ku4nHPT+cnaziSekXkP zNk5(DSeTw;!bSP?kMA24kGV|wPNBb7@PnNOKPdS31lRNFYHzfxH2n3PYOl~=CV0E( zPt|)n1lRM*+J8R=oaAjp?2R$ubHC8vcO8yG4F_1n&}D&#Asv z@Cv`O6p}$9PJ^zkuuxmO?>RIeJh!eY2g6nzlwSxBxuIER!KCcm6&yPY~ zv-^bd+3NE9Ap<*K7F^FI^IlhkJukSPU)TD7Jvt`UkJD@TXxRq^*K^H35&k=YD~&+x zeNg#G`)WJhA3*=F0r=C5W4!)|$q#>fMhY(re?2#@_SIQX)Rcd{FG1_g1#Zjd`-NWb z9l*4O-Ccs~d4A>pIB@FkbB~P}Lizo-0RB%0;0;jtw)~t0oaJ+m%beBzyj<{+HiN%Q z(!D|OR%tJYAG>=5*ZTle{tpT6AKx#9z*D+~EiT9BJy&?j3ht~iIPZl8ej9M=kF6gw zxGUrJ7ldBVt-l4dsj{{_jv_M@)}AH8=%`|tk%r+T2xtCjB5a||C>$~7Z_7Xqhz z?ptN(wY+Z=+&>T9CAi<-`v`Db{eN2M^?rnur2EYP`u`Goy~pEPp+D($CZGO(dV%12 zzlP=^0bHMj1$);9;5Q2&y_X>*ksc19p9{c$EqwI8l1~bs*FjLJp7nl;ZoxMR9{1_V zDZ%wVl~qFjgy4E_0Ppoe*lz^adn?3jcV2y-r0e6q0l2MxJ|Of*-{5kJDsKk@=zq*O z6vMp68INsJoDOZdoN|8uvEIH}n7Ncybh^vO8L%=$> zC5zdVm&fPYOmaMy8up@zZZEze)mBJOmfE`h;OMFRM6pzu%oo~vE_S?ujhol3-spM7 zscl}&8&75@nDNG5&&%fddoA5L72EvB=Q!d|%-b4?yam@*Wrp*^NFzRq!ez5gYO+); zO^uGWr5w*&-#6rKj3tIW&q1|iwr8mhz0$ZxcYon}E(_#&>Abf+m*1Amd1>4zR`iln zGjyF;A%{E0(rrDR-QD19A@b;(wE3Rw;T!SfF2{V6PUF_GU5HZg=26vDk$PVL(CST5 zFFLpunKe0=XQ|`O{r!pPus6JV&BiEERIx}`1Gv~~vgCMc-#WN@Q>+h7G8n17U}|k7 zI_=Q6ZO@S74y^Hv>ke#;t?BbR+B(`$&aK0Vz^kvowWHoORFSUjnNnP?*&5oEKx!nj z9qpcn%e3-4y_Gf4?C_j;|$xJa*>dz%h z{n<<|&B+B)_ZFs#V{|!KbSMr<**Y8x;%fry-YFtTUb&2kNhq;h3d$ArQIjwv?ADPO)H8S{AZMHyTOr|4QHjy68|g5rKT! zsr{LXAX6ef)V7fxZQE_hw3kX2OVgGrvEMyLDcWvji*!{~1c?!3-JUr=VHVp@Kw7s6 zda}heQ^j2nmf=``e_t+NL>E90pf|6pt^Jq`SjzixUrikP4%!8rNOb9Pxyk%4e^b@Z zOfT9ZMi=T)BeZ+wALk(8`Z=zf=tO36`>vtn1mdMl<@j%(TiZ^(!P|~25lw56j#)qL z4T(zZ%O#6N=tFNSHo*8>td59c97tt~Meav&>fXV$eq|7+oThW}8UY!!-37R|g9<>} zC|G!EHC&D8aH4}75~c1`X;&eG!7|c`OC5J+)0qiYsL|@|v_&xNQ0Kk`RN0AQX0jBE zZJKdR?G2IoT7gw0& z7E3B-JUf%gMYCniY8T>T7PgbcWojF;tJK_ypl%KDz4U^%SH$bvs${A3p^v8=o{%lZ2aq|YgME1)Q~B>RcIoe z33Mw>AN3R<%aI<`XNjYP9&p#r!VpTPQAxQP?g?MFCx$gdmaalV zJb)W?ODuG;VY&xbZRg=l+>eHNf8`%IDLI}zplRf)R&81ok!1^dt%!?fC%ma*Ce7pc z0O>y1fUj&QhZ2li3J*~uYXCxJX)n6k5Nef~Ab@ow{mN;js||?6EA}x&&n%R*y@c3?lUhunzQ2?(Un=Pq-5g zZyAhW@{pRDNp8zdMdSj%24iKOkt48T! zVG?HEsIT3QH5hkRX-`$2LGlryqLa=pj^#0rIkpiS>7X9J&@C_#WJ-0R7-sC_30AYC@c;VlI;@ z_^U5b>cTaiIy;qVe1JV$+j7S?Bt+ZORJ~rslk6U+v{{G0x*fyJ_Q(hs8_=d>cwGT?RFL%_Ha^*sDojmgyy@&@X)1{_R;z1U59wxevXAfc=b?bB z+`9Hh6GU6BQsEM%9&?Rie$OXi#wIU=REC7GKdTN#aub))h8 zw8bc}HeDyorwaSI?!2CC=%LLqnv))f4Eg3vW@HPf>x&EdLjPn2Ya!I&9**CX7n+VC zsnx@$WUJ0-M}rQFNWOa%%j{zoo49V*sLmlEb%-LJWJD}XvPc$>XUN8^79SRj8noye-)xY^-~ip%K|S&b^}B)(u15;%?epy zV1?WI4>v~;r*8BYU2a{_Awd?-!s;Z|nWc7V^|ngO8fmbR(ozZ6%fc!z7J6x|m!z)% z?NlRub)`O1Glu_(#Ea#{R?%3qKwFC~7Es|jA+gnny;TEd0YOWx7Uf$f&5lmA!f+mc zcf{Ban@1}=4xdR>bD zMohId(F~dvJ^VBBNC!mYt3I2 zN|YwE1?rmBHAHqdeK{GF+OnhAy7GrxlvXYK_cAIUYN{a=R~I33K)W_%+Sgl$wRDmqD;~d zq|I#~Ufb!z{`Q#}uYhCmdHhAFY-yJ_-R_$vLGmnp#E4!U1J5@npuAb@Tdl>K;z+>6 z&e~6eHsSM`VwX?s#3DXEIpVLdkt*x!!*_oEzG@YcY9*9T8?(_SPMY7qnr>fMi800k zJ^(`3)7LL-k)efqdxB%FOUFBrkUKzu7NWLRT z;}3^i4jUGrmERVA8W$e>cxc9djejKM;tltg32@UQes##BvTJfC-p@310zf#u zkIB3l9{NoF+wq_7Leltp|H#{MzKz@-A&V* zVCNlppgPs~b1g0h&B?~f5q_(!KP~_D2&edKgih~0apmBG-Y3XY$q4-2w>Ar)P zW%y5hh@9S|vG1$q-Ish>eoUz2{1^`;>Kb3~p?Of^`}6C6{~6*^o;1GRS2KH@@qdZ9 zB)b}4@3q-4@!O?gRlYTSm7|w^>FfPBFMkaIa7u697Ty*A9Re!zU%t0w*2FmCt8k5{ z{vw1|#y|4N2uD2fh>7o)7au!^?xi7KT7EqT{G7y}^F?SiY|Fo=Ae>0r%ZPc<_*;Fr zvQwYbOgZfe+Wtq)%b7|jVZK*RTl-!`bZ6u-^EM^}Adb~IJ@~KnrgBN$)0Y1Jr%n79 LRvBM^{ptT72MY~X literal 0 HcmV?d00001 diff --git a/slsDetectorCalibration/jungfrauExecutables/jungfrauClusterFinder.cpp b/slsDetectorCalibration/jungfrauExecutables/jungfrauClusterFinder.cpp new file mode 100644 index 000000000..f7b2f916a --- /dev/null +++ b/slsDetectorCalibration/jungfrauExecutables/jungfrauClusterFinder.cpp @@ -0,0 +1,161 @@ +//#include "sls/ansi.h" +#include + + +//#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