mirror of
https://github.com/slsdetectorgroup/slsDetectorPackage.git
synced 2025-04-23 15:00:02 +02:00
370 lines
8.3 KiB
C++
370 lines
8.3 KiB
C++
#ifndef ETA_INTERPOLATION_BASE_H
|
|
#define ETA_INTERPOLATION_BASE_H
|
|
|
|
#ifdef MYROOT1
|
|
#include <TObject.h>
|
|
#include <TTree.h>
|
|
#include <TH2D.h>
|
|
#include <TH2F.h>
|
|
#endif
|
|
|
|
#include "slsInterpolation.h"
|
|
#include "tiffIO.h"
|
|
|
|
class etaInterpolationBase : public slsInterpolation {
|
|
|
|
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) {
|
|
// cout << "eb " << nb << " " << emin << " " << emax << endl;
|
|
// cout << nb << " " << etamin << " " << etamax << endl;
|
|
if (nbeta<=0) {
|
|
//cout << "aaa:" <<endl;
|
|
nbeta=nSubPixels*10;
|
|
}
|
|
if (etamin>=etamax) {
|
|
etamin=-1;
|
|
etamax=2;
|
|
}
|
|
etastep=(etamax-etamin)/nbeta;
|
|
heta=new int[nbeta*nbeta];
|
|
hhx=new float[nbeta*nbeta];
|
|
hhy=new float[nbeta*nbeta];
|
|
rangeMin=etamin;
|
|
rangeMax=etamax;
|
|
flat= new double[nSubPixels*nSubPixels];
|
|
hintcorr=new int [nSubPixels*nSubPixels*nPixelsX*nPixelsY];
|
|
|
|
};
|
|
|
|
etaInterpolationBase(etaInterpolationBase *orig): slsInterpolation(orig){
|
|
nbeta=orig->nbeta;
|
|
etamin=orig->etamin;
|
|
etamax=orig->etamax;
|
|
rangeMin=orig->rangeMin;
|
|
rangeMax=orig->rangeMax;
|
|
|
|
|
|
etastep=(etamax-etamin)/nbeta;
|
|
heta=new int[nbeta*nbeta];
|
|
memcpy(heta,orig->heta,nbeta*nbeta*sizeof(int));
|
|
hhx=new float[nbeta*nbeta];
|
|
memcpy(hhx,orig->hhx,nbeta*nbeta*sizeof(float));
|
|
hhy=new float[nbeta*nbeta];
|
|
memcpy(hhy,orig->hhy,nbeta*nbeta*sizeof(float));
|
|
hintcorr=new int [nSubPixels*nSubPixels*nPixelsX*nPixelsY];
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
|
|
virtual void resetFlatField() {
|
|
for (int ibx=0; ibx<nbeta*nbeta; ibx++) {
|
|
heta[ibx]=0;
|
|
hhx[ibx]=0;
|
|
hhy[ibx]=0;
|
|
}
|
|
|
|
|
|
};
|
|
|
|
int *setEta(int *h, int nb=-1, double emin=1, double emax=0)
|
|
{
|
|
if (h) {
|
|
if (heta) delete [] heta;
|
|
heta=h;
|
|
nbeta=nb;
|
|
if (nb<=0) nbeta=nSubPixels*10;
|
|
etamin=emin;
|
|
etamax=emax;
|
|
if (etamin>=etamax) {
|
|
etamin=-1;
|
|
etamax=2;
|
|
}
|
|
rangeMin=etamin;
|
|
rangeMax=etamax;
|
|
etastep=(etamax-etamin)/nbeta;
|
|
}
|
|
return heta;
|
|
};
|
|
|
|
int *setFlatField(int *h, int nb=-1, double emin=1, double emax=0)
|
|
{
|
|
return setEta(h, nb, emin, emax);
|
|
};
|
|
|
|
|
|
|
|
int *getFlatField(){return setEta(NULL);};
|
|
|
|
int *getFlatField(int &nb, double &emin, double &emax){
|
|
nb=nbeta;
|
|
emin=etamin;
|
|
emax=etamax;
|
|
return getFlatField();
|
|
};
|
|
|
|
|
|
void *writeFlatField(const char * imgname) {
|
|
float *gm=NULL;
|
|
gm=new float[nbeta*nbeta];
|
|
for (int ix=0; ix<nbeta; ix++) {
|
|
for (int iy=0; iy<nbeta; iy++) {
|
|
gm[iy*nbeta+ix]=heta[iy*nbeta+ix];
|
|
}
|
|
}
|
|
WriteToTiff(gm, imgname, nbeta, nbeta);
|
|
delete [] gm;
|
|
return NULL;
|
|
};
|
|
|
|
int readFlatField(const char * imgname, double emin=1, double emax=0) {
|
|
if (emax>=1) etamax=emax;
|
|
if (emin<=0) etamin=emin;
|
|
|
|
if (etamin>=etamax) {
|
|
etamin=-1;
|
|
etamax=2;
|
|
}
|
|
|
|
etastep=(etamax-etamin)/nbeta;
|
|
uint32 nnx;
|
|
uint32 nny;
|
|
float *gm=ReadFromTiff(imgname, nnx, nny);
|
|
if (nnx!=nny) {
|
|
cout << "different number of bins in x " << nnx << " and y " << nny<< " !"<< endl;
|
|
cout << "Aborting read"<< endl;
|
|
return 0;
|
|
}
|
|
nbeta=nnx;
|
|
if (gm) {
|
|
if (heta) {
|
|
delete [] heta;
|
|
delete [] hhx;
|
|
delete [] hhy;
|
|
}
|
|
|
|
heta=new int[nbeta*nbeta];
|
|
hhx=new float[nbeta*nbeta];
|
|
hhy=new float[nbeta*nbeta];
|
|
|
|
for (int ix=0; ix<nbeta; ix++) {
|
|
for (int iy=0; iy<nbeta; iy++) {
|
|
heta[iy*nbeta+ix]=gm[iy*nbeta+ix];
|
|
}
|
|
}
|
|
delete [] gm;
|
|
return 1;
|
|
}
|
|
return 0;
|
|
};
|
|
|
|
|
|
|
|
|
|
|
|
float *gethhx()
|
|
{
|
|
// hhx->Scale((double)nSubPixels);
|
|
return hhx;
|
|
};
|
|
|
|
float *gethhy()
|
|
{
|
|
// hhy->Scale((double)nSubPixels);
|
|
return hhy;
|
|
};
|
|
virtual int addToFlatField(double etax, double etay){
|
|
int ex,ey;
|
|
ex=(etax-etamin)/etastep;
|
|
ey=(etay-etamin)/etastep;
|
|
if (ey<nbeta && ex<nbeta && ex>=0 && ey>=0)
|
|
heta[ey*nbeta+ex]++;
|
|
return 0;
|
|
};
|
|
|
|
|
|
// virtual void prepareInterpolation(int &ok)=0;
|
|
|
|
|
|
void debugSaveAll(int ind=0) {
|
|
int ibx, iby;
|
|
char tit[10000];
|
|
|
|
float tot_eta=0;
|
|
|
|
float *etah=new float[nbeta*nbeta];
|
|
int etabins=nbeta;
|
|
int ibb=0;
|
|
|
|
for (int ii=0; ii<etabins*etabins; ii++) {
|
|
|
|
etah[ii]=heta[ii];
|
|
tot_eta+=heta[ii];
|
|
}
|
|
sprintf(tit,"/scratch/eta_%d.tiff",ind);
|
|
WriteToTiff(etah, tit, etabins, etabins);
|
|
|
|
|
|
for (int ii=0; ii<etabins*etabins; ii++) {
|
|
ibb=(hhx[ii]*nSubPixels);
|
|
etah[ii]=ibb;
|
|
}
|
|
sprintf(tit,"/scratch/eta_hhx_%d.tiff",ind);
|
|
WriteToTiff(etah, tit, etabins, etabins);
|
|
|
|
for (int ii=0; ii<etabins*etabins; ii++) {
|
|
ibb=hhy[ii]*nSubPixels;
|
|
etah[ii]=ibb;
|
|
}
|
|
sprintf(tit,"/scratch/eta_hhy_%d.tiff",ind);
|
|
WriteToTiff(etah, tit, etabins, etabins);
|
|
|
|
|
|
float *ftest=new float[nSubPixels*nSubPixels];
|
|
|
|
for (int ib=0; ib<nSubPixels*nSubPixels; ib++) ftest[ib]=0;
|
|
|
|
|
|
//int ibx=0, iby=0;
|
|
|
|
for (int ii=0; ii<nbeta*nbeta; ii++) {
|
|
|
|
ibx=nSubPixels*hhx[ii];
|
|
iby=nSubPixels*hhy[ii];
|
|
if (ibx<0) ibx=0;
|
|
if (iby<0) iby=0;
|
|
if (ibx>=nSubPixels) ibx=nSubPixels-1;
|
|
if (iby>=nSubPixels) iby=nSubPixels-1;
|
|
|
|
|
|
if (ibx>=0 && ibx<nSubPixels && iby>=0 && iby<nSubPixels) {
|
|
//
|
|
// if (ibx>0 && iby>0) cout << ibx << " " << iby << " " << ii << endl;
|
|
ftest[ibx+iby*nSubPixels]+=heta[ii];
|
|
} else
|
|
cout << "Bad interpolation "<< ii << " " << ibx << " " << iby<< endl;
|
|
|
|
}
|
|
|
|
sprintf(tit,"/scratch/ftest_%d.tiff",ind);
|
|
WriteToTiff(ftest, tit, nSubPixels, nSubPixels);
|
|
|
|
//int ibx=0, iby=0;
|
|
tot_eta/=nSubPixels*nSubPixels;
|
|
int nbad=0;
|
|
for (int ii=0; ii<etabins*etabins; ii++) {
|
|
ibx=nSubPixels*hhx[ii];
|
|
iby=nSubPixels*hhy[ii];
|
|
if (ftest[ibx+iby*nSubPixels]<tot_eta*0.5) {
|
|
etah[ii]=1;
|
|
nbad++;
|
|
} else if(ftest[ibx+iby*nSubPixels]>tot_eta*2.){
|
|
etah[ii]=2;
|
|
nbad++;
|
|
} else
|
|
etah[ii]=0;
|
|
}
|
|
sprintf(tit,"/scratch/eta_bad_%d.tiff",ind);
|
|
WriteToTiff(etah, tit, etabins, etabins);
|
|
// cout << "Index: " << ind << "\t Bad bins: "<< nbad << endl;
|
|
//int ibx=0, iby=0;
|
|
|
|
delete [] ftest;
|
|
delete [] etah;
|
|
|
|
}
|
|
|
|
|
|
protected:
|
|
|
|
|
|
double calcDiff(double avg, float *hx, float *hy) {
|
|
//double p_tot=0;
|
|
double diff=0, d;
|
|
//double bsize=1./nSubPixels;
|
|
int nbad=0;
|
|
double p_tot_x[nSubPixels], p_tot_y[nSubPixels], p_tot[nSubPixels*nSubPixels];
|
|
double maxdiff=0, mindiff=avg*nSubPixels*nSubPixels;
|
|
|
|
int ipx, ipy;
|
|
for (ipy=0; ipy<nSubPixels; ipy++) {
|
|
for (ipx=0; ipx<nSubPixels; ipx++) {
|
|
p_tot[ipx+ipy*nSubPixels]=0;
|
|
}
|
|
p_tot_y[ipy]=0;
|
|
p_tot_x[ipy]=0;
|
|
}
|
|
|
|
for (int ibx=0; ibx<nbeta; ibx++) {
|
|
for (int iby=0; iby<nbeta; iby++) {
|
|
ipx=hx[ibx+iby*nbeta]*nSubPixels;
|
|
if (ipx<0) ipx=0;
|
|
if (ipx>=nSubPixels) ipx=nSubPixels-1;
|
|
|
|
ipy=hy[ibx+iby*nbeta]*nSubPixels;
|
|
if (ipy<0) ipy=0;
|
|
if (ipy>=nSubPixels) ipy=nSubPixels-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];
|
|
|
|
|
|
}
|
|
}
|
|
|
|
|
|
// cout << endl << endl;
|
|
for (ipy=0; ipy<nSubPixels; ipy++) {
|
|
cout.width(5);
|
|
//flat_y[ipy]=p_tot_y[ipy];//avg/nSubPixels;
|
|
for (ipx=0; ipx<nSubPixels; ipx++) {
|
|
|
|
// flat_x[ipx]=p_tot_x[ipx];///avg/nSubPixels;
|
|
flat[ipx+nSubPixels*ipy]=p_tot[ipx+nSubPixels*ipy];///avg;
|
|
d=p_tot[ipx+nSubPixels*ipy]-avg;
|
|
if (d<0) d*=-1.;
|
|
if (d>5*sqrt(avg) )
|
|
nbad++;
|
|
diff+=d*d;
|
|
if (d<mindiff) mindiff=d;
|
|
if (d>maxdiff) maxdiff=d;
|
|
// cout << setprecision(4) << p_tot[ipx+nSubPixels*ipy] << " ";
|
|
}
|
|
|
|
/* cout << "** " << setprecision(4) << flat_y[ipy]; */
|
|
//cout << "\n";
|
|
}
|
|
/* cout << "**" << endl; cout.width(5); */
|
|
/* for (ipx=0; ipx<nSubPixels; ipx++) { */
|
|
/* cout << setprecision(4) << flat_x[ipx] << " "; */
|
|
/* } */
|
|
//cout << "**" << endl; cout.width(5);
|
|
//cout << "Min diff: " << mindiff/sqrt(avg) << " Max diff: " << maxdiff/sqrt(avg) << " Nbad: " << nbad << endl;
|
|
|
|
|
|
|
|
// cout << "Bad pixels: " << 100.*(float)nbad/((float)(nSubPixels*nSubPixels)) << " %" << endl;
|
|
return sqrt(diff);
|
|
}
|
|
|
|
float *hhx;
|
|
float *hhy;
|
|
int *heta;
|
|
int nbeta;
|
|
double etamin, etamax, etastep;
|
|
double rangeMin, rangeMax;
|
|
|
|
|
|
|
|
double *flat;
|
|
int *hintcorr;
|
|
|
|
|
|
};
|
|
|
|
#endif
|