minor modifications for interpolation and mythe data structure

This commit is contained in:
2019-02-06 16:22:17 +01:00
parent 9e5ec6a57b
commit dcad6c80ce
17 changed files with 582 additions and 158 deletions

View File

@@ -0,0 +1,45 @@
class Stat
{
public:
Stat() : n(0), m(0.), m2(0.) {}
void Clear()
{
n = 0;
m=0;
m2=0;
}
void Push(double x)
{
m+=x;
m2+=x*x;
n++;
}
int NumDataValues() const
{
return n;
}
double Mean() const
{
return (n > 0) ? m/n : 0.0;
}
double Variance() const
{
return ( (n >0 ) ? (m2/n-m*m/(n*n)) : 0.0 );
}
double StandardDeviation() const
{
return sqrt( Variance() );
}
private:
int n;
double m, m2;
};

View File

@@ -0,0 +1,116 @@
#ifndef COMMONMODESUBTRACTION_H
#define COMMONMODESUBTRACTION_H
#include <cmath>
class commonModeSubtraction {
/** @short class to calculate the common mode of the pedestals based on an approximated moving average*/
public:
/** constructor
\param nn number of samples for the moving average to calculate the average common mode
\param iroi number of regions on which one can calculate the common mode separately. Defaults to 1 i.e. whole detector
*/
commonModeSubtraction(int iroi=1, int ns=3) : nROI(iroi), nsigma(ns) {
mean=new double[nROI];
mean2=new double[nROI];
nCm=new double[nROI];
};
/** destructor - deletes the moving average(s) and the sum of pedestals calculator(s) */
virtual ~commonModeSubtraction() {delete [] mean; delete [] mean2; delete [] nCm;};
/** clears the moving average and the sum of pedestals calculation - virtual func*/
virtual void Clear(){
for (int i=0; i<nROI; i++) {
mean[i]=0;
nCm[i]=0;
mean2[i]=0;
}};
/** adds the average of pedestals to the moving average and reinitializes the calculation of the sum of pedestals for all ROIs. - virtual func*/
virtual void newFrame(){
for (int i=0; i<nROI; i++) {
// if (nCm[i]>0) cmStat[i].Calc(cmPed[i]/nCm[i]);
nCm[i]=0;
mean[i]=0;
mean2[i]=0;
}};
/** adds the pixel to the sum of pedestals -- virtual func must be overloaded to define the regions of interest
\param val value to add
\param ix pixel x coordinate
\param iy pixel y coordinate
*/
virtual void addToCommonMode(double val, int ix=0, int iy=0) {
int iroi=getROI(ix,iy);
// if (iroi==0) val=100;
// else val=-100;
// if (isc>=0 && isc<nROI) {
if (iroi>=0 && iroi<nROI) {
mean[iroi]+=val;
mean2[iroi]+=val*val;
nCm[iroi]++;
}
};
/** gets the common mode i.e. the difference between the current average sum of pedestals mode and the average pedestal
\param ix pixel x coordinate
\param iy pixel y coordinate
\return the difference between the current average sum of pedestals and the average pedestal
*/
virtual double getCommonMode(int ix=0, int iy=0) {
int iroi=getROI(ix,iy);
/* if (iroi==0) */
/* return 100; */
/* else */
/* return -100; */
if (iroi>=0 && iroi<nROI) {
if (nCm[iroi]>0)
return mean[iroi]/nCm[iroi];
}
return 0;
};
/** gets the common mode i.e. the difference between the current average sum of pedestals mode and the average pedestal
\param ix pixel x coordinate
\param iy pixel y coordinate
\return the difference between the current average sum of pedestals and the average pedestal
*/
virtual double getCommonModeRMS(int ix=0, int iy=0) {
int iroi=getROI(ix,iy);
if (iroi>=0 && iroi<nROI) {
if (nCm[iroi]>0)
return sqrt(mean2[iroi]/nCm[iroi]-(mean[iroi]/nCm[iroi])*(mean[iroi]/nCm[iroi]));
}
return 0;
};
/**
gets the common mode ROI for pixel ix, iy -should be overloaded!
*/
virtual int getROI(int ix, int iy){ (void) ix; (void) iy; return 0;};
protected:
double *mean; /**<array of moving average of the pedestal average per region of interest */
double *mean2; /**< array storing the sum of pedestals per region of interest */
double *nCm; /**< array storing the number of pixels currently contributing to the pedestals */
int nsigma; /** number of rms above which the pedestal should be considered as a photon */
const int nROI; /**< constant parameter for number of regions on which the common mode should be calculated separately e.g. supercolumns */
};
#endif

View File

@@ -74,7 +74,7 @@ class mythen3_01_jctbData : public slsDetectorData<short unsigned int> {
for (ib=0; ib<nb; ib++) {
if (word&(1<<bit[ib])) {
cout << "+" ;
val[iw+nch*(ib/nb)]|=(1<<idr);
val[iw+nch/nb*(ib)]|=(1<<idr);
} else {
cout << "-" ;
}
@@ -96,8 +96,8 @@ class mythen3_01_jctbData : public slsDetectorData<short unsigned int> {
ii++;
}//end for
cout << "Decoded "<<ii << " samples"<< endl;
cout << "Should be "<< nch/nb*dr+off << " samples"<< endl;
cout << "M3.01 Decoded "<<ii << " samples"<< endl;
cout << "M3.01 Should be "<< nch/nb*dr+off << " samples"<< endl;
return val;
}

View File

@@ -57,8 +57,9 @@ class mythen3_02_jctbData : public mythen3_01_jctbData {
int ioff=0;
int idr=0;
int ib=0;
int iw=0;
int ich=0;
int ii=0;
int iw=0;
bit[0]=17;//19;
bit[1]=6;//8;
idr=0;
@@ -68,7 +69,8 @@ class mythen3_02_jctbData : public mythen3_01_jctbData {
wp=(int64_t*)ptr;
for (iw=0; iw<nch/nb; iw) {
word=*wp;;
word=*wp;
if (ioff<off) {
ioff++;
cout <<"*";
@@ -78,10 +80,11 @@ class mythen3_02_jctbData : public mythen3_01_jctbData {
for (ib=0; ib<nb; ib++) {
if (word&(1<<bit[ib])) {
cout << "+" ;
val[iw+nch*(ib/nb)]|=(1<<idr);
val[iw+nch/nb*(ib)]|=(1<<idr);
} else {
cout << "-" ;
}
}
// cout << iw+nch/nb*(ib)<< " " ;
}//end for()
}
@@ -94,14 +97,14 @@ class mythen3_02_jctbData : public mythen3_01_jctbData {
cout <<dec << iw<<endl;
iw++;
}//end if()
}//end else()
wp+=1;
ii++;
}//end for
cout << "Decoded "<<ii << " samples"<< endl;
cout << "Should be "<< nch/nb*dr+off << " samples"<< endl;
cout << "M3.02 Decoded "<<ii << " samples"<< endl;
cout << "M3.02 Should be "<< nch/nb*dr+off << " samples"<< endl;
return val;
}

View File

@@ -125,7 +125,11 @@ class etaInterpolationAdaptiveBins : public etaInterpolationPosXY {
virtual void prepareInterpolation(int &ok, int nint=1000)
virtual void prepareInterpolation(int &ok) {
prepareInterpolation(ok, 1000);
}
virtual void prepareInterpolation(int &ok, int nint)
{
ok=1;

View File

@@ -48,7 +48,7 @@ int main(int argc, char *argv[]) {
}
int p=10000;
int fifosize=1000;
int nthreads=1;
int nthreads=24;
int nsubpix=25;
int etabins=nsubpix*10;
double etamin=-1, etamax=2;

View File

@@ -92,16 +92,23 @@ int main(int argc, char *argv[]) {
#endif
int nSubPixels=nsubpix;
//eta2InterpolationPosXY *interp=new eta2InterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax);
eta2InterpolationCleverAdaptiveBins *interp=new eta2InterpolationCleverAdaptiveBins(NC, NR, nsubpix, etabins, etamin, etamax);
#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);
interp->prepareInterpolation(ok);//, MAX_ITERATIONS);
#endif
// return 0;
#endif
#ifdef FF
@@ -146,8 +153,9 @@ int main(int argc, char *argv[]) {
if (nframes==0) f0=lastframe;
nframes++;
}
quad=interp->calcQuad(cl.get_cluster(), sum, totquad, sDum);
if (sum>cmin && totquad/sum>0.8 && totquad/sum<1.2 && sum<cmax ) {
//quad=interp->calcQuad(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 && sum<cmax ) {
nph++;
// if (sum>200 && sum<580) {
// interp->getInterpolatedPosition(cl.x,cl.y, totquad,quad,cl.get_cluster(),int_x, int_y);
@@ -155,16 +163,32 @@ int main(int argc, char *argv[]) {
// 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, 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;
interp->addToImage(int_x, int_y);
// 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);
// 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
// }
@@ -179,12 +203,12 @@ int main(int argc, char *argv[]) {
interp->writeFlatField(outfname);
#endif
}
}
}
}
}
}
fclose(f);
#ifdef FF
interp->writeFlatField(outfname);
@@ -203,7 +227,7 @@ int main(int argc, char *argv[]) {
}
}
}
cout << "Read " << nframes << " frames (first frame: " << f0 << " last frame: " << lastframe << " delta:" << lastframe-f0 << ")"<<endl;
cout << "Read " << nframes << " frames (first frame: " << f0 << " last frame: " << lastframe << " delta:" << lastframe-f0 << ") nph="<< nph <<endl;
interp->clearInterpolatedImage();
#endif

