fixed slsDeetctorcalibration - including interpolation for rectangular pixels

This commit is contained in:
bergamaschi 2020-03-19 15:25:00 +01:00
parent 0e3cd00579
commit 78dd96357d
15 changed files with 673 additions and 479 deletions

View File

@ -1,37 +1,37 @@
#initialchecks 0 initialchecks 0
############################################# #############################################
### edit with hostname or IP address of your detector ### edit with hostname or IP address of your detector
############################################ ############################################
hostname bchip181+ hostname bchip181+
#############################################
### edit with 10 Gbs IP of your server
############################################
udp_dstip 10.1.1.102
#############################################
### edit with any number in the subnet of your server (first 3 numbers as above)
############################################
udp_srcip 10.1.1.19
udp_dstport 33410
#############################################
### edit with 10 Gbs IP of your server
############################################
#rx_zmqip 10.1.1.102
#rx_zmqport 30001
#############################################
### edit with 1 Gbs IP of PC where you will run the GUI
############################################
#zmqip 129.129.202.136
#zmqport 40001
############################################# #############################################
### edit with hostname or 1Gbs IP address of your server ### edit with hostname or 1Gbs IP address of your server
############################################ ############################################
rx_hostname mpc2011 rx_hostname mpc2608
rx_tcpport 1954
#############################################
### edit with 10 Gbs IP of your server
############################################
udp_dstip 10.1.2.117
#############################################
### edit with any number in the subnet of your server (first 3 numbers as above)
############################################
udp_srcip 10.1.2.19
udp_dstport 32411
#############################################
### edit with 10 Gbs IP of your server
############################################
rx_zmqip 10.1.2.117
rx_zmqport 50003
#############################################
### edit with 1 Gbs IP of PC where you will run the GUI
############################################
zmqip 129.129.202.86
zmqport 50001
tengiga 1 tengiga 1

View File