View File

@@ -0,0 +1,269 @@
#include "ansi.h"
#include <iostream>
//#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 "etaInterpolationCleverAdaptiveBins.h"
//#include "etaInterpolationRandomBins.h"
using namespace std;
#define NC 400
#define NR 400
#define MAX_ITERATIONS (nSubPixels*100)
#define MAX_EBINS 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 n_ebins=1;
if (argc>iarg)
n_ebins=atoi(argv[iarg++]);
//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;
int iebin=0;
double eb_size=(cmax-cmin)/n_ebins;
#ifndef NOINTERPOLATION
// eta2InterpolationPosXY *interp[MAX_EBINS];
eta2InterpolationCleverAdaptiveBins *interp[MAX_EBINS];
for (int i=0; i< n_ebins; i++) {
//interp[i]=new eta2InterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax);
interp[i]=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;
for (int i=0; i< n_ebins; i++) {
sprintf(fname,argv[2],i);
interp[i]->readFlatField(fname);
interp[i]->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; ix<NC; ix++) {
for (iy=0; iy<NR; iy++) {
for (isx=0; isx<nsubpix; isx++) {
for (isy=0; isy<nsubpix; isy++) {
totimg[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)]=0;
}
}
}
}
#ifdef FF
sprintf(outfname,argv[2]);
#endif
int irun;
for (irun=runmin; irun<runmax; irun++) {
sprintf(infname,argv[1],irun);
#ifndef FF
sprintf(outfname,argv[3],irun);
#endif
f=fopen(infname,"r");
if (f) {
cout << infname << endl;
nframes=0;
f0=-1;
while (cl.read(f)) {
totph++;
if (lastframe!=cl.iframe) {
lastframe=cl.iframe;
// cout << cl.iframe << endl;
// f0=cl.iframe;
if (nframes==0) f0=lastframe;
nframes++;
}
//quad=interp->calcQuad(cl.get_cluster(), sum, totquad, sDum);
quad=interp[0]->calcEta(cl.get_cluster(), etax, etay, sum, totquad, sDum);
if (sum>cmin && totquad/sum>0.8 && totquad/sum<1.2 && sum<cmax ) {
nph++;
iebin=(sum-cmin)/eb_size;
if (iebin>=0 && iebin<n_ebins) {
// if (sum>200 && 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[iebin]->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[iebin]->addToImage(int_x, int_y);
#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[iebin]->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
for (int i=0; i<n_ebins; i++) {
sprintf(outfname,argv[3],i);
interp[i]->writeInterpolatedImage(outfname);
}
#endif
#ifdef FF
for (int i=0; i<n_ebins; i++) {
sprintf(outfname,argv[2],i);
cout << outfname << " " << argv[2] << " " << i << endl;
interp[i]->writeFlatField(outfname);
}
#endif
}
}
}
}
fclose(f);
#ifdef FF
for (int i=0; i<n_ebins; i++) {
sprintf(outfname,argv[2],i);
cout << outfname << " " << argv[2] << " " << i << endl;
interp[i]->writeFlatField(outfname);
}
#endif
#ifndef FF
for (int i=0; i<n_ebins; i++) {
sprintf(outfname,argv[3],i,irun);
interp[i]->writeInterpolatedImage(outfname);
img=interp[i]->getInterpolatedImage();
for (ix=0; ix<NC; ix++) {
for (iy=0; iy<NR; iy++) {
for (isx=0; isx<nsubpix; isx++) {
for (isy=0; isy<nsubpix; isy++) {
totimg[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)]+=img[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)];
}
}
}
}
//interp[i]->clearInterpolatedImage();
}
cout << "Read " << nframes << " frames (first frame: " << f0 << " last frame: " << lastframe << " delta:" << lastframe-f0 << ")"<<endl;
#endif
} else
cout << "could not open file " << infname << endl;
}
#ifndef FF
sprintf(outfname,argv[3], 11111);
WriteToTiff(totimg, outfname,NC*nsubpix,NR*nsubpix);
#endif
#ifdef FF
interp[iebin]->writeFlatField(outfname);
#endif
cout << "Filled " << nph << " (/"<< totph <<") " << endl;
return 0;
}

View File

@@ -69,7 +69,6 @@ int main(int argc, char *argv[]) {
time_t begin,end,finished;
if (argc > 4) {
socketip2 = argv[3];
portnum2 = atoi(argv[4]);
@@ -113,6 +112,7 @@ int main(int argc, char *argv[]) {
int multisize=size;
int dataSize=size;
char dummybuff[size];
@@ -685,6 +685,10 @@ int main(int argc, char *argv[]) {
mt->pushData(buff);
mt->nextThread();
mt->popFree(buff);
} else {
cprintf(RED, "Incomplete frame: received only %d packet\n", packetNumber);
length = zmqsocket->ReceiveData(0, dummybuff, size);
}

View File

@@ -562,10 +562,10 @@ int *getClusters(char *data, int *ph=NULL) {
static void writeClusters(FILE *f, single_photon_hit *cl, int nph, int fn=0){
#ifndef OLDFORMAT
if (fwrite((void*)&fn, 1, sizeof(int), f))
if (fwrite((void*)&nph, 1, sizeof(int), f))
#endif
/* #ifndef OLDFORMAT */
/* if (fwrite((void*)&fn, 1, sizeof(int), f)) */
/* if (fwrite((void*)&nph, 1, sizeof(int), f)) */
/* #endif */
for (int i=0; i<nph; i++) (cl+i)->write(f);
};
void writeClusters(FILE *f, int fn=0){

View File

@@ -39,14 +39,16 @@ class single_photon_hit {
//fwrite((void*)this, 1, 3*sizeof(int)+4*sizeof(double)+sizeof(quad), myFile);
// if (fwrite((void*)this, 1, sizeof(int)+2*sizeof(int16_t), myFile))
#ifdef OLDFORMAT
if (fwrite((void*)&iframe, 1, sizeof(int), myFile))
#endif
//#ifdef OLDFORMAT
if (fwrite((void*)&iframe, 1, sizeof(int), myFile)) {};
//#endif
#ifndef WRITE_QUAD
//printf("no quad ");
if (fwrite((void*)&x, 2, sizeof(int16_t), myFile))
return fwrite((void*)data, 1, dx*dy*sizeof(int), myFile);
#endif
#ifdef WRITE_QUAD
// printf("quad ");
int qq[4];
switch(quad) {
case TOP_LEFT:
@@ -102,18 +104,43 @@ class single_photon_hit {
size_t read(FILE *myFile) {
//fread((void*)this, 1, 3*sizeof(int)+4*sizeof(double)+sizeof(quad), myFile);
#ifdef OLDFORMAT
if (fread((void*)&iframe, 1, sizeof(int), myFile))
#endif
//#ifdef OLDFORMAT
if (fread((void*)&iframe, 1, sizeof(int), myFile)) {}
//#endif
#ifndef WRITE_QUAD
// printf( "no quad \n");
if (fread((void*)&x, 2, sizeof(int16_t), myFile))
return fread((void*)data, 1, dx*dy*sizeof(int), myFile);
#endif
#ifdef WRITE_QUAD
int qq[4];
// printf( "quad \n");
if (fread((void*)&x, 2, sizeof(int16_t), myFile))
if (fread((void*)qq, 1, 4*sizeof(int), myFile)) {
quad=TOP_RIGHT;
int mm=qq[0];
for (int i=1; i<4; i++) {
if (qq[i]<mm) {
switch (i) {
case 1:
quad=TOP_LEFT;
break;
case 2:
quad=BOTTOM_RIGHT;
break;
case 3:
quad=BOTTOM_LEFT;
break;
default:
;
}
}
}
switch(quad) {
case TOP_LEFT:
data[0]=0;