@ -199,7 +199,7 @@ template <class dataType> class analogDetector {
\param nns reference to number of subpixels for interpolating detector, will always be 1 in this case \param nns reference to number of subpixels for interpolating detector, will always be 1 in this case
\returns number of pixels of the detector image \returns number of pixels of the detector image
*/ */
virtual int getImageSize(int &nnx, int &nny, int &nns){nnx=nx; nny=ny; nns=1; return nx*ny;}; virtual int getImageSize(int &nnx, int &nny, int &nnsx, int &nnsy){nnx=nx; nny=ny; nnsx=1; nnsy=1; return nx*ny;};
/** /**
Returns data size of the detector image matrix Returns data size of the detector image matrix
\param nnx reference to pixel size in x \param nnx reference to pixel size in x
@ -423,7 +423,7 @@ template <class dataType> class analogDetector {
if (g==0) g=-1.; if (g==0) g=-1.;
} }
return stat[iy][ix].getPedestalRMS();//divide by gain? return stat[iy][ix].getPedestalRMS()/g;//divide by gain?
} }
return -1; return -1;
}; };
@ -441,8 +441,9 @@ template <class dataType> class analogDetector {
\returns pedestal value \returns pedestal value
*/ */
virtual double* getPedestal(double *ped){ virtual double* getPedestal(double *ped){
if (ped==NULL) if (ped==NULL) {
ped=new double[nx*ny]; ped=new double[nx*ny];
}
for (int iy=0; iy<ny; iy++) { for (int iy=0; iy<ny; iy++) {
for (int ix=0; ix<nx; ix++) { for (int ix=0; ix<nx; ix++) {
ped[iy*nx+ix]=stat[iy][ix].getPedestal(); ped[iy*nx+ix]=stat[iy][ix].getPedestal();
@ -459,8 +460,9 @@ template <class dataType> class analogDetector {
\returns pedestal rms \returns pedestal rms
*/ */
virtual double* getPedestalRMS(double *ped=NULL){ virtual double* getPedestalRMS(double *ped=NULL){
if (ped==NULL) if (ped==NULL) {
ped=new double[nx*ny]; ped=new double[nx*ny];
}
for (int iy=0; iy<ny; iy++) { for (int iy=0; iy<ny; iy++) {
for (int ix=0; ix<nx; ix++) { for (int ix=0; ix<nx; ix++) {
ped[iy*nx+ix]=stat[iy][ix].getPedestalRMS(); ped[iy*nx+ix]=stat[iy][ix].getPedestalRMS();
@ -675,7 +677,7 @@ template <class dataType> class analogDetector {
delete [] gm; delete [] gm;
return 1; return 1;
} }
return NULL; return 0;
} }
/** /**
@ -700,7 +702,7 @@ template <class dataType> class analogDetector {
delete [] gm; delete [] gm;
return 1; return 1;
} }
return NULL; return 0;
} }

View File

@ -39,9 +39,9 @@ class slsDetectorData {
*/ */
slsDetectorData(int npx, int npy, int dsize, int **dMap=NULL, dataType **dMask=NULL, int **dROI=NULL): nx(npx), ny(npy), dataSize(dsize) { slsDetectorData(int npx, int npy, int dsize, int **dMap=NULL, dataType **dMask=NULL, int **dROI=NULL): nx(npx), ny(npy), dataSize(dsize) {
int el=dsize/sizeof(dataType);
xmap=new int[dsize/sizeof(dataType)]; xmap=new int[el];
ymap=new int[dsize/sizeof(dataType)]; ymap=new int[el];
// if (dataMask==NULL) { // if (dataMask==NULL) {
@ -65,7 +65,8 @@ class slsDetectorData {
dataROIMask[i][j]=1; dataROIMask[i][j]=1;
} }
// } // }
for (int ip=0; ip<dsize/sizeof(dataType); ip++){
for (int ip=0; ip<el; ip++){
xmap[ip]=-1; xmap[ip]=-1;
ymap[ip]=-1; ymap[ip]=-1;
} }
@ -205,7 +206,7 @@ class slsDetectorData {
virtual void getPixel(int ip, int &x, int &y) {x=xmap[ip]; y=ymap[ip];}; virtual void getPixel(int ip, int &x, int &y) {x=xmap[ip]; y=ymap[ip];};
virtual dataType **getData(char *ptr, int dsize=-1) { virtual dataType **getData(char *ptr, int dsize=-1) {
int el=dsize/sizeof(dataType);
dataType **data; dataType **data;
int ix,iy; int ix,iy;
data=new dataType*[ny]; data=new dataType*[ny];
@ -213,7 +214,7 @@ class slsDetectorData {
data[i]=new dataType[nx]; data[i]=new dataType[nx];
} }
if (dsize<=0 || dsize>dataSize) dsize=dataSize; if (dsize<=0 || dsize>dataSize) dsize=dataSize;
for (int ip=0; ip<(dsize/sizeof(dataType)); ip++) { for (int ip=0; ip<(el); ip++) {
getPixel(ip,ix,iy); getPixel(ip,ix,iy);
if (ix>=0 && ix<nx && iy>=0 && iy<ny) { if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
data[iy][ix]=getChannel(ptr,ix,iy); data[iy][ix]=getChannel(ptr,ix,iy);
@ -231,8 +232,9 @@ class slsDetectorData {
for(int i = 0; i < ny; i++) { for(int i = 0; i < ny; i++) {
data[i]=new double[nx]; data[i]=new double[nx];
} }
int el=dsize/sizeof(dataType);
if (dsize<=0 || dsize>dataSize) dsize=dataSize; if (dsize<=0 || dsize>dataSize) dsize=dataSize;
for (int ip=0; ip<(dsize/sizeof(dataType)); ip++) { for (int ip=0; ip<el; ip++) {
getPixel(ip,ix,iy); getPixel(ip,ix,iy);
if (ix>=0 && ix<nx && iy>=0 && iy<ny) { if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
data[iy][ix]=getValue(ptr,ix,iy); data[iy][ix]=getValue(ptr,ix,iy);

View File

@ -13,7 +13,7 @@ template <class dataType> class ghostSummation {
/** constructor /** constructor
\param xt crosstalk \param xt crosstalk
*/ */
ghostSummation(slsDetectorData<dataType> *d, double xt) : xtalk(xt),det(d) { ghostSummation(slsDetectorData<dataType> *d, double xt) : xtalk(xt),det(d), nx(1), ny(1) {
if (det) if (det)
det->getDetectorSize(nx,ny); det->getDetectorSize(nx,ny);
ghost=new double[nx*ny]; ghost=new double[nx*ny];
@ -22,6 +22,9 @@ template <class dataType> class ghostSummation {
ghostSummation(ghostSummation *orig) { ghostSummation(ghostSummation *orig) {
xtalk=orig->xtalk; xtalk=orig->xtalk;
det=orig->det; det=orig->det;
nx=1;
ny=1;
det->getDetectorSize(nx,ny);
ghost=new double[nx*ny]; ghost=new double[nx*ny];
} }

View File

@ -107,11 +107,11 @@ class interpolatingDetector : public singlePhotonDetector {
singlePhotonDetector::clearImage(); singlePhotonDetector::clearImage();
}; };
int getImageSize(int &nnx, int &nny, int &ns) { int getImageSize(int &nnx, int &nny, int &nsx, int &nsy) {
if (interp) if (interp)
return interp->getImageSize(nnx, nny, ns); return interp->getImageSize(nnx, nny, nsx, nsy);
else else
return analogDetector<uint16_t>::getImageSize(nnx, nny, ns); return analogDetector<uint16_t>::getImageSize(nnx, nny, nsx, nsy);
}; };
@ -251,6 +251,7 @@ int addFrame(char *data, int *ph=NULL, int ff=0) {
} }
virtual int getNSubPixels(){ if (interp) return interp->getNSubPixels(); else return 1;} virtual int getNSubPixels(){ if (interp) return interp->getNSubPixels(); else return 1;}
virtual int setNSubPixels(int ns) { virtual int setNSubPixels(int ns) {
if (interp) { if (interp) {
pthread_mutex_lock(fi); pthread_mutex_lock(fi);

View File

@ -13,7 +13,7 @@
class eta2InterpolationBase : public virtual etaInterpolationBase { class eta2InterpolationBase : public virtual etaInterpolationBase {
public: public:
eta2InterpolationBase(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny, ns, nb, emin, emax) { eta2InterpolationBase(int nx=400, int ny=400, int ns=25, int nsy=25, int nb=-1, int nby=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny, ns, nsy, nb, nby, emin, emax) {
/* if (etamin>=etamax) { */ /* if (etamin>=etamax) { */
/* etamin=-1; */ /* etamin=-1; */
@ -37,7 +37,7 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
int corner; int corner;
corner=calcQuad(data, tot, totquad, sDum); corner=calcQuad(data, tot, totquad, sDum);
if (nSubPixels>2) if (nSubPixelsX>2 || nSubPixelsY>2)
calcEta(totquad, sDum, etax, etay); calcEta(totquad, sDum, etax, etay);
getInterpolatedPosition(x,y,etax,etay,corner,int_x,int_y); getInterpolatedPosition(x,y,etax,etay,corner,int_x,int_y);
@ -53,7 +53,7 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
int corner; int corner;
corner=calcQuad(data, tot, totquad, sDum); corner=calcQuad(data, tot, totquad, sDum);
if (nSubPixels>2) if (nSubPixelsX>2 || nSubPixelsY>2 )
calcEta(totquad, sDum, etax, etay); calcEta(totquad, sDum, etax, etay);
getInterpolatedPosition(x,y,etax,etay,corner,int_x,int_y); getInterpolatedPosition(x,y,etax,etay,corner,int_x,int_y);
@ -93,7 +93,7 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
; ;
} }
double etax=0, etay=0; double etax=0, etay=0;
if (nSubPixels>2) { if (nSubPixelsX>2 || nSubPixelsY>2) {
cc[0][0]=cl[xoff+3*yoff]; cc[0][0]=cl[xoff+3*yoff];
cc[1][0]=cl[xoff+3*(yoff+1)]; cc[1][0]=cl[xoff+3*(yoff+1)];
cc[0][1]=cl[xoff+1+3*yoff]; cc[0][1]=cl[xoff+1+3*yoff];
@ -133,7 +133,7 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
; ;
} }
double etax=0, etay=0; double etax=0, etay=0;
if (nSubPixels>2) { if (nSubPixelsX>2 || nSubPixelsY>2) {
cc[0][0]=cl[xoff+3*yoff]; cc[0][0]=cl[xoff+3*yoff];
cc[1][0]=cl[xoff+3*(yoff+1)]; cc[1][0]=cl[xoff+3*(yoff+1)];
cc[0][1]=cl[xoff+1+3*yoff]; cc[0][1]=cl[xoff+1+3*yoff];
@ -182,31 +182,31 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
} }
if (nSubPixels>2) { if (nSubPixelsX>2 || nSubPixelsY>2 ) {
ex=(etax-etamin)/etastep; ex=(etax-etamin)/etastepX;
ey=(etay-etamin)/etastep; ey=(etay-etamin)/etastepY;
if (ex<0) { if (ex<0) {
cout << "x*"<< ex << endl; cout << "x*"<< ex << endl;
ex=0; ex=0;
} }
if (ex>=nbeta) { if (ex>=nbetaX) {
cout << "x?"<< ex << endl; cout << "x?"<< ex << endl;
ex=nbeta-1; ex=nbetaX-1;
} }
if (ey<0) { if (ey<0) {
cout << "y*"<< ey << endl; cout << "y*"<< ey << " " << nbetaY << endl;
ey=0; ey=0;
} }
if (ey>=nbeta) { if (ey>=nbetaY) {
cout << "y?"<< ey << endl; cout << "y?"<< ey << " " << nbetaY << endl;
ey=nbeta-1; ey=nbetaY-1;
} }
xpos_eta=(((double)hhx[(ey*nbeta+ex)]))+dX ;///((double)nSubPixels); xpos_eta=(((double)hhx[(ey*nbetaX+ex)]))+dX ;///((double)nSubPixels);
ypos_eta=(((double)hhy[(ey*nbeta+ex)]))+dY ;///((double)nSubPixels); ypos_eta=(((double)hhy[(ey*nbetaX+ex)]))+dY ;///((double)nSubPixels);
} else { } else {
xpos_eta=0.5*dX+0.25; xpos_eta=0.5*dX+0.25;
@ -347,10 +347,10 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
#endif #endif
#ifndef MYROOT1 #ifndef MYROOT1
int ex,ey; int ex,ey;
ex=(etax-etamin)/etastep; ex=(etax-etamin)/etastepX;
ey=(etay-etamin)/etastep; ey=(etay-etamin)/etastepY;
if (ey<nbeta && ex<nbeta && ex>=0 && ey>=0) if (ey<nbetaY && ex<nbetaX && ex>=0 && ey>=0)
heta[ey*nbeta+ex]++; heta[ey*nbetaX+ex]++;
#endif #endif
return 0; return 0;
}; };
@ -360,22 +360,22 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
// cout << "ff" << endl; // cout << "ff" << endl;
calcDiff(1, hhx, hhy); //get flat calcDiff(1, hhx, hhy); //get flat
double avg=0; double avg=0;
for (ipx=0; ipx<nSubPixels; ipx++) for (ipx=0; ipx<nSubPixelsX; ipx++)
for (ipy=0; ipy<nSubPixels; ipy++) for (ipy=0; ipy<nSubPixelsY; ipy++)
avg+=flat[ipx+ipy*nSubPixels]; avg+=flat[ipx+ipy*nSubPixelsX];
avg/=nSubPixels*nSubPixels; avg/=nSubPixelsY*nSubPixelsX;
for (int ibx=0 ; ibx<nSubPixels*nPixelsX; ibx++) { for (int ibx=0 ; ibx<nSubPixelsX*nPixelsX; ibx++) {
ipx=ibx%nSubPixels-nSubPixels/2; ipx=ibx%nSubPixelsX-nSubPixelsX/2;
if (ipx<0) ipx=nSubPixels+ipx; if (ipx<0) ipx=nSubPixelsX+ipx;
for (int iby=0 ; iby<nSubPixels*nPixelsY; iby++) { for (int iby=0 ; iby<nSubPixelsY*nPixelsY; iby++) {
ipy=iby%nSubPixels-nSubPixels/2; ipy=iby%nSubPixelsY-nSubPixelsY/2;
if (ipy<0) ipy=nSubPixels+ipy; if (ipy<0) ipy=nSubPixelsY+ipy;
// cout << ipx << " " << ipy << " " << ibx << " " << iby << endl; // cout << ipx << " " << ipy << " " << ibx << " " << iby << endl;
if (flat[ipx+ipy*nSubPixels]>0) if (flat[ipx+ipy*nSubPixelsX]>0)
hintcorr[ibx+iby*nSubPixels*nPixelsX]=hint[ibx+iby*nSubPixels*nPixelsX]*(avg/flat[ipx+ipy*nSubPixels]); hintcorr[ibx+iby*nSubPixelsX*nPixelsX]=hint[ibx+iby*nSubPixelsX*nPixelsX]*(avg/flat[ipx+ipy*nSubPixelsX]);
else else
hintcorr[ibx+iby*nSubPixels*nPixelsX]=hint[ibx+iby*nSubPixels*nPixelsX]; hintcorr[ibx+iby*nSubPixelsX*nPixelsX]=hint[ibx+iby*nSubPixelsX*nPixelsX];
} }

View File

@ -13,35 +13,12 @@
class eta3InterpolationBase : public virtual etaInterpolationBase { class eta3InterpolationBase : public virtual etaInterpolationBase {
public: public:
eta3InterpolationBase(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationBase(nx, ny, ns, nb, emin, emax) { eta3InterpolationBase(int nx=400, int ny=400, int ns=25, int nsy=25, int nb=-1, int nby=-1, double emin=1, double emax=0) : etaInterpolationBase(nx, ny, ns, nsy, nb, nby, emin, emax) {
// cout << "e3ib " << nb << " " << emin << " " << emax << endl; // cout << "e3ib " << nb << " " << emin << " " << emax << endl;
/* if (nbeta<=0) { */ /* if (nbeta<=0) { */
/* nbeta=nSubPixels*10; */ /* nbeta=nSubPixels*10; */
/* } */ /* } */
if (etamin>=etamax) {
etamin=-1;
etamax=1;
}
etastep=(etamax-etamin)/nbeta;
#ifdef MYROOT1
delete heta;
delete hhx;
delete hhy;
heta=new TH2D("heta","heta",nbeta,etamin,etamax,nbeta,etamin,etamax);
hhx=new TH2D("hhx","hhx",nbeta,etamin,etamax,nbeta,etamin,etamax);
hhy=new TH2D("hhy","hhy",nbeta,etamin,etamax,nbeta,etamin,etamax);
#endif
#ifndef MYROOT1
/* delete [] heta; */
/* delete [] hhx; */
/* delete [] hhy; */
/* heta=new int[nbeta*nbeta]; */
/* hhx=new float[nbeta*nbeta]; */
/* hhy=new float[nbeta*nbeta]; */
#endif
// cout << nbeta << " " << etamin << " " << etamax << endl; // cout << nbeta << " " << etamin << " " << etamax << endl;
}; };
@ -88,7 +65,7 @@ class eta3InterpolationBase : public virtual etaInterpolationBase {
double etax, etay; double etax, etay;
if (nSubPixels>2) { if (nSubPixelsX>2 || nSubPixelsY>2 ) {
calcEta3(cl,etax,etay, totquad); calcEta3(cl,etax,etay, totquad);
} }
return getInterpolatedPosition(x,y,etax, etay,quad,int_x,int_y); return getInterpolatedPosition(x,y,etax, etay,quad,int_x,int_y);
@ -101,7 +78,7 @@ class eta3InterpolationBase : public virtual etaInterpolationBase {
double etax, etay; double etax, etay;
if (nSubPixels>2) { if (nSubPixelsX>2 || nSubPixelsY>2 ) {
calcEta3(cl,etax,etay, totquad); calcEta3(cl,etax,etay, totquad);
} }
return getInterpolatedPosition(x,y,etax, etay,quad,int_x,int_y); return getInterpolatedPosition(x,y,etax, etay,quad,int_x,int_y);
@ -117,38 +94,38 @@ class eta3InterpolationBase : public virtual etaInterpolationBase {
double xpos_eta=0,ypos_eta=0; double xpos_eta=0,ypos_eta=0;
int ex,ey; int ex,ey;
if (nSubPixels>2) { if (nSubPixelsX>2 || nSubPixelsY>2 ) {
#ifdef MYROOT1 #ifdef MYROOT1
xpos_eta=(hhx->GetBinContent(hhx->GetXaxis()->FindBin(etax),hhy->GetYaxis()->FindBin(etay)))/((double)nSubPixels); xpos_eta=(hhx->GetBinContent(hhx->GetXaxis()->FindBin(etax),hhy->GetYaxis()->FindBin(etay)))/((double)nSubPixelsX);
ypos_eta=(hhy->GetBinContent(hhx->GetXaxis()->FindBin(etax),hhy->GetYaxis()->FindBin(etay)))/((double)nSubPixels); ypos_eta=(hhy->GetBinContent(hhx->GetXaxis()->FindBin(etax),hhy->GetYaxis()->FindBin(etay)))/((double)nSubPixelsY);
#endif #endif
#ifndef MYROOT1 #ifndef MYROOT1
ex=(etax-etamin)/etastep; ex=(etax-etamin)/etastepX;
ey=(etay-etamin)/etastep; ey=(etay-etamin)/etastepY;
if (ex<0) { if (ex<0) {
/* cout << etax << " " << etamin << " "; */ /* cout << etax << " " << etamin << " "; */
/* cout << "3x*"<< ex << endl; */ /* cout << "3x*"<< ex << endl; */
ex=0; ex=0;
} }
if (ex>=nbeta) { if (ex>=nbetaX) {
/* cout << etax << " " << etamin << " "; */ /* cout << etax << " " << etamin << " "; */
/* cout << "3x?"<< ex << endl; */ /* cout << "3x?"<< ex << endl; */
ex=nbeta-1; ex=nbetaX-1;
} }
if (ey<0) { if (ey<0) {
/* cout << etay << " " << etamin << " "; */ /* cout << etay << " " << etamin << " "; */
/* cout << "3y*"<< ey << endl; */ /* cout << "3y*"<< ey << endl; */
ey=0; ey=0;
} }
if (ey>=nbeta) { if (ey>=nbetaY) {
/* cout << etay << " " << etamin << " "; */ /* cout << etay << " " << etamin << " "; */
/* cout << "3y?"<< ey << endl; */ /* cout << "3y?"<< ey << endl; */
ey=nbeta-1; ey=nbetaY-1;
} }
xpos_eta=(((double)hhx[(ey*nbeta+ex)]));///((double)nSubPixels); xpos_eta=(((double)hhx[(ey*nbetaX+ex)]));///((double)nSubPixels);
ypos_eta=(((double)hhy[(ey*nbeta+ex)]));///((double)nSubPixels); ypos_eta=(((double)hhy[(ey*nbetaX+ex)]));///((double)nSubPixels);
#endif #endif
@ -265,10 +242,10 @@ class eta3InterpolationBase : public virtual etaInterpolationBase {
#endif #endif
#ifndef MYROOT1 #ifndef MYROOT1
int ex,ey; int ex,ey;
ex=(etax-etamin)/etastep; ex=(etax-etamin)/etastepX;
ey=(etay-etamin)/etastep; ey=(etay-etamin)/etastepY;
if (ey<nbeta && ex<nbeta && ex>=0 && ey>=0) if (ey<nbetaY && ex<nbetaX && ex>=0 && ey>=0)
heta[ey*nbeta+ex]++; heta[ey*nbetaX+ex]++;
#endif #endif
return 0; return 0;
}; };

View File

@ -7,7 +7,7 @@
#include <TH2D.h> #include <TH2D.h>
#include <TH2F.h> #include <TH2F.h>
#endif #endif
#include <cmath>
#include "slsInterpolation.h" #include "slsInterpolation.h"
#include "tiffIO.h" #include "tiffIO.h"
@ -15,44 +15,51 @@ class etaInterpolationBase : public slsInterpolation {
public: public:
etaInterpolationBase(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : slsInterpolation(nx,ny,ns), hhx(NULL), hhy(NULL), heta(NULL), nbeta(nb), etamin(emin), etamax(emax) { etaInterpolationBase(int nx=400, int ny=400, int ns=25, int nsy=25, int nb=-1, int nby=-1, double emin=1, double emax=0) : slsInterpolation(nx,ny,ns,nsy), hhx(NULL), hhy(NULL), heta(NULL), nbetaX(nb), nbetaY(nby), etamin(emin), etamax(emax) {
// cout << "eb " << nb << " " << emin << " " << emax << endl; // cout << "eb " << nb << " " << emin << " " << emax << endl;
// cout << nb << " " << etamin << " " << etamax << endl; // cout << nb << " " << etamin << " " << etamax << endl;
if (nbeta<=0) { if (nbetaX<=0) {
//cout << "aaa:" <<endl; //cout << "aaa:" <<endl;
nbeta=nSubPixels*10; nbetaX=nSubPixelsX*10;
}
if (nbetaY<=0) {
//cout << "aaa:" <<endl;
nbetaY=nSubPixelsY*10;
} }
if (etamin>=etamax) { if (etamin>=etamax) {
etamin=-1; etamin=-1;
etamax=2; etamax=2;
} }
etastep=(etamax-etamin)/nbeta; etastepX=(etamax-etamin)/nbetaX;
heta=new int[nbeta*nbeta]; etastepY=(etamax-etamin)/nbetaY;
hhx=new float[nbeta*nbeta]; heta=new int[nbetaX*nbetaY];
hhy=new float[nbeta*nbeta]; hhx=new float[nbetaX*nbetaY];
hhy=new float[nbetaX*nbetaY];
rangeMin=etamin; rangeMin=etamin;
rangeMax=etamax; rangeMax=etamax;
flat= new double[nSubPixels*nSubPixels]; flat= new double[nSubPixelsX*nSubPixelsY];
hintcorr=new int [nSubPixels*nSubPixels*nPixelsX*nPixelsY]; hintcorr=new int [nSubPixelsX*nSubPixelsY*nPixelsX*nPixelsY];
}; };
etaInterpolationBase(etaInterpolationBase *orig): slsInterpolation(orig){ etaInterpolationBase(etaInterpolationBase *orig): slsInterpolation(orig){
nbeta=orig->nbeta; nbetaX=orig->nbetaX;
nbetaY=orig->nbetaY;
etamin=orig->etamin; etamin=orig->etamin;
etamax=orig->etamax; etamax=orig->etamax;
rangeMin=orig->rangeMin; rangeMin=orig->rangeMin;
rangeMax=orig->rangeMax; rangeMax=orig->rangeMax;
etastep=(etamax-etamin)/nbeta; etastepX=(etamax-etamin)/nbetaX;
heta=new int[nbeta*nbeta]; etastepY=(etamax-etamin)/nbetaY;
memcpy(heta,orig->heta,nbeta*nbeta*sizeof(int)); heta=new int[nbetaX*nbetaY];
hhx=new float[nbeta*nbeta]; memcpy(heta,orig->heta,nbetaX*nbetaY*sizeof(int));
memcpy(hhx,orig->hhx,nbeta*nbeta*sizeof(float)); hhx=new float[nbetaX*nbetaY];
hhy=new float[nbeta*nbeta]; memcpy(hhx,orig->hhx,nbetaX*nbetaY*sizeof(float));
memcpy(hhy,orig->hhy,nbeta*nbeta*sizeof(float)); hhy=new float[nbetaX*nbetaY];
hintcorr=new int [nSubPixels*nSubPixels*nPixelsX*nPixelsY]; memcpy(hhy,orig->hhy,nbetaX*nbetaY*sizeof(float));
hintcorr=new int [nSubPixelsX*nSubPixelsY*nPixelsX*nPixelsY];
}; };
@ -61,7 +68,7 @@ class etaInterpolationBase : public slsInterpolation {
virtual void resetFlatField() { virtual void resetFlatField() {
for (int ibx=0; ibx<nbeta*nbeta; ibx++) { for (int ibx=0; ibx<nbetaX*nbetaY; ibx++) {
heta[ibx]=0; heta[ibx]=0;
hhx[ibx]=0; hhx[ibx]=0;
hhy[ibx]=0; hhy[ibx]=0;
@ -70,13 +77,16 @@ class etaInterpolationBase : public slsInterpolation {
}; };
int *setEta(int *h, int nb=-1, double emin=1, double emax=0) int *setEta(int *h, int nb=-1, int nby=-1, double emin=1, double emax=0)
{ {
if (h) { if (h) {
if (heta) delete [] heta; if (heta) delete [] heta;
heta=h; heta=h;
nbeta=nb; nbetaX=nb;
if (nb<=0) nbeta=nSubPixels*10; nbetaY=nby;
if (nbetaX<=0) nbetaX=nSubPixelsX*10;
if (nbetaY<=0) nbetaY=nSubPixelsY*10;
etamin=emin; etamin=emin;
etamax=emax; etamax=emax;
if (etamin>=etamax) { if (etamin>=etamax) {
@ -85,22 +95,24 @@ class etaInterpolationBase : public slsInterpolation {
} }
rangeMin=etamin; rangeMin=etamin;
rangeMax=etamax; rangeMax=etamax;
etastep=(etamax-etamin)/nbeta; etastepX=(etamax-etamin)/nbetaX;
etastepY=(etamax-etamin)/nbetaY;
} }
return heta; return heta;
}; };
int *setFlatField(int *h, int nb=-1, double emin=1, double emax=0) int *setFlatField(int *h, int nb=-1, int nby=-1, double emin=1, double emax=0)
{ {
return setEta(h, nb, emin, emax); return setEta(h, nb, nby, emin, emax);
}; };
int *getFlatField(){return setEta(NULL);}; int *getFlatField(){return setEta(NULL);};
int *getFlatField(int &nb, double &emin, double &emax){ int *getFlatField(int &nb, int &nby, double &emin, double &emax){
nb=nbeta; nb=nbetaX;
nby=nbetaY;
emin=etamin; emin=etamin;
emax=etamax; emax=etamax;
return getFlatField(); return getFlatField();
@ -109,13 +121,13 @@ class etaInterpolationBase : public slsInterpolation {
void *writeFlatField(const char * imgname) { void *writeFlatField(const char * imgname) {
float *gm=NULL; float *gm=NULL;
gm=new float[nbeta*nbeta]; gm=new float[nbetaX*nbetaY];
for (int ix=0; ix<nbeta; ix++) { for (int ix=0; ix<nbetaX; ix++) {
for (int iy=0; iy<nbeta; iy++) { for (int iy=0; iy<nbetaY; iy++) {
gm[iy*nbeta+ix]=heta[iy*nbeta+ix]; gm[iy*nbetaX+ix]=heta[iy*nbetaX+ix];
} }
} }
WriteToTiff(gm, imgname, nbeta, nbeta); WriteToTiff(gm, imgname, nbetaX, nbetaY);
delete [] gm; delete [] gm;
return NULL; return NULL;
}; };
@ -129,16 +141,18 @@ class etaInterpolationBase : public slsInterpolation {
etamax=2; etamax=2;
} }
etastep=(etamax-etamin)/nbeta; etastepX=(etamax-etamin)/nbetaX;
etastepY=(etamax-etamin)/nbetaY;
uint32 nnx; uint32 nnx;
uint32 nny; uint32 nny;
float *gm=ReadFromTiff(imgname, nnx, nny); float *gm=ReadFromTiff(imgname, nnx, nny);
if (nnx!=nny) { /* if (nnx!=nny) { */
cout << "different number of bins in x " << nnx << " and y " << nny<< " !"<< endl; /* cout << "different number of bins in x " << nnx << " and y " << nny<< " !"<< endl; */
cout << "Aborting read"<< endl; /* cout << "Aborting read"<< endl; */
return 0; /* return 0; */
} /* } */
nbeta=nnx; nbetaX=nnx;
nbetaY=nny;
if (gm) { if (gm) {
if (heta) { if (heta) {
delete [] heta; delete [] heta;
@ -146,13 +160,13 @@ class etaInterpolationBase : public slsInterpolation {
delete [] hhy; delete [] hhy;
} }
heta=new int[nbeta*nbeta]; heta=new int[nbetaX*nbetaY];
hhx=new float[nbeta*nbeta]; hhx=new float[nbetaX*nbetaY];
hhy=new float[nbeta*nbeta]; hhy=new float[nbetaX*nbetaY];
for (int ix=0; ix<nbeta; ix++) { for (int ix=0; ix<nbetaX; ix++) {
for (int iy=0; iy<nbeta; iy++) { for (int iy=0; iy<nbetaY; iy++) {
heta[iy*nbeta+ix]=gm[iy*nbeta+ix]; heta[iy*nbetaX+ix]=gm[iy*nbetaX+ix];
} }
} }
delete [] gm; delete [] gm;
@ -178,10 +192,10 @@ float *gethhx()
}; };
virtual int addToFlatField(double etax, double etay){ virtual int addToFlatField(double etax, double etay){
int ex,ey; int ex,ey;
ex=(etax-etamin)/etastep; ex=(etax-etamin)/etastepX;
ey=(etay-etamin)/etastep; ey=(etay-etamin)/etastepY;
if (ey<nbeta && ex<nbeta && ex>=0 && ey>=0) if (ey<nbetaY && ex<nbetaX && ex>=0 && ey>=0)
heta[ey*nbeta+ex]++; heta[ey*nbetaX+ex]++;
return 0; return 0;
}; };
@ -195,80 +209,80 @@ float *gethhx()
float tot_eta=0; float tot_eta=0;
float *etah=new float[nbeta*nbeta]; float *etah=new float[nbetaX*nbetaY];
int etabins=nbeta; // int etabins=nbeta;
int ibb=0; int ibb=0;
for (int ii=0; ii<etabins*etabins; ii++) { for (int ii=0; ii<nbetaX*nbetaY; ii++) {
etah[ii]=heta[ii]; etah[ii]=heta[ii];
tot_eta+=heta[ii]; tot_eta+=heta[ii];
} }
sprintf(tit,"/scratch/eta_%d.tiff",ind); sprintf(tit,"/scratch/eta_%d.tiff",ind);
WriteToTiff(etah, tit, etabins, etabins); WriteToTiff(etah, tit, nbetaX, nbetaY);
for (int ii=0; ii<etabins*etabins; ii++) { for (int ii=0; ii<nbetaX*nbetaY; ii++) {
ibb=(hhx[ii]*nSubPixels); ibb=(hhx[ii]*nSubPixelsX);
etah[ii]=ibb; etah[ii]=ibb;
} }
sprintf(tit,"/scratch/eta_hhx_%d.tiff",ind); sprintf(tit,"/scratch/eta_hhx_%d.tiff",ind);
WriteToTiff(etah, tit, etabins, etabins); WriteToTiff(etah, tit, nbetaX, nbetaY);
for (int ii=0; ii<etabins*etabins; ii++) { for (int ii=0; ii<nbetaX*nbetaY; ii++) {
ibb=hhy[ii]*nSubPixels; ibb=hhy[ii]*nSubPixelsY;
etah[ii]=ibb; etah[ii]=ibb;
} }
sprintf(tit,"/scratch/eta_hhy_%d.tiff",ind); sprintf(tit,"/scratch/eta_hhy_%d.tiff",ind);
WriteToTiff(etah, tit, etabins, etabins); WriteToTiff(etah, tit, nbetaX, nbetaY);
float *ftest=new float[nSubPixels*nSubPixels]; float *ftest=new float[nSubPixelsX*nSubPixelsY];
for (int ib=0; ib<nSubPixels*nSubPixels; ib++) ftest[ib]=0; for (int ib=0; ib<nSubPixelsX*nSubPixelsY; ib++) ftest[ib]=0;
//int ibx=0, iby=0; //int ibx=0, iby=0;
for (int ii=0; ii<nbeta*nbeta; ii++) { for (int ii=0; ii<nbetaX*nbetaY; ii++) {
ibx=nSubPixels*hhx[ii]; ibx=nSubPixelsX*hhx[ii];
iby=nSubPixels*hhy[ii]; iby=nSubPixelsY*hhy[ii];
if (ibx<0) ibx=0; if (ibx<0) ibx=0;
if (iby<0) iby=0; if (iby<0) iby=0;
if (ibx>=nSubPixels) ibx=nSubPixels-1; if (ibx>=nSubPixelsX) ibx=nSubPixelsX-1;
if (iby>=nSubPixels) iby=nSubPixels-1; if (iby>=nSubPixelsY) iby=nSubPixelsY-1;
if (ibx>=0 && ibx<nSubPixels && iby>=0 && iby<nSubPixels) { if (ibx>=0 && ibx<nSubPixelsX && iby>=0 && iby<nSubPixelsY) {
// //
// if (ibx>0 && iby>0) cout << ibx << " " << iby << " " << ii << endl; // if (ibx>0 && iby>0) cout << ibx << " " << iby << " " << ii << endl;
ftest[ibx+iby*nSubPixels]+=heta[ii]; ftest[ibx+iby*nSubPixelsX]+=heta[ii];
} else } else
cout << "Bad interpolation "<< ii << " " << ibx << " " << iby<< endl; cout << "Bad interpolation "<< ii << " " << ibx << " " << iby<< endl;
} }
sprintf(tit,"/scratch/ftest_%d.tiff",ind); sprintf(tit,"/scratch/ftest_%d.tiff",ind);
WriteToTiff(ftest, tit, nSubPixels, nSubPixels); WriteToTiff(ftest, tit, nSubPixelsX, nSubPixelsY);
//int ibx=0, iby=0; //int ibx=0, iby=0;
tot_eta/=nSubPixels*nSubPixels; tot_eta/=nSubPixelsX*nSubPixelsY;
int nbad=0; int nbad=0;
for (int ii=0; ii<etabins*etabins; ii++) { for (int ii=0; ii<nbetaX*nbetaY; ii++) {
ibx=nSubPixels*hhx[ii]; ibx=nSubPixelsX*hhx[ii];
iby=nSubPixels*hhy[ii]; iby=nSubPixelsY*hhy[ii];
if (ftest[ibx+iby*nSubPixels]<tot_eta*0.5) { if (ftest[ibx+iby*nSubPixelsX]<tot_eta*0.5) {
etah[ii]=1; etah[ii]=1;
nbad++; nbad++;
} else if(ftest[ibx+iby*nSubPixels]>tot_eta*2.){ } else if(ftest[ibx+iby*nSubPixelsX]>tot_eta*2.){
etah[ii]=2; etah[ii]=2;
nbad++; nbad++;
} else } else
etah[ii]=0; etah[ii]=0;
} }
sprintf(tit,"/scratch/eta_bad_%d.tiff",ind); sprintf(tit,"/scratch/eta_bad_%d.tiff",ind);
WriteToTiff(etah, tit, etabins, etabins); WriteToTiff(etah, tit, nbetaX, nbetaY);
// cout << "Index: " << ind << "\t Bad bins: "<< nbad << endl; // cout << "Index: " << ind << "\t Bad bins: "<< nbad << endl;
//int ibx=0, iby=0; //int ibx=0, iby=0;
@ -286,46 +300,44 @@ float *gethhx()
double diff=0, d; double diff=0, d;
//double bsize=1./nSubPixels; //double bsize=1./nSubPixels;
int nbad=0; int nbad=0;
double p_tot_x[nSubPixels], p_tot_y[nSubPixels], p_tot[nSubPixels*nSubPixels]; double p_tot_x[nSubPixelsX], p_tot_y[nSubPixelsY], p_tot[nSubPixelsX*nSubPixelsY];
double maxdiff=0, mindiff=avg*nSubPixels*nSubPixels; double maxdiff=0, mindiff=avg*nSubPixelsX*nSubPixelsY;
int ipx, ipy; int ipx, ipy;
for (ipy=0; ipy<nSubPixels; ipy++) { for (ipy=0; ipy<nSubPixelsY; ipy++) {
for (ipx=0; ipx<nSubPixels; ipx++) { for (ipx=0; ipx<nSubPixelsX; ipx++) {
p_tot[ipx+ipy*nSubPixels]=0; p_tot[ipx+ipy*nSubPixelsX]=0;
} }
p_tot_y[ipy]=0; p_tot_y[ipy]=0;
p_tot_x[ipy]=0; p_tot_x[ipy]=0;
} }
for (int ibx=0; ibx<nbeta; ibx++) { for (int ibx=0; ibx<nbetaX; ibx++) {
for (int iby=0; iby<nbeta; iby++) { for (int iby=0; iby<nbetaY; iby++) {
ipx=hx[ibx+iby*nbeta]*nSubPixels; ipx=hx[ibx+iby*nbetaX]*nSubPixelsX;
if (ipx<0) ipx=0; if (ipx<0) ipx=0;
if (ipx>=nSubPixels) ipx=nSubPixels-1; if (ipx>=nSubPixelsX) ipx=nSubPixelsX-1;
ipy=hy[ibx+iby*nbeta]*nSubPixels; ipy=hy[ibx+iby*nbetaX]*nSubPixelsY;
if (ipy<0) ipy=0; if (ipy<0) ipy=0;
if (ipy>=nSubPixels) ipy=nSubPixels-1; if (ipy>=nSubPixelsY) ipy=nSubPixelsY-1;
p_tot[ipx+ipy*nSubPixels]+=heta[ibx+iby*nbeta];
p_tot_y[ipy]+=heta[ibx+iby*nbeta];
p_tot_x[ipx]+=heta[ibx+iby*nbeta];
p_tot[ipx+ipy*nSubPixelsX]+=heta[ibx+iby*nbetaX];
p_tot_y[ipy]+=heta[ibx+iby*nbetaX];
p_tot_x[ipx]+=heta[ibx+iby*nbetaX];
} }
} }
// cout << endl << endl; // cout << endl << endl;
for (ipy=0; ipy<nSubPixels; ipy++) { for (ipy=0; ipy<nSubPixelsY; ipy++) {
cout.width(5); cout.width(5);
//flat_y[ipy]=p_tot_y[ipy];//avg/nSubPixels; //flat_y[ipy]=p_tot_y[ipy];//avg/nSubPixels;
for (ipx=0; ipx<nSubPixels; ipx++) { for (ipx=0; ipx<nSubPixelsX; ipx++) {
// flat_x[ipx]=p_tot_x[ipx];///avg/nSubPixels; // flat_x[ipx]=p_tot_x[ipx];///avg/nSubPixels;
flat[ipx+nSubPixels*ipy]=p_tot[ipx+nSubPixels*ipy];///avg; flat[ipx+nSubPixelsX*ipy]=p_tot[ipx+nSubPixelsX*ipy];///avg;
d=p_tot[ipx+nSubPixels*ipy]-avg; d=p_tot[ipx+nSubPixelsX*ipy]-avg;
if (d<0) d*=-1.; if (d<0) d*=-1.;
if (d>5*sqrt(avg) ) if (d>5*sqrt(avg) )
nbad++; nbad++;
@ -354,8 +366,8 @@ float *gethhx()
float *hhx; float *hhx;
float *hhy; float *hhy;
int *heta; int *heta;
int nbeta; int nbetaX, nbetaY;
double etamin, etamax, etastep; double etamin, etamax, etastepX, etastepY;
double rangeMin, rangeMax; double rangeMin, rangeMax;

View File

@ -9,7 +9,7 @@
class etaInterpolationPosXY : public virtual etaInterpolationBase{ class etaInterpolationPosXY : public virtual etaInterpolationBase{
public: public:
etaInterpolationPosXY(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny,ns, nb, emin,emax){ etaInterpolationPosXY(int nx=400, int ny=400, int ns=25, int nsy=25, int nb=-1, int nby=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny, ns, nsy, nb, nby, emin, emax){
// cout << "epxy " << nb << " " << emin << " " << emax << endl; cout << nbeta << " " << etamin << " " << etamax << endl; // cout << "epxy " << nb << " " << emin << " " << emax << endl; cout << nbeta << " " << etamin << " " << etamax << endl;
}; };
@ -24,14 +24,6 @@ class etaInterpolationPosXY : public virtual etaInterpolationBase{
virtual void prepareInterpolation(int &ok) virtual void prepareInterpolation(int &ok)
{ {
ok=1; ok=1;
#ifdef MYROOT1
if (hhx) delete hhx;
if (hhy) delete hhy;
hhx=new TH2D("hhx","hhx",heta->GetNbinsX(),heta->GetXaxis()->GetXmin(),heta->GetXaxis()->GetXmax(), heta->GetNbinsY(),heta->GetYaxis()->GetXmin(),heta->GetYaxis()->GetXmax());
hhy=new TH2D("hhy","hhy",heta->GetNbinsX(),heta->GetXaxis()->GetXmin(),heta->GetXaxis()->GetXmax(), heta->GetNbinsY(),heta->GetYaxis()->GetXmin(),heta->GetYaxis()->GetXmax());
#endif
///*Eta Distribution Rebinning*/// ///*Eta Distribution Rebinning*///
@ -40,106 +32,128 @@ class etaInterpolationPosXY : public virtual etaInterpolationBase{
double tot_eta=0; double tot_eta=0;
double tot_eta_x=0; double tot_eta_x=0;
double tot_eta_y=0; double tot_eta_y=0;
for (int ip=0; ip<nbeta*nbeta; ip++) for (int ip=0; ip<nbetaX*nbetaY; ip++)
tot_eta+=heta[ip]; tot_eta+=heta[ip];
cout << "total eta entries is :"<< tot_eta << endl; cout << "total eta entries is :"<< tot_eta << endl;
if (tot_eta<=0) {ok=0; return;}; if (tot_eta<=0) {ok=0; return;};
double hx[nbeta]; //profile x double hx[nbetaX]; //profile x
double hy[nbeta]; //profile y double hy[nbetaY]; //profile y
double hix[nbeta]; //integral of projection x double hix[nbetaX]; //integral of projection x
double hiy[nbeta]; //integral of projection y double hiy[nbetaY]; //integral of projection y
// int ii=0; // int ii=0;
double etax;//, etay; double etax, etay;
for (int ib=0; ib<nbeta; ib++) { for (int ib=0; ib<nbetaX; ib++) {
tot_eta_x=0;
tot_eta_y=0;
for (int iby=0; iby<nbeta; iby++) { //tot_eta_y=0;
etax=etamin+iby*etastep;
for (int iby=0; iby<nbetaY; iby++) {
etay=etamin+iby*etastepY;
//cout << etax << endl; //cout << etax << endl;
if (etax>=0 && etax<=1)
hx[iby]=heta[iby+ib*nbeta];
else {
hx[iby]=0;
}
// tot_eta_x+=hx[iby]; // tot_eta_x+=hx[iby];
if (etax>=0 && etax<=1) if (etay>=0 && etay<=1)
hy[iby]=heta[ib+iby*nbeta]; hy[iby]=heta[ib+iby*nbetaX];
else else
hy[iby]=0; hy[iby]=0;
// tot_eta_y+=hy[iby]; // tot_eta_y+=hy[iby];
} }
hix[0]=hx[0];
hiy[0]=hy[0]; hiy[0]=hy[0];
for (int iby=1; iby<nbeta; iby++) { for (int iby=1; iby<nbetaY; iby++) {
hix[iby]=hix[iby-1]+hx[iby];
hiy[iby]=hiy[iby-1]+hy[iby]; hiy[iby]=hiy[iby-1]+hy[iby];
} }
// ii=0; tot_eta_y=hiy[nbetaY-1]+1;
tot_eta_x=hix[nbeta-1]+1;
tot_eta_y=hiy[nbeta-1]+1;
for (int ibx=0; ibx<nbeta; ibx++) { for (int iby=0; iby<nbetaY; iby++) {
if (tot_eta_x<=0) {
hhx[ibx+ib*nbeta]=-1;
//ii=(ibx)/nbeta;
} else //if (hix[ibx]>(ii+1)*tot_eta_x*bsize)
{
//if (hix[ibx]>tot_eta_x*(ii+1)/nSubPixels) ii++;
hhx[ibx+ib*nbeta]=hix[ibx]/tot_eta_x;
}
}
/* if (ii!=(nSubPixels-1)) */
/* cout << ib << " x " << tot_eta_x << " " << (ii+1)*tot_eta_x*bsize << " " << ii << " " << hix[nbeta-1]<< endl; */
//ii=0;
for (int ibx=0; ibx<nbeta; ibx++) {
if (tot_eta_y<=0) { if (tot_eta_y<=0) {
hhy[ib+ibx*nbeta]=-1; hhy[ib+iby*nbetaX]=-1;
//ii=(ibx*nSubPixels)/nbeta; //ii=(ibx*nSubPixels)/nbeta;
} else { } else {
//if (hiy[ibx]>tot_eta_y*(ii+1)/nSubPixels) ii++; //if (hiy[ibx]>tot_eta_y*(ii+1)/nSubPixels) ii++;
hhy[ib+ibx*nbeta]=hiy[ibx]/tot_eta_y; hhy[ib+iby*nbetaX]=hiy[iby]/tot_eta_y;
} }
} }
} }
for (int ib=0; ib<nbetaY; ib++) {
for (int ibx=0; ibx<nbetaX; ibx++) {
etax=etamin+ibx*etastepX;
//cout << etax << endl;
if (etax>=0 && etax<=1)
hx[ibx]=heta[ibx+ib*nbetaX];
else {
hx[ibx]=0;
}
}
hix[0]=hx[0];
for (int ibx=1; ibx<nbetaX; ibx++) {
hix[ibx]=hix[ibx-1]+hx[ibx];
}
tot_eta_x=hix[nbetaX-1]+1;
for (int ibx=0; ibx<nbetaX; ibx++) {
if (tot_eta_x<=0) {
hhx[ibx+ib*nbetaX]=-1;
}
else {
hhx[ibx+ib*nbetaX]=hix[ibx]/tot_eta_x;
}
}
for (int ibx=0; ibx<nbetaX; ibx++) {
if (tot_eta_x<=0) {
hhx[ibx+ib*nbetaX]=-1;
} else {
//if (hix[ibx]>tot_eta_x*(ii+1)/nSubPixels) ii++;
hhx[ibx+ib*nbetaX]=hix[ibx]/tot_eta_x;
}
}
}
int ibx, iby, ib; int ibx, iby, ib;
iby=0; iby=0;
while (hhx[iby*nbeta+nbeta/2]<0) iby++; while (hhx[iby*nbetaY+nbetaY/2]<0) iby++;
for (ib=0; ib<iby;ib++) { for (ib=0; ib<iby;ib++) {
for (ibx=0; ibx<nbeta;ibx++) for (ibx=0; ibx<nbetaX;ibx++)
hhx[ibx+nbeta*ib]=hhx[ibx+nbeta*iby]; hhx[ibx+nbetaX*ib]=hhx[ibx+nbetaX*iby];
} }
iby=nbeta-1; iby=nbetaY-1;
while (hhx[iby*nbeta+nbeta/2]<0) iby--; while (hhx[iby*nbetaY+nbetaY/2]<0) iby--;
for (ib=iby+1; ib<nbeta;ib++) { for (ib=iby+1; ib<nbetaY;ib++) {
for (ibx=0; ibx<nbeta;ibx++) for (ibx=0; ibx<nbetaX;ibx++)
hhx[ibx+nbeta*ib]=hhx[ibx+nbeta*iby]; hhx[ibx+nbetaX*ib]=hhx[ibx+nbetaX*iby];
} }
iby=0; iby=0;
while (hhy[nbeta/2*nbeta+iby]<0) iby++; while (hhy[nbetaX/2*nbetaX+iby]<0) iby++;
for (ib=0; ib<iby;ib++) { for (ib=0; ib<iby;ib++) {
for (ibx=0; ibx<nbeta;ibx++) for (ibx=0; ibx<nbetaY;ibx++)
hhy[ib+nbeta*ibx]=hhy[iby+nbeta*ibx]; hhy[ib+nbetaX*ibx]=hhy[iby+nbetaX*ibx];
} }
iby=nbeta-1; iby=nbetaX-1;
while (hhy[nbeta/2*nbeta+iby]<0) iby--; while (hhy[nbetaX/2*nbetaX+iby]<0) iby--;
for (ib=iby+1; ib<nbeta;ib++) { for (ib=iby+1; ib<nbetaX;ib++) {
for (ibx=0; ibx<nbeta;ibx++) for (ibx=0; ibx<nbetaY;ibx++)
hhy[ib+nbeta*ibx]=hhy[iby+nbeta*ibx]; hhy[ib+nbetaX*ibx]=hhy[iby+nbetaX*ibx];
} }
@ -156,9 +170,16 @@ class etaInterpolationPosXY : public virtual etaInterpolationBase{
}; };
class eta2InterpolationPosXY : public virtual eta2InterpolationBase, public virtual etaInterpolationPosXY { class eta2InterpolationPosXY : public virtual eta2InterpolationBase, public virtual etaInterpolationPosXY {
public: public:
eta2InterpolationPosXY(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny,ns, nb, emin,emax),eta2InterpolationBase(nx,ny,ns, nb, emin,emax),etaInterpolationPosXY(nx,ny,ns, nb, emin,emax){ eta2InterpolationPosXY(int nx=400, int ny=400, int ns=25, int nsy=25, int nb=-1, int nby=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny, ns, nsy, nb, nby, emin, emax),eta2InterpolationBase(nx,ny, ns, nsy, nb, nby, emin, emax),etaInterpolationPosXY(nx,ny, ns, nsy, nb, nby, emin, emax){
// cout << "e2pxy " << nb << " " << emin << " " << emax << endl; // cout << "e2pxy " << nb << " " << emin << " " << emax << endl;
}; };
@ -172,8 +193,8 @@ class eta2InterpolationPosXY : public virtual eta2InterpolationBase, public virt
class eta3InterpolationPosXY : public virtual eta3InterpolationBase, public virtual etaInterpolationPosXY { class eta3InterpolationPosXY : public virtual eta3InterpolationBase, public virtual etaInterpolationPosXY {
public: public:
eta3InterpolationPosXY(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny,ns, nb, emin,emax),eta3InterpolationBase(nx,ny,ns, nb, emin,emax), etaInterpolationPosXY(nx,ny,ns, nb, emin,emax){ eta3InterpolationPosXY(int nx=400, int ny=400, int ns=25, int nsy=25, int nb=-1, int nby=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny, ns, nsy, nb, nby, emin, emax),eta3InterpolationBase(nx,ny, ns, nsy, nb, nby, emin, emax), etaInterpolationPosXY(nx,ny, ns, nsy, nb, nby, emin, emax){
cout << "e3pxy " << nbeta << " " << etamin << " " << etamax << " " << nSubPixels<< endl; // cout << "e3pxy " << nbeta << " " << etamin << " " << etamax << " " << nSubPixels<< endl;
}; };
eta3InterpolationPosXY(eta3InterpolationPosXY *orig): etaInterpolationBase(orig), etaInterpolationPosXY(orig) {}; eta3InterpolationPosXY(eta3InterpolationPosXY *orig): etaInterpolationBase(orig), etaInterpolationPosXY(orig) {};

View File

@ -29,19 +29,21 @@ class slsInterpolation
{ {
public: public:
slsInterpolation(int nx=400, int ny=400, int ns=25) :nPixelsX(nx), nPixelsY(ny), nSubPixels(ns), id(0) { slsInterpolation(int nx=400, int ny=400, int ns=25, int nsy=-1) :nPixelsX(nx), nPixelsY(ny), nSubPixelsX(ns), nSubPixelsY(nsy),id(0) {
hint=new int[ns*nx*ns*ny]; if (nSubPixelsY<=0) nSubPixelsY=nSubPixelsX;
hint=new int[nSubPixelsX*nx*nSubPixelsY*ny];
}; };
slsInterpolation(slsInterpolation *orig){ slsInterpolation(slsInterpolation *orig){
nPixelsX=orig->nPixelsX; nPixelsX=orig->nPixelsX;
nPixelsY=orig->nPixelsY; nPixelsY=orig->nPixelsY;
nSubPixels=orig->nSubPixels; nSubPixelsX=orig->nSubPixelsX;
nSubPixelsY=orig->nSubPixelsY;
hint=new int[nSubPixels*nPixelsX*nSubPixels*nPixelsY]; hint=new int[nSubPixelsX*nPixelsX*nSubPixelsY*nPixelsY];
memcpy(hint, orig->hint,nSubPixels*nPixelsX*nSubPixels*nPixelsY*sizeof(int)); memcpy(hint, orig->hint,nSubPixelsX*nPixelsX*nSubPixelsY*nPixelsY*sizeof(int));
}; };
@ -51,23 +53,37 @@ class slsInterpolation
return new slsInterpolation(this); return new slsInterpolation(this);
}*/ }*/
int getNSubPixels() {return nSubPixels;}; int getNSubPixelsX() {return nSubPixelsX;};
int getNSubPixelsY() {return nSubPixelsY;};
int getNSubPixels() {if (nSubPixelsX==nSubPixelsY) return nSubPixelsX; else return 0;};
void getNSubPixels(int &nsx, int &nsy) {nsx=nSubPixelsX; nsy=nsx=nSubPixelsY;}
void setNSubPixels(int ns, int nsy=-1) {
int setNSubPixels(int ns) {
if (ns>0 && ns!=nSubPixels) {
delete [] hint; delete [] hint;
nSubPixels=ns; nSubPixelsX=ns;
hint=new int[nSubPixels*nPixelsX*nSubPixels*nPixelsY]; if (nsy>0) nSubPixelsY=nsy;
} else nSubPixelsY=ns;
return nSubPixels;
hint=new int[nSubPixelsX*nPixelsX*nSubPixelsY*nPixelsY];
//return nSubPixels;
} }
int getImageSize(int &nnx, int &nny, int &ns) {
nnx=nSubPixels*nPixelsX;
nny=nSubPixels*nPixelsY;
ns=nSubPixels; int getImageSize(int &nnx, int &nny, int &nsx, int &nsy) {
return nSubPixels*nSubPixels*nPixelsX*nPixelsY; nnx=nSubPixelsX*nPixelsX;
nny=nSubPixelsY*nPixelsY;
nsx=nSubPixelsX;
nsy=nSubPixelsY;
return nSubPixelsX*nSubPixelsY*nPixelsX*nPixelsY;
};
int getImageSize() {
return nSubPixelsX*nSubPixelsY*nPixelsX*nPixelsY;
}; };
@ -92,14 +108,14 @@ class slsInterpolation
//cout << "!" <<endl; //cout << "!" <<endl;
float *gm=NULL; float *gm=NULL;
int *dummy=getInterpolatedImage(); int *dummy=getInterpolatedImage();
gm=new float[ nSubPixels* nSubPixels* nPixelsX*nPixelsY]; gm=new float[ nSubPixelsX* nSubPixelsY* nPixelsX*nPixelsY];
if (gm) { if (gm) {
for (int ix=0; ix<nPixelsX*nSubPixels; ix++) { for (int ix=0; ix<nPixelsX*nSubPixelsX; ix++) {
for (int iy=0; iy<nPixelsY*nSubPixels; iy++) { for (int iy=0; iy<nPixelsY*nSubPixelsY; iy++) {
gm[iy*nPixelsX*nSubPixels+ix]=dummy[iy*nPixelsX*nSubPixels+ix]; gm[iy*nPixelsX*nSubPixelsX+ix]=dummy[iy*nPixelsX*nSubPixelsX+ix];
} }
} }
WriteToTiff(gm, imgname,nSubPixels* nPixelsX ,nSubPixels* nPixelsY); WriteToTiff(gm, imgname,nSubPixelsY* nPixelsX ,nSubPixelsY* nPixelsY);
delete [] gm; delete [] gm;
} else cout << "Could not allocate float image " << endl; } else cout << "Could not allocate float image " << endl;
return NULL; return NULL;
@ -120,9 +136,9 @@ class slsInterpolation
virtual void clearInterpolatedImage() { virtual void clearInterpolatedImage() {
for (int ix=0; ix<nPixelsX*nSubPixels; ix++) { for (int ix=0; ix<nPixelsX*nSubPixelsX; ix++) {
for (int iy=0; iy<nPixelsY*nSubPixels; iy++) { for (int iy=0; iy<nPixelsY*nSubPixelsY; iy++) {
hint[iy*nPixelsX*nSubPixels+ix]=0; hint[iy*nPixelsX*nSubPixelsX+ix]=0;
} }
} }
@ -133,11 +149,11 @@ class slsInterpolation
virtual int *addToImage(double int_x, double int_y){ virtual int *addToImage(double int_x, double int_y){
int iy=((double)nSubPixels)*int_y; int iy=((double)nSubPixelsY)*int_y;
int ix=((double)nSubPixels)*int_x; int ix=((double)nSubPixelsX)*int_x;
if (ix>=0 && ix<(nPixelsX*nSubPixels) && iy<(nSubPixels*nPixelsY) && iy>=0 ){ if (ix>=0 && ix<(nPixelsX*nSubPixelsX) && iy<(nSubPixelsY*nPixelsY) && iy>=0 ){
// cout << int_x << " " << int_y << " " << " " << ix << " " << iy << " " << ix+iy*nPixelsX*nSubPixels << " " << hint[ix+iy*nPixelsX*nSubPixels]; // cout << int_x << " " << int_y << " " << " " << ix << " " << iy << " " << ix+iy*nPixelsX*nSubPixels << " " << hint[ix+iy*nPixelsX*nSubPixels];
(*(hint+ix+iy*nPixelsX*nSubPixels))+=1; (*(hint+ix+iy*nPixelsX*nSubPixelsX))+=1;
// cout << " " << hint[ix+iy*nPixelsX*nSubPixels] << endl; // cout << " " << hint[ix+iy*nPixelsX*nSubPixels] << endl;
}// else }// else
// cout << "bad! "<< int_x << " " << int_y << " " << " " << ix << " " << iy << " " << ix+iy*nPixelsX*nSubPixels << endl; // cout << "bad! "<< int_x << " " << int_y << " " << " " << ix << " " << iy << " " << ix+iy*nPixelsX*nSubPixels << endl;
@ -180,11 +196,12 @@ class slsInterpolation
/* cluster[2]=cl+6; */ /* cluster[2]=cl+6; */
sum=0; sum=0;
int xoff=0, yoff=0;
#ifndef WRITE_QUAD
double sumBL=0; double sumBL=0;
double sumTL=0; double sumTL=0;
double sumBR=0; double sumBR=0;
double sumTR=0; double sumTR=0;
int xoff=0, yoff=0;
for (int ix=0; ix<3; ix++) { for (int ix=0; ix<3; ix++) {
for (int iy=0; iy<3; iy++) { for (int iy=0; iy<3; iy++) {
sum+=cl[ix+3*iy]; sum+=cl[ix+3*iy];
@ -204,34 +221,95 @@ class slsInterpolation
if(sumTL >= totquad){ if(sumTL >= totquad){
/* sDum[0][0] = cluster[1][0]; sDum[1][0] = cluster[2][0]; */ /* #ifdef WRITE_QUAD */
/* sDum[0][1] = cluster[1][1]; sDum[1][1] = cluster[2][1]; */ /* /\* sDum[0][0] = cluster[1][0]; sDum[1][0] = cluster[2][0]; *\/ */
/* /\* sDum[0][1] = cluster[1][1]; sDum[1][1] = cluster[2][1]; *\/ */
/* if (sumTL ==sum) { */
/* #endif */
corner = TOP_LEFT; corner = TOP_LEFT;
totquad=sumTL; totquad=sumTL;
xoff=0; xoff=0;
yoff=1; yoff=1;
/* #ifdef WRITE_QUAD */
/* } */
/* #endif */
} }
if(sumBR >= totquad){ if(sumBR >= totquad){
/* sDum[0][0] = cluster[0][1]; sDum[1][0] = cluster[1][1]; */ /* sDum[0][0] = cluster[0][1]; sDum[1][0] = cluster[1][1]; */
/* sDum[0][1] = cluster[0][2]; sDum[1][1] = cluster[1][2]; */ /* sDum[0][1] = cluster[0][2]; sDum[1][1] = cluster[1][2]; */
/* #ifdef WRITE_QUAD */
/* /\* sDum[0][0] = cluster[1][0]; sDum[1][0] = cluster[2][0]; *\/ */
/* /\* sDum[0][1] = cluster[1][1]; sDum[1][1] = cluster[2][1]; *\/ */
/* if (sumBR ==sum) { */
/* #endif */
corner = BOTTOM_RIGHT;
xoff=1; xoff=1;
yoff=0; yoff=0;
corner = BOTTOM_RIGHT;
totquad=sumBR; totquad=sumBR;
/* #ifdef WRITE_QUAD */
/* } */
/* #endif */
} }
if(sumTR >= totquad){ if(sumTR >= totquad){
/* #ifdef WRITE_QUAD */
/* /\* sDum[0][0] = cluster[1][0]; sDum[1][0] = cluster[2][0]; *\/ */
/* /\* sDum[0][1] = cluster[1][1]; sDum[1][1] = cluster[2][1]; *\/ */
/* if (sumTR ==sum) { */
/* #endif */
xoff=1; xoff=1;
yoff=1; yoff=1;
/* sDum[0][0] = cluster[1][1]; sDum[1][0] = cluster[2][1]; */ /* sDum[0][0] = cluster[1][1]; sDum[1][0] = cluster[2][1]; */
/* sDum[0][1] = cluster[1][2]; sDum[1][1] = cluster[2][2]; */ /* sDum[0][1] = cluster[1][2]; sDum[1][1] = cluster[2][2]; */
corner = TOP_RIGHT; corner = TOP_RIGHT;
totquad=sumTR; totquad=sumTR;
/* #ifdef WRITE_QUAD */
/* } */
/* #endif */
} }
#endif
#ifdef WRITE_QUAD
double sumB=0;
double sumT=0;
double sumR=0;
double sumL=0;
for (int ix=0; ix<3; ix++) {
for (int iy=0; iy<3; iy++) {
sum+=cl[ix+3*iy];
if (ix<1 ) sumL+=cl[ix+iy*3];
if (ix>1) sumR+=cl[ix+iy*3];
if (iy<1) sumB=cl[ix+iy*3];
if (iy>1) sumT+=cl[ix+iy*3];
}
}
totquad=sum;
if ( sumT==0 && sumR==0) {
corner = BOTTOM_LEFT;
xoff=0;
yoff=0;
} else if ( sumB==0 && sumR==0 ) {
corner = TOP_LEFT;
xoff=0;
yoff=1;
} else if ( sumT==0 && sumL==0) {
corner = BOTTOM_RIGHT;
xoff=1;
yoff=0;
} else if ( sumB==0 && sumL==0) {
xoff=1;
yoff=1;
corner = TOP_RIGHT;
} else
printf("** bad 2x2 cluster!\n");
#endif
for (int ix=0; ix<2; ix++) { for (int ix=0; ix<2; ix++) {
for (int iy=0; iy<2; iy++) { for (int iy=0; iy<2; iy++) {
@ -495,7 +573,7 @@ class slsInterpolation
protected: protected:
int nPixelsX, nPixelsY; int nPixelsX, nPixelsY;
int nSubPixels; int nSubPixelsX, nSubPixelsY;
int id; int id;
int *hint; int *hint;
}; };

View File

@ -1,6 +1,6 @@
INCDIR= -I. -I../dataStructures ../tiffIO.cpp -I../ -I../interpolations/ -I../../slsSupportLib/include/ -I../../slsReceiverSoftware/include/ -I../../libs/rapidjson/ INCDIR= -I. -I../dataStructures ../tiffIO.cpp -I../ -I../interpolations/ -I../../slsSupportLib/include/ -I../../slsReceiverSoftware/include/ -I../../libs/rapidjson/
LDFLAG= -L/usr/lib64/ -lpthread -lm -lstdc++ -lzmq -pthread -lrt -ltiff -O3 -g -std=c++11 -Wall LDFLAG= -L/usr/lib64/ -lpthread -lm -lstdc++ -lzmq -pthread -lrt -ltiff -O3 -g -std=c++11 -Wall -L../../build/bin/ -lSlsSupport
#-L../../bin -lhdf5 -L. #-L../../bin -lhdf5 -L.
#DESTDIR?=../bin #DESTDIR?=../bin

View File

@ -7,7 +7,12 @@
#include "sls_detector_defs.h" #include "sls_detector_defs.h"
#include "ZmqSocket.h" #include "ZmqSocket.h"
#ifndef RECT
#include "moench03T1ZmqDataNew.h" #include "moench03T1ZmqDataNew.h"
#endif
#ifdef RECT
#include "moench03T1ZmqDataNewRect.h"
#endif
#include "moench03GhostSummation.h" #include "moench03GhostSummation.h"
#include "moench03CommonMode.h" #include "moench03CommonMode.h"
#include <vector> #include <vector>
@ -49,10 +54,11 @@ int main(int argc, char *argv[]) {
*/ */
FILE *of=NULL; FILE *of=NULL;
int fifosize=5000; int fifosize=5000;
int etabins=1000;//nsubpix*2*100; int etabins=1000, etabinsy=1000;//nsubpix*2*100;
double etamin=-1, etamax=2; double etamin=-1, etamax=2;
int nSubPixels=2; int nSubPixelsX=2;
// int emin, emax; // int emin, emax;
int nSubPixelsY=2;
// help // help
if (argc < 3 ) { 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"); cprintf(RED, "Help: ./trial [receive socket ip] [receive starting port number] [send_socket ip] [send starting port number] [nthreads] [nsubpix] [gainmap] [etafile]\n");
@ -76,12 +82,6 @@ int main(int argc, char *argv[]) {
time_t begin,end,finished; time_t begin,end,finished;
int rms=0; int rms=0;
int nped=1000, nped0=100;
#ifdef PTC
nped=10000;
nped0=10000;
#endif
if (argc > 4) { if (argc > 4) {
@ -101,10 +101,14 @@ int main(int argc, char *argv[]) {
nthreads=atoi(argv[5]); nthreads=atoi(argv[5]);
cout << "Number of threads is: " << nthreads << endl; cout << "Number of threads is: " << nthreads << endl;
if (argc>6) if (argc>6) {
nSubPixels=atoi(argv[6]); nSubPixelsX=atoi(argv[6]);
cout << "Number of subpixels is: " << nSubPixels << endl; nSubPixelsY=nSubPixelsX;
#ifdef RECT
nSubPixelsX=2;
#endif
}
cout << "Number of subpixels is: " << nSubPixelsX << " " << nSubPixelsY << endl;
char *gainfname=NULL; char *gainfname=NULL;
if (argc>7) { if (argc>7) {
@ -119,12 +123,12 @@ int main(int argc, char *argv[]) {
} }
//slsDetectorData *det=new moench03T1ZmqDataNew(); //slsDetectorData *det=new moench03T1ZmqDataNew();
moench03T1ZmqDataNew *det=new moench03T1ZmqDataNew(5000,sizeof(int)); moench03T1ZmqDataNew *det=new moench03T1ZmqDataNew();
cout << endl << " det" <<endl; cout << endl << " det" <<endl;
int npx, npy; int npx, npy;
det->getDetectorSize(npx, npy); det->getDetectorSize(npx, npy);
int send_something=0;
int maxSize = npx*npy*2;//32*2*8192;//5000;//atoi(argv[3]); int maxSize = npx*npy*2;//32*2*8192;//5000;//atoi(argv[3]);
@ -137,15 +141,8 @@ int main(int argc, char *argv[]) {
int ncol_cm=CM_ROWS; int ncol_cm=CM_ROWS;
double xt_ghost=C_GHOST; double xt_ghost=C_GHOST;
moench03CommonMode *cm=new moench03CommonMode(ncol_cm);
moench03GhostSummation *gs=new moench03GhostSummation(det, xt_ghost);
moench03CommonMode *cm=NULL;
moench03GhostSummation *gs=NULL;
#ifdef CORR
cm=new moench03CommonMode(ncol_cm);
gs=new moench03GhostSummation(det, xt_ghost);
#endif
double *gainmap=NULL; double *gainmap=NULL;
float *gm; float *gm;
double *gmap=NULL; double *gmap=NULL;
@ -173,18 +170,18 @@ int main(int argc, char *argv[]) {
//analogDetector<uint16_t> *filter=new analogDetector<uint16_t>(det,1,NULL,1000); //analogDetector<uint16_t> *filter=new analogDetector<uint16_t>(det,1,NULL,1000);
#ifndef INTERP #ifndef INTERP
singlePhotonDetector *filter=new singlePhotonDetector(det,3, nSigma, 1, cm, nped, nped0, -1, -1, gainmap, gs); singlePhotonDetector *filter=new singlePhotonDetector(det,3, nSigma, 1, cm, 1000, 10, -1, -1, gainmap, gs);
multiThreadedCountingDetector *mt=new multiThreadedCountingDetector(filter,nthreads,fifosize); multiThreadedCountingDetector *mt=new multiThreadedCountingDetector(filter,nthreads,fifosize);
// multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize); // multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize);
#endif #endif
#ifdef INTERP #ifdef INTERP
eta2InterpolationPosXY *interp=new eta2InterpolationPosXY(npx, npy, nSubPixels, etabins, etamin, etamax); eta2InterpolationPosXY *interp=new eta2InterpolationPosXY(npx, npy, nSubPixelsX,nSubPixelsY, etabins, etabinsy, etamin, etamax);
if (etafname) interp->readFlatField(etafname); if (etafname) interp->readFlatField(etafname);
interpolatingDetector *filter=new interpolatingDetector(det,interp, nSigma, 1, cm, nped, nped0, -1, -1, gainmap, gs); interpolatingDetector *filter=new interpolatingDetector(det,interp, nSigma, 1, cm, 1000, 10, -1, -1, gainmap, gs);
multiThreadedInterpolatingDetector *mt=new multiThreadedInterpolatingDetector(filter,nthreads,fifosize); multiThreadedInterpolatingDetector *mt=new multiThreadedInterpolatingDetector(filter,nthreads,fifosize);
#endif #endif
@ -275,7 +272,17 @@ int main(int argc, char *argv[]) {
// header variables // header variables
uint64_t acqIndex = -1; uint64_t acqIndex = -1;
uint64_t frameIndex = -1; uint64_t frameIndex = -1;
//uint32_t subFrameIndex = -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; uint64_t fileindex = -1;
string filename = ""; string filename = "";
// char* image = new char[size]; // char* image = new char[size];
@ -285,10 +292,10 @@ int main(int argc, char *argv[]) {
int iframe=0; int iframe=0;
char ofname[10000]; char ofname[10000];
char fname[10000]; string fname;
// int length; // int length;
int *detimage; int *detimage=NULL;
int nnx, nny,nns; int nnx, nny,nnsx, nnsy;
//uint32_t imageSize = 0, nPixelsX = 0, nPixelsY = 0, //uint32_t imageSize = 0, nPixelsX = 0, nPixelsY = 0,
//uint32_t dynamicRange = 0; //uint32_t dynamicRange = 0;
// infinite loop // infinite loop
@ -304,12 +311,11 @@ int main(int argc, char *argv[]) {
//int16_t *dout;//=new int16_t [nnx*nny]; //int16_t *dout;//=new int16_t [nnx*nny];
uint32_t dr = 32; uint32_t dr = 32;
int32_t *dout=NULL;//=new int32_t [nnx*nny]; int32_t *dout=NULL;//=new int32_t [nnx*nny];
float *doutf=NULL;//=new int32_t [nnx*nny];
uint16_t roundRNumber = 0; uint16_t roundRNumber = 0;
uint8_t detType = 0; uint8_t detType = 0;
uint8_t version = 0; uint8_t version = 0;
// int* flippedData = 0; string additionalJsonHeader="" ;
string* additionalJsonHeader = 0;
//char* additionalJsonHeader = 0;
int32_t threshold=0; int32_t threshold=0;
@ -327,7 +333,7 @@ int main(int argc, char *argv[]) {
frameMode fMode=eFrame; frameMode fMode=eFrame;
double *ped; double *ped;
filter->getImageSize(nnx, nny,nns); filter->getImageSize(nnx, nny,nnsx, nnsy);
@ -368,27 +374,60 @@ int main(int argc, char *argv[]) {
cprintf(RED, "Sent Dummy\n"); cprintf(RED, "Sent Dummy\n");
} }
} else { } else {
send_something=0;
if (fMode==ePedestal) { if (fMode==ePedestal) {
sprintf(ofname,"%s_%ld_ped.tiff",fname,fileindex); sprintf(ofname,"%s_%ld_ped.tiff",fname.c_str(),fileindex);
mt->writePedestal(ofname); mt->writePedestal(ofname);
cout << "Writing pedestal to " << ofname << endl; cout << "Writing pedestal to " << ofname << endl;
if (rms){ if (rms){
sprintf(ofname,"%s_%ld_var.tiff",fname,fileindex); sprintf(ofname,"%s_%ld_var.tiff",fname.c_str(),fileindex);
mt->writePedestalRMS(ofname); mt->writePedestalRMS(ofname);
cout << "Writing pedestal variance to " << ofname << endl;
} }
send_something=1;
} }
#ifdef INTERP #ifdef INTERP
else if (fMode==eFlat) { else if (fMode==eFlat) {
mt->prepareInterpolation(ok); mt->prepareInterpolation(ok);
sprintf(ofname,"%s_%ld_eta.tiff",fname,fileindex); sprintf(ofname,"%s_%ld_eta.tiff",fname.c_str(),fileindex);
mt->writeFlatField(ofname); mt->writeFlatField(ofname);
cout << "Writing eta to " << ofname << endl; cout << "Writing eta to " << ofname << endl;
send_something=1;
} }
#endif #endif
else { else {
sprintf(ofname,"%s_%ld.tiff",fname,fileindex); 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; ix<nnx*nny; ix++) {
doutf[ix]=detimage[ix]*subnorm;
if (doutf[ix]<0) doutf[ix]=0;
}
cout << "Writing image to " << ofname << endl;
WriteToTiff(doutf,ofname ,nnx, nny);
if (doutf)
delete [] doutf;
doutf=NULL;
nsubframes++;
insubframe=0;
send_something=1;
}
} else {
sprintf(ofname,"%s_%ld.tiff",fname.c_str(),fileindex);
mt->writeImage(ofname); mt->writeImage(ofname);
send_something=1;
}
cout << "Writing image to " << ofname << endl; cout << "Writing image to " << ofname << endl;
} }
// cout << nns*nnx*nny*nns*dr/8 << " " << length << endl; // cout << nns*nnx*nny*nns*dr/8 << " " << length << endl;
@ -397,11 +436,13 @@ int main(int argc, char *argv[]) {
if (fMode==ePedestal) { if (fMode==ePedestal) {
cprintf(MAGENTA,"Get pedestal!\n"); cprintf(MAGENTA,"Get pedestal!\n");
nns=1; nnsx=1;
nnsy=1;
nnx=npx; nnx=npx;
nny=npy; nny=npy;
//dout= new int16_t[nnx*nny*nns*nns]; //dout= new int16_t[nnx*nny*nns*nns];
dout= new int32_t[nnx*nny*nns*nns]; dout= new int32_t[nnx*nny*nnsx*nnsy];
// cout << "get pedestal " << endl; // cout << "get pedestal " << endl;
ped=mt->getPedestal(); ped=mt->getPedestal();
// cout << "got pedestal " << endl; // cout << "got pedestal " << endl;
@ -409,7 +450,7 @@ int main(int argc, char *argv[]) {
dout[ix]=ped[ix]; dout[ix]=ped[ix];
// if (ix<100*400) // if (ix<100*400)
// cout << ix << " " << ped[ix] << " "<< dout[ix] << endl; // cout << ix << " " << ped[ix] << endl;
} }
} }
@ -427,15 +468,19 @@ int main(int argc, char *argv[]) {
} }
#endif #endif
else { else {
detimage=mt->getImage(nnx,nny,nns); detimage=mt->getImage(nnx,nny,nnsx, nnsy);
cprintf(MAGENTA,"Get image!\n"); cprintf(MAGENTA,"Get image!\n");
cout << nnx << " " << nny << " " << nns << endl; cout << nnx << " " << nny << " " << nnsx << " " << nnsy << endl;
// nns=1; // nns=1;
// nnx=npx; // nnx=npx;
// nny=npy; // nny=npy;
// nnx=nnx*nns; // nnx=nnx*nns;
//nny=nny*nns; //nny=nny*nns;
dout= new int32_t[nnx*nny]; 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; ix<nnx*nny; ix++) { for (int ix=0; ix<nnx*nny; ix++) {
// for (int iy=0; iy<nny*nns; iy++) { // for (int iy=0; iy<nny*nns; iy++) {
// for (int isx=0; isx<nns; isx++) { // for (int isx=0; isx<nns; isx++) {
@ -447,46 +492,22 @@ int main(int argc, char *argv[]) {
// } // }
// } // }
dout[ix]=detimage[ix]; dout[ix]=detimage[ix]*subnorm;
if (dout[ix]<0) dout[ix]=0; if (dout[ix]<0) dout[ix]=0;
// cout << ix << " " << dout[ix] << endl; // cout << ix << " " << dout[ix] << endl;
// } // }
} }
} }
//if ((insubframe>0 && subframes>0) || (subframes<=0) ){
if(send_something) {
//// int SendHeaderData ( int index, bool dummy, uint32_t jsonversion, uint32_t dynamicrange = 0, uint64_t fileIndex = 0, 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);
// uint32_t ndetx = 0, uint32_t ndety = 0, uint32_t npixelsx = 0, uint32_t npixelsy = 0, uint32_t imageSize = 0,
// uint64_t acqIndex = 0, uint64_t fIndex = 0, const char* fname = NULL,
// uint64_t frameNumber = 0, uint32_t expLength = 0, uint32_t packetNumber = 0,
// uint64_t bunchId = 0, uint64_t timestamp = 0,
// uint16_t modId = 0, uint16_t row = 0, uint16_t column = 0, uint16_t reserved = 0,
// uint32_t debug = 0, uint16_t roundRNumber = 0,
// uint8_t detType = 0, uint8_t version = 0, int gapPixelsEnable = 0, int flippedDataX = 0,
// char* additionalJsonHeader = 0) {
// cout << "Sending image size " << nnx << " " << nny << endl;
// #ifndef DEVELOPER
// #ifndef MOENCH_BRANCH
// 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, additionalJsonHeader);
// #endif
// #endif
// #ifdef DEVELOPER
// #ifdef CTBGUI
// zmqsocket2->SendHeaderData (0, false,SLS_DETECTOR_JSON_HEADER_VERSION , dr, fileindex, 0,0,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);
// #endif
// #ifndef CTBGUI
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);
//#endif
zmqsocket2->SendData((char*)dout,nnx*nny*dr/8); zmqsocket2->SendData((char*)dout,nnx*nny*dr/8);
cprintf(GREEN, "Sent Data\n"); cprintf(GREEN, "Sent Data\n");
}
zmqsocket2->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION); zmqsocket2->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION);
cprintf(RED, "Sent Dummy\n"); cprintf(RED, "Sent Dummy\n");
@ -541,8 +562,8 @@ int main(int argc, char *argv[]) {
//dataSize=size; //dataSize=size;
strcpy(fname,filename.c_str()); //strcpy(fname,filename.c_str());
fname=filename;
// cprintf(BLUE, "Header Info:\n" // cprintf(BLUE, "Header Info:\n"
// "size: %u\n" // "size: %u\n"
// "multisize: %u\n" // "multisize: %u\n"
@ -585,7 +606,7 @@ int main(int argc, char *argv[]) {
if (frameMode_s == "pedestal"){ if (frameMode_s == "pedestal"){
fMode=ePedestal; fMode=ePedestal;
//isPedestal=1; //isPedestal=1;
} else if (frameMode_s == "newpedestal"){ } else if (frameMode_s == "newPedestal"){
mt->newDataSet(); //resets pedestal mt->newDataSet(); //resets pedestal
// cprintf(MAGENTA, "Resetting pedestal\n"); // cprintf(MAGENTA, "Resetting pedestal\n");
fMode=ePedestal; fMode=ePedestal;
@ -594,8 +615,8 @@ int main(int argc, char *argv[]) {
mt->newDataSet(); //resets pedestal mt->newDataSet(); //resets pedestal
// cprintf(MAGENTA, "Resetting pedestal\n"); // cprintf(MAGENTA, "Resetting pedestal\n");
fMode=ePedestal; fMode=ePedestal;
rms=1;
//isPedestal=1; //isPedestal=1;
rms=1;
} }
#ifdef INTERP #ifdef INTERP
else if (frameMode_s == "flatfield") { else if (frameMode_s == "flatfield") {
@ -616,13 +637,11 @@ int main(int argc, char *argv[]) {
frameMode_s="frame"; frameMode_s="frame";
} }
} }
}
cprintf(MAGENTA, "%s\n" , frameMode_s.c_str()); cprintf(MAGENTA, "%s\n" , frameMode_s.c_str());
} else
cprintf(RED, "%s\n" , frameMode_s.c_str());
mt->setFrameMode(fMode); mt->setFrameMode(fMode);
threshold=0; // threshold=0;
cprintf(MAGENTA, "Threshold: "); cprintf(MAGENTA, "Threshold: ");
if (doc.HasMember("threshold")) { if (doc.HasMember("threshold")) {
if (doc["threshold"].IsInt()) { if (doc["threshold"].IsInt()) {
@ -690,11 +709,10 @@ int main(int argc, char *argv[]) {
} }
} }
cprintf(MAGENTA, "%s\n" , detectorMode_s.c_str()); }
} else
cprintf(RED, "%s\n" , frameMode_s.c_str());
mt->setDetectorMode(dMode); mt->setDetectorMode(dMode);
cprintf(MAGENTA, "%s\n" , detectorMode_s.c_str());
// cout << "done " << endl; // cout << "done " << endl;
@ -737,6 +755,21 @@ int main(int argc, char *argv[]) {
// mt->setNSubPixels(nSubPixels); // 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; newFrame=0;
/* zmqsocket->CloseHeaderMessage();*/ /* zmqsocket->CloseHeaderMessage();*/
@ -745,7 +778,7 @@ int main(int argc, char *argv[]) {
// cout << "file" << endl; // cout << "file" << endl;
// cout << "data " << endl; // cout << "data " << endl;
if (of==NULL && dMode!=eAnalog && fMode!=ePedestal && threshold<=0) { if (of==NULL) {
#ifdef WRITE_QUAD #ifdef WRITE_QUAD
sprintf(ofname,"%s_%ld.clust2",filename.c_str(),fileindex); sprintf(ofname,"%s_%ld.clust2",filename.c_str(),fileindex);
#endif #endif
@ -762,10 +795,13 @@ int main(int argc, char *argv[]) {
} }
// cout << "data" << endl; // cout << "data" << endl;
// get data // get data
// acqIndex = doc["acqIndex"].GetUint64(); // acqIndex = doc["acqIndex"].GetUint64();
frameIndex = doc["fIndex"].GetUint64(); frameIndex = doc["fIndex"].GetUint64();
// subFrameIndex = doc["expLength"].GetUint(); // subFrameIndex = doc["expLength"].GetUint();
// bunchId=doc["bunchId"].GetUint(); // bunchId=doc["bunchId"].GetUint();
@ -774,19 +810,81 @@ int main(int argc, char *argv[]) {
// cout << acqIndex << " " << frameIndex << " " << subFrameIndex << " "<< bunchId << " " << timestamp << " " << packetNumber << endl; // cout << acqIndex << " " << frameIndex << " " << subFrameIndex << " "<< bunchId << " " << timestamp << " " << packetNumber << endl;
if (packetNumber>=40) { if (packetNumber>=40) {
//*((int*)buff)=frameIndex; //*((int*)buff)=frameIndex;
if (insubframe==0) f0=frameIndex;
memcpy(buff,&frameIndex,sizeof(int)); memcpy(buff,&frameIndex,sizeof(int));
//length = //length =
zmqsocket->ReceiveData(0, buff+sizeof(int), size); zmqsocket->ReceiveData(0, buff+sizeof(int), size);
mt->pushData(buff); mt->pushData(buff);
mt->nextThread(); mt->nextThread();
mt->popFree(buff); mt->popFree(buff);
cprintf(GREEN, "Frame\n"); insubframe++;
nsubframes=frameIndex+1-f0;
} else { } else {
cprintf(RED, "Incomplete frame: received only %d packet\n", packetNumber); cprintf(RED, "Incomplete frame: received only %d packet\n", packetNumber);
//length = //length =
zmqsocket->ReceiveData(0, dummybuff, size); 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; ix<nnx*nny; ix++) {
dout[ix]=detimage[ix]*subnorm;
if (dout[ix]<0) dout[ix]=0;
doutf[ix]=dout[ix];
}
sprintf(ofname,"%s_sf%ld_%ld.tiff",fname.c_str(),nnsubframe,fileindex);
cout << "Writing image to " << ofname << endl;
WriteToTiff(doutf,ofname ,nnx, nny);
nsubframes++;
insubframe=0;
nnsubframe++;
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);
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++; iframe++;
} // exiting infinite loop } // exiting infinite loop

View File

@ -95,7 +95,7 @@ public:
virtual int *getImage() { virtual int *getImage() {
return det->getImage(); return det->getImage();
} }
virtual int getImageSize(int &nnx, int &nny, int &ns) {return det->getImageSize(nnx, nny, ns);}; virtual int getImageSize(int &nnx, int &nny, int &ns, int &nsy) {return det->getImageSize(nnx, nny, ns, nsy);};
virtual int getDetectorSize(int &nnx, int &nny) {return det->getDetectorSize(nnx, nny);}; virtual int getDetectorSize(int &nnx, int &nny) {return det->getDetectorSize(nnx, nny);};
virtual ~threadedAnalogDetector() {StopThread(); delete fifoFree; delete fifoData;} virtual ~threadedAnalogDetector() {StopThread(); delete fifoFree; delete fifoData;}
@ -221,10 +221,10 @@ FILE *getFilePointer(){return det->getFilePointer();};
} }
virtual int setNSubPixels(int ns) { virtual int setNSubPixels(int ns, int nsy) {
slsInterpolation *interp=(det)->getInterpolation(); slsInterpolation *interp=(det)->getInterpolation();
if (interp) return interp->setNSubPixels(ns); if (interp) interp->setNSubPixels(ns, nsy);
else return 1;}; return 1;};
virtual slsInterpolation *setInterpolation(slsInterpolation *f){ virtual slsInterpolation *setInterpolation(slsInterpolation *f){
@ -314,11 +314,11 @@ public:
virtual void newDataSet(){for (int i=0; i<nThreads; i++) dets[i]->newDataSet();}; virtual void newDataSet(){for (int i=0; i<nThreads; i++) dets[i]->newDataSet();};
virtual int *getImage(int &nnx, int &nny, int &ns) { virtual int *getImage(int &nnx, int &nny, int &ns, int &nsy) {
int *img; int *img;
// int nnx, nny, ns; // int nnx, nny, ns;
// int nnx, nny, ns; // int nnx, nny, ns;
int nn=dets[0]->getImageSize(nnx, nny,ns); int nn=dets[0]->getImageSize(nnx, nny,ns, nsy);
if (image) { if (image) {
delete image; delete image;
image=NULL; image=NULL;
@ -363,10 +363,10 @@ public:
/* dets[ii]->writeImage(tit); */ /* dets[ii]->writeImage(tit); */
/* } */ /* } */
/* #endif */ /* #endif */
int nnx, nny, ns; int nnx, nny, ns, nsy;
getImage(nnx, nny, ns); getImage(nnx, nny, ns,nsy);
//int nnx, nny, ns; //int nnx, nny, ns;
int nn=dets[0]->getImageSize(nnx, nny, ns); int nn=dets[0]->getImageSize(nnx, nny, ns, nsy);
float *gm=new float[nn]; float *gm=new float[nn];
if (gm) { if (gm) {
for (int ix=0; ix<nn; ix++) { for (int ix=0; ix<nn; ix++) {

View File

@ -46,7 +46,7 @@ public:
}; };
virtual int setNSubPixels(int ns) { return (dets[0])->setNSubPixels(ns);}; /* virtual int setNSubPixels(int ns) { return (dets[0])->setNSubPixels(ns);}; */
virtual void resetFlatField() {(dets[0])->resetFlatField();}; virtual void resetFlatField() {(dets[0])->resetFlatField();};
@ -69,13 +69,13 @@ public:
virtual int *getImage(int &nnx, int &nny, int &ns) { virtual int *getImage(int &nnx, int &nny, int &nsx, int &nsy) {
if (getInterpolation()==NULL) return multiThreadedAnalogDetector::getImage(nnx,nny,ns); if (getInterpolation()==NULL) return multiThreadedAnalogDetector::getImage(nnx,nny,nsx, nsy);
//if one interpolates, the whole image is stored in detector 0; //if one interpolates, the whole image is stored in detector 0;
int *img; int *img;
// int nnx, nny, ns; // int nnx, nny, ns;
// int nnx, nny, ns; // int nnx, nny, ns;
int nn=dets[0]->getImageSize(nnx, nny,ns); int nn=dets[0]->getImageSize(nnx, nny,nsx, nsy);
if (image) { if (image) {
delete image; delete image;
image=NULL; image=NULL;

View File

@ -366,9 +366,9 @@ int *getClusters(char *data, int *ph=NULL) {
if (cm) if (cm) {
addToCommonMode(data); addToCommonMode(data);
}
for (int iy=ymin; iy<ymax; iy++) { for (int iy=ymin; iy<ymax; iy++) {
for (int ix=xmin; ix<xmax; ix++) { for (int ix=xmin; ix<xmax; ix++) {