detectors/slsDetectorPackage
25
0
Fork 0
mirror of https://github.com/slsdetectorgroup/slsDetectorPackage.git synced 2025-04-25 07:40:03 +02:00
Code Issues Actions Packages Releases Activity

Improved the interpolation with adaptive bins and implemented flat field correction when saving the interpolated images

Browse Source
This commit is contained in:
bergamaschi 2018-12-14 18:11:12 +01:00
parent 9d8aa5fe91
commit 563d644d2f
9 changed files with 885 additions and 592 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

181
slsDetectorCalibration/interpolations/eta2InterpolationBase.h
View File

@ -14,47 +14,18 @@ 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 nb=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny, ns, nb, emin, emax) {
// cout << "e2ib " << nb << " " << emin << " " << emax << endl;
/* if (nbeta<=0) { */ /* if (etamin>=etamax) { */
/* nbeta=nSubPixels*10; */ /* etamin=-1; */
/* etamax=2; */
/* // cout << ":" <<endl; */
/* } */ /* } */
if (etamin>=etamax) { /* etastep=(etamax-etamin)/nbeta; */
etamin=-1;
etamax=2;
// cout << ":" <<endl;
}
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;
}; };
eta2InterpolationBase(eta2InterpolationBase *orig): etaInterpolationBase(orig){ }; eta2InterpolationBase(eta2InterpolationBase *orig): etaInterpolationBase(orig){ };
/* virtual eta2InterpolationBase* Clone()=0; {
return new eta2InterpolationBase(this);
};
*/
////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////
//////////// /*It return position hit for the event in input */ ////////////// //////////// /*It return position hit for the event in input */ //////////////
@ -100,12 +71,7 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cl,double &int_x, double &int_y) { virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cl,double &int_x, double &int_y) {
double cc[2][2]; double cc[2][2];
double *cluster[3];
int xoff, yoff; int xoff, yoff;
cluster[0]=cl;
cluster[1]=cl+3;
cluster[2]=cl+6;
switch (quad) { switch (quad) {
case BOTTOM_LEFT: case BOTTOM_LEFT:
xoff=0; xoff=0;
@ -128,10 +94,10 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
} }
double etax, etay; double etax, etay;
if (nSubPixels>2) { if (nSubPixels>2) {
cc[0][0]=cluster[yoff][xoff]; cc[0][0]=cl[xoff+3*yoff];
cc[1][0]=cluster[yoff+1][xoff]; cc[1][0]=cl[xoff+3*(yoff+1)];
cc[0][1]=cluster[yoff][xoff+1]; cc[0][1]=cl[xoff+1+3*yoff];
cc[1][1]=cluster[yoff+1][xoff+1]; cc[1][1]=cl[xoff+1+3*(yoff+1)];
calcEta(totquad,cc,etax,etay); calcEta(totquad,cc,etax,etay);
} }
return getInterpolatedPosition(x,y,etax, etay,quad,int_x,int_y); return getInterpolatedPosition(x,y,etax, etay,quad,int_x,int_y);
@ -143,11 +109,7 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,int *cl,double &int_x, double &int_y) { virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,int *cl,double &int_x, double &int_y) {
double cc[2][2]; double cc[2][2];
int *cluster[3];
int xoff, yoff; int xoff, yoff;
cluster[0]=cl;
cluster[1]=cl+3;
cluster[2]=cl+6;
switch (quad) { switch (quad) {
case BOTTOM_LEFT: case BOTTOM_LEFT:
@ -171,10 +133,10 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
} }
double etax, etay; double etax, etay;
if (nSubPixels>2) { if (nSubPixels>2) {
cc[0][0]=cluster[yoff][xoff]; cc[0][0]=cl[xoff+3*yoff];
cc[1][0]=cluster[yoff+1][xoff]; cc[1][0]=cl[xoff+3*(yoff+1)];
cc[0][1]=cluster[yoff][xoff+1]; cc[0][1]=cl[xoff+1+3*yoff];
cc[1][1]=cluster[yoff+1][xoff+1]; cc[1][1]=cl[xoff+1+3*(xoff+1)];
calcEta(totquad,cc,etax,etay); calcEta(totquad,cc,etax,etay);
} }
return getInterpolatedPosition(x,y,etax, etay,quad,int_x,int_y); return getInterpolatedPosition(x,y,etax, etay,quad,int_x,int_y);
@ -221,60 +183,46 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
if (nSubPixels>2) { if (nSubPixels>2) {
#ifdef MYROOT1 ex=(etax-etamin)/etastep;
xpos_eta=(hhx->GetBinContent(hhx->GetXaxis()->FindBin(etax),hhy->GetYaxis()->FindBin(etay)))/((double)nSubPixels); ey=(etay-etamin)/etastep;
ypos_eta=(hhy->GetBinContent(hhx->GetXaxis()->FindBin(etax),hhy->GetYaxis()->FindBin(etay)))/((double)nSubPixels); if (ex<0) {
#endif cout << "x*"<< ex << endl;
#ifndef MYROOT1 ex=0;
ex=(etax-etamin)/etastep; }
ey=(etay-etamin)/etastep; if (ex>=nbeta) {
if (ex<0) { cout << "x?"<< ex << endl;
cout << "x*"<< ex << endl; ex=nbeta-1;
ex=0; }
} if (ey<0) {
if (ex>=nbeta) { cout << "y*"<< ey << endl;
cout << "x?"<< ex << endl; ey=0;
ex=nbeta-1; }
if (ey>=nbeta) {
} cout << "y?"<< ey << endl;
if (ey<0) { ey=nbeta-1;
cout << "y*"<< ey << endl; }
ey=0;
}
if (ey>=nbeta) {
cout << "y?"<< ey << endl;
ey=nbeta-1;
}
xpos_eta=(((double)hhx[(ey*nbeta+ex)]))+dX ;///((double)nSubPixels); xpos_eta=(((double)hhx[(ey*nbeta+ex)]))+dX ;///((double)nSubPixels);
ypos_eta=(((double)hhy[(ey*nbeta+ex)]))+dY ;///((double)nSubPixels); ypos_eta=(((double)hhy[(ey*nbeta+ex)]))+dY ;///((double)nSubPixels);
//else
//return 0;
#endif
} else { } else {
xpos_eta=0.5*dX+0.25; xpos_eta=0.5*dX+0.25;
ypos_eta=0.5*dY+0.25; ypos_eta=0.5*dY+0.25;
} }
int_x=((double)x) + xpos_eta+0.5; int_x=((double)x) + xpos_eta+0.5;
int_y=((double)y) + ypos_eta+0.5; int_y=((double)y) + ypos_eta+0.5;
} }
virtual int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay) { virtual int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay) {
double cc[2][2]; double cc[2][2];
int *cluster[3];
int xoff, yoff; int xoff, yoff;
cluster[0]=cl;
cluster[1]=cl+3;
cluster[2]=cl+6;
switch (quad) { switch (quad) {
case BOTTOM_LEFT: case BOTTOM_LEFT:
@ -296,17 +244,11 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
default: default:
; ;
} }
cc[0][0]=cluster[yoff][xoff]; cc[0][0]=cl[xoff+3*yoff];
cc[1][0]=cluster[yoff+1][xoff]; cc[1][0]=cl[xoff+3*(yoff+1)];
cc[0][1]=cluster[yoff][xoff+1]; cc[0][1]=cl[xoff+1+3*yoff];
cc[1][1]=cluster[yoff+1][xoff+1]; cc[1][1]=cl[xoff+1+3*(yoff+1)];
/* cout << cl[0] << " " << cl[1] << " " << cl[2] << endl; */
/* cout << cl[3] << " " << cl[4] << " " << cl[5] << endl; */
/* cout << cl[6] << " " << cl[7] << " " << cl[8] << endl; */
/* cout <<"******"<<totquad << " " << quad << endl; */
/* cout << cc[0][0]<< " " << cc[0][1] << endl; */
/* cout << cc[1][0]<< " " << cc[1][1] << endl; */
//calcMyEta(totquad,quad,cl,etax, etay); //calcMyEta(totquad,quad,cl,etax, etay);
calcEta(totquad, cc,etax, etay); calcEta(totquad, cc,etax, etay);
@ -318,11 +260,7 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
virtual int addToFlatField(double totquad,int quad,double *cl,double &etax, double &etay) { virtual int addToFlatField(double totquad,int quad,double *cl,double &etax, double &etay) {
double cc[2][2]; double cc[2][2];
double *cluster[3];
int xoff, yoff; int xoff, yoff;
cluster[0]=cl;
cluster[1]=cl+3;
cluster[2]=cl+6;
switch (quad) { switch (quad) {
case BOTTOM_LEFT: case BOTTOM_LEFT:
@ -344,10 +282,10 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
default: default:
; ;
} }
cc[0][0]=cluster[yoff][xoff]; cc[0][0]=cl[xoff+3*yoff];
cc[1][0]=cluster[yoff+1][xoff]; cc[1][0]=cl[(yoff+1)*3+xoff];
cc[0][1]=cluster[yoff][xoff+1]; cc[0][1]=cl[yoff*3+xoff+1];
cc[1][1]=cluster[yoff+1][xoff+1]; cc[1][1]=cl[(yoff+1)*3+xoff+1];
/* cout << cl[0] << " " << cl[1] << " " << cl[2] << endl; */ /* cout << cl[0] << " " << cl[1] << " " << cl[2] << endl; */
/* cout << cl[3] << " " << cl[4] << " " << cl[5] << endl; */ /* cout << cl[3] << " " << cl[4] << " " << cl[5] << endl; */
@ -414,6 +352,37 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
return 0; return 0;
}; };
virtual int *getInterpolatedImage(){
int ipx, ipy;
// cout << "ff" << endl;
calcDiff(1, hhx, hhy); //get flat
double avg=0;
for (ipx=0; ipx<nSubPixels; ipx++)
for (ipy=0; ipy<nSubPixels; ipy++)
avg+=flat[ipx+ipy*nSubPixels];
avg/=nSubPixels*nSubPixels;
for (int ibx=0 ; ibx<nSubPixels*nPixelsX; ibx++) {
ipx=ibx%nSubPixels-nSubPixels/2;
if (ipx<0) ipx=nSubPixels+ipx;
for (int iby=0 ; iby<nSubPixels*nPixelsY; iby++) {
ipy=iby%nSubPixels-nSubPixels/2;
if (ipy<0) ipy=nSubPixels+ipy;
// cout << ipx << " " << ipy << " " << ibx << " " << iby << endl;
if (flat[ipx+ipy*nSubPixels]>0)
hintcorr[ibx+iby*nSubPixels*nPixelsX]=hint[ibx+iby*nSubPixels*nPixelsX]*(avg/flat[ipx+ipy*nSubPixels]);
else
hintcorr[ibx+iby*nSubPixels*nPixelsX]=hint[ibx+iby*nSubPixels*nPixelsX];
}
}
return hintcorr;
};
/* protected: */ /* protected: */
/* #ifdef MYROOT1 */ /* #ifdef MYROOT1 */

308
slsDetectorCalibration/interpolations/etaInterpolationAdaptiveBins.h
View File

@ -1,131 +1,133 @@
#ifndef ETA_INTERPOLATION_ADAPTIVEBINS_H #ifndef ETA_INTERPOLATION_ADAPTIVEBINS_H
#define ETA_INTERPOLATION_ADAPTIVEBINS_H #define ETA_INTERPOLATION_ADAPTIVEBINS_H
#include <cmath>
#include "tiffIO.h" #include "tiffIO.h"
//#include "etaInterpolationBase.h" //#include "etaInterpolationBase.h"
#include "etaInterpolationPosXY.h" #include "etaInterpolationPosXY.h"
class etaInterpolationAdaptiveBins : public etaInterpolationPosXY { class etaInterpolationAdaptiveBins : public etaInterpolationPosXY {
// protected:
private: private:
double calcDiff(double avg, float *hx, float *hy) {
double p_tot=0;
double diff=0;
double bsize=1./nSubPixels;
for (int ipx=0; ipx<nSubPixels; ipx++) { virtual void iterate(float *newhhx, float *newhhy) {
for (int ipy=0; ipy<nSubPixels; ipy++) {
p_tot=0;
for (int ibx=0; ibx<nbeta; ibx++) {
for (int iby=0; iby<nbeta; iby++) {
if ( hx[ibx+iby*nbeta]>=((ipx)*bsize) && hx[ibx+iby*nbeta]<((ipx+1)*bsize) && hy[ibx+iby*nbeta]>=((ipy)*bsize) && hy[ibx+iby*nbeta]<((ipy+1)*bsize)) {
p_tot+=heta[ibx+iby*nbeta];
}
}
}
cout << p_tot << " \t ";
diff+=(p_tot-avg)*(p_tot-avg);
}
cout << "\n";
}
return diff;
}
void iterate(float *newhhx, float *newhhy) {
double bsize=1./nSubPixels; double bsize=1./nSubPixels;
double hy[nbeta]; //profile y double hy[nSubPixels][nbeta]; //profile y
double hx[nbeta]; //profile x double hx[nSubPixels][nbeta]; //profile x
double hix[nbeta]; //integral of projection x double hix[nSubPixels][nbeta]; //integral of projection x
double hiy[nbeta]; //integral of projection y double hiy[nSubPixels][nbeta]; //integral of projection y
int ipy, ipx;
double tot_eta_x[nSubPixels];
double tot_eta_y[nSubPixels];
//for (int ipy=0; ipy<nSubPixels; ipy++) {
double tot_eta_x=0; for (ipy=0; ipy<nSubPixels; ipy++) {
double tot_eta_y=0; for (int ibx=0; ibx<nbeta; ibx++) {
for (int ipy=0; ipy<nSubPixels; ipy++) { hx[ipy][ibx]=0;
hy[ipy][ibx]=0;
for (int ibx=0; ibx<nbeta; ibx++) {
hx[ibx]=0;
hy[ibx]=0;
} }
}
tot_eta_x=0;
tot_eta_y=0;
// cout << ipy << " " << ((ipy)*bsize) << " " << ((ipy+1)*bsize) << endl; // cout << ipy << " " << ((ipy)*bsize) << " " << ((ipy+1)*bsize) << endl;
for (int ibx=0; ibx<nbeta; ibx++) { for (int ibx=0; ibx<nbeta; ibx++) {
for (int iby=0; iby<nbeta; iby++) { for (int iby=0; iby<nbeta; iby++) {
ipy=hhy[ibx+iby*nbeta]*nSubPixels;
if (hhy[ibx+iby*nbeta]>=((ipy)*bsize) && hhy[ibx+iby*nbeta]<=((ipy+1)*bsize)) { if (ipy<0) ipy=0;
hx[ibx]+=heta[ibx+iby*nbeta]; if (ipy>=nSubPixels) ipy=nSubPixels-1;
tot_eta_x+=heta[ibx+iby*nbeta]; hx[ipy][ibx]+=heta[ibx+iby*nbeta];
}
if (hhx[ibx+iby*nbeta]>=((ipy)*bsize) && hhx[ibx+iby*nbeta]<=((ipy+1)*bsize)) { ipx=hhx[ibx+iby*nbeta]*nSubPixels;
hy[iby]+=heta[ibx+iby*nbeta]; if (ipx<0) ipx=0;
tot_eta_y+=heta[ibx+iby*nbeta]; if (ipx>=nSubPixels) ipx=nSubPixels-1;
} hy[ipx][iby]+=heta[ibx+iby*nbeta];
} }
} }
hix[0]=hx[0];
hiy[0]=hy[0];
for (int ib=1; ib<nbeta; ib++) {
hix[ib]=hix[ib-1]+hx[ib];
hiy[ib]=hiy[ib-1]+hy[ib];
}
// tot_eta_x=hix[nbeta-1];
// tot_eta_y=hiy[nbeta-1];
/* cout << "ipx " << ipy << " x: " << tot_eta_x << " " << hix[10]<< " " << hix[nbeta-1] << endl; */
/* cout << "ipy " << ipy << " y: " << tot_eta_y << " " << hiy[10]<< " " << hiy[nbeta-1] << endl; */
// for (int ipy=0; ipy<nSubPixels; ipy++) { for (ipy=0; ipy<nSubPixels; ipy++) {
hix[ipy][0]=hx[ipy][0];
hiy[ipy][0]=hy[ipy][0];
for (int ib=1; ib<nbeta; ib++) {
hix[ipy][ib]=hix[ipy][ib-1]+hx[ipy][ib];
hiy[ipy][ib]=hiy[ipy][ib-1]+hy[ipy][ib];
}
tot_eta_x[ipy]=hix[ipy][nbeta-1]+1;
tot_eta_y[ipy]=hiy[ipy][nbeta-1]+1;
// cout << ipy << " " << tot_eta_x[ipy] << " " << tot_eta_y[ipy] << endl;
}
// for (int ipy=0; ipy<nSubPixels; ipy++) {
for (int ibx=0; ibx<nbeta; ibx++) { for (int ibx=0; ibx<nbeta; ibx++) {
for (int iby=0; iby<nbeta; iby++) { for (int iby=0; iby<nbeta; iby++) {
if ( hhy[ibx+iby*nbeta]>=((ipy)*bsize) && hhy[ibx+iby*nbeta]<=((ipy+1)*bsize)) { // if ( hhy[ibx+iby*nbeta]>=((ipy)*bsize) && hhy[ibx+iby*nbeta]<=((ipy+1)*bsize)) {
newhhx[ibx+iby*nbeta]=hix[ibx]/((double)tot_eta_x); ipy=hhy[ibx+iby*nbeta]*nSubPixels;
if (newhhx[ibx+iby*nbeta]>1) cout << "***"<< ibx << " " << iby << newhhx[ibx+iby*nbeta] << endl;
if (ipy<0) ipy=0;
if (ipy>=nSubPixels) ipy=nSubPixels-1;
if (ipy>=0 && ipy<nSubPixels)
if (tot_eta_x[ipy]>0)
newhhx[ibx+iby*nbeta]=hix[ipy][ibx]/(tot_eta_x[ipy]);
else
cout << "Bad tot_etax " << ipy << " " << tot_eta_x[ipy] << endl;
else
cout << "** Bad value ipy " << ibx << " " << iby << " "<< ipy << " " << hhy[ibx+iby*nbeta]*nSubPixels << endl;
// if (newhhx[ibx+iby*nbeta]>=1 || newhhx[ibx+iby*nbeta]<0 ) cout << "***"<< ibx << " " << iby << newhhx[ibx+iby*nbeta] << endl;
// if (ipy==3 && ibx==10) cout << newhhx[ibx+iby*nbeta] << " " << hix[ibx] << " " << ibx+iby*nbeta << endl; // if (ipy==3 && ibx==10) cout << newhhx[ibx+iby*nbeta] << " " << hix[ibx] << " " << ibx+iby*nbeta << endl;
} // }
if (hhx[ibx+iby*nbeta]>=((ipy)*bsize) && hhx[ibx+iby*nbeta]<=((ipy+1)*bsize)) { ipy=hhx[ibx+iby*nbeta]*nSubPixels;
newhhy[ibx+iby*nbeta]=hiy[iby]/((double)tot_eta_y); //if (hhx[ibx+iby*nbeta]>=((ipy)*bsize) && hhx[ibx+iby*nbeta]<=((ipy+1)*bsize)) {
if (newhhy[ibx+iby*nbeta]>1) cout << "***"<< ibx << " " << iby << newhhy[ibx+iby*nbeta] << endl; if (ipy<0) ipy=0;
if (ipy>=nSubPixels) ipy=nSubPixels-1;
if (ipy>=0 && ipy<nSubPixels)
if (tot_eta_y[ipy]>0)
newhhy[ibx+iby*nbeta]=hiy[ipy][iby]/(tot_eta_y[ipy]);
else
cout << "Bad tot_etay " << ipy << " " << tot_eta_y[ipy] << endl;
else
cout << "** Bad value ipx " << ibx << " " << iby << " "<< ipy << " " << hhx[ibx+iby*nbeta]*nSubPixels << endl;
// if (newhhy[ibx+iby*nbeta]>=1 || newhhy[ibx+iby*nbeta]<0 ) cout << "***"<< ibx << " " << iby << newhhy[ibx+iby*nbeta] << endl;
// if (ipy==3 && iby==10) cout << newhhy[ibx+iby*nbeta] << " " << hiy[iby] << " " << ibx+iby*nbeta << endl; // if (ipy==3 && iby==10) cout << newhhy[ibx+iby*nbeta] << " " << hiy[iby] << " " << ibx+iby*nbeta << endl;
} // }
} }
} }
} // }
} }
public: public:
etaInterpolationAdaptiveBins(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationPosXY(nx,ny,ns, nb, emin,emax){}; etaInterpolationAdaptiveBins(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationPosXY(nx,ny,ns, nb, emin,emax){
// flat=new double[nSubPixels*nSubPixels]; flat_x=new double[nSubPixels]; flat_y=new double[nSubPixels];
// flat=new double[nSubPixels*nSubPixels];
};
etaInterpolationAdaptiveBins(etaInterpolationAdaptiveBins *orig): etaInterpolationPosXY(orig){}; etaInterpolationAdaptiveBins(etaInterpolationAdaptiveBins *orig): etaInterpolationPosXY(orig){hintcorr=new int[nPixelsX*nPixelsY*nSubPixels];};
virtual etaInterpolationAdaptiveBins* Clone() { virtual etaInterpolationAdaptiveBins* Clone()=0;
return new etaInterpolationAdaptiveBins(this); /* return new etaInterpolationAdaptiveBins(this); */
}; /* }; */
virtual void prepareInterpolation(int &ok) virtual void prepareInterpolation(int &ok, int nint=1000)
{ {
ok=1; ok=1;
cout << "Adaptive bins" << endl;
///*Eta Distribution Rebinning*/// ///*Eta Distribution Rebinning*///
double bsize=1./nSubPixels; //precision double bsize=1./nSubPixels; //precision
@ -153,65 +155,36 @@ class etaInterpolationAdaptiveBins : public etaInterpolationPosXY {
// } // }
// } // }
etaInterpolationPosXY::prepareInterpolation(ok); etaInterpolationPosXY::prepareInterpolation(ok);
#ifdef SAVE_ALL
char tit[10000];
float *etah=new float[nbeta*nbeta];
int etabins=nbeta;
for (int ii=0; ii<etabins*etabins; ii++) {
etah[ii]=hhx[ii];
}
sprintf(tit,"/scratch/start_hhx.tiff");
WriteToTiff(etah, tit, etabins, etabins);
for (int ii=0; ii<etabins*etabins; ii++) {
etah[ii]=hhy[ii];
}
sprintf(tit,"/scratch/start_hhy.tiff");
WriteToTiff(etah, tit, etabins, etabins);
#endif
int nint=1000;
double thr=1./((double)nSubPixels); double thr=1./((double)nSubPixels);
double avg=tot_eta/((double)(nSubPixels*nSubPixels)); double avg=tot_eta/((double)(nSubPixels*nSubPixels));
cout << "total eta entries is :"<< tot_eta << " avg: "<< avg << endl; double rms=sqrt(tot_eta);
cout << "Start " << endl; cout << "total eta entries is :"<< tot_eta << " avg: "<< avg << " rms: " << sqrt(tot_eta) << endl;
double old_diff=calcDiff(avg, hhx, hhy), new_diff=old_diff+1, best_diff=old_diff; double old_diff=calcDiff(avg, hhx, hhy), new_diff=old_diff+1, best_diff=old_diff;
cout << " diff= " << old_diff << endl; // cout << " chi2= " << old_diff << " (rms= " << sqrt(tot_eta) << ")" << endl;
cout << endl;
cout << endl;
debugSaveAll(0);
int iint=0; int iint=0;
float *newhhx=new float[nbeta*nbeta]; //profile x float *newhhx=new float[nbeta*nbeta]; //profile x
float *newhhy=new float[nbeta*nbeta]; //profile y float *newhhy=new float[nbeta*nbeta]; //profile y
float *besthhx=hhx; //profile x float *besthhx=hhx; //profile x
float *besthhy=hhy; //profile y float *besthhy=hhy; //profile y
while (iint<nint) {
cout << "Iteration " << iint << endl; cout << "Iteration "<< iint << " Chi2: " << old_diff << endl; //" Best: "<< best_diff << " RMS: "<< rms<< endl;
while (iint<nint && best_diff > rms) {
/* #ifdef SAVE_ALL */
/* if (iint%10==0) */
/* debugSaveAll(iint); */
/* #endif */
// cout << "Iteration " << iint << endl;
iterate(newhhx,newhhy); iterate(newhhx,newhhy);
new_diff=calcDiff(avg, newhhx, newhhy); new_diff=calcDiff(avg, newhhx, newhhy);
cout << " diff= " << new_diff << endl; // cout << " chi2= " << new_diff << " (rms= " << sqrt(tot_eta) << ")"<<endl;
/* #ifdef SAVE_ALL */
/* for (int ii=0; ii<etabins*etabins; ii++) { */
/* etah[ii]=newhhx[ii]; */
/* if (etah[ii]>1 || etah[ii]<0 ) cout << "***"<< ii << etah[ii] << endl; */
/* } */
/* sprintf(tit,"/scratch/neweta_hhx_%d.tiff",iint); */
/* WriteToTiff(etah, tit, etabins, etabins); */
/* for (int ii=0; ii<etabins*etabins; ii++) { */
/* etah[ii]=newhhy[ii]; */
/* if (etah[ii]>1 || etah[ii]<0 ) cout << "***"<< ii << etah[ii] << endl; */
/* } */
/* sprintf(tit,"/scratch/neweta_hhy_%d.tiff",iint); */
/* WriteToTiff(etah, tit, etabins, etabins); */
/* #endif */
if (new_diff<best_diff) { if (new_diff<best_diff) {
best_diff=new_diff; best_diff=new_diff;
besthhx=newhhx; besthhx=newhhx;
@ -225,19 +198,33 @@ class etaInterpolationAdaptiveBins : public etaInterpolationPosXY {
hhx=newhhx; hhx=newhhx;
hhy=newhhy; hhy=newhhy;
#ifdef SAVE_ALL
if (new_diff<=best_diff) {
debugSaveAll(iint);
}
#endif
newhhx=new float[nbeta*nbeta]; //profile x newhhx=new float[nbeta*nbeta]; //profile x
newhhy=new float[nbeta*nbeta]; //profile y newhhy=new float[nbeta*nbeta]; //profile y
old_diff=new_diff;
//} /* else { */
/* cout << "Difference not decreasing after "<< iint << " iterations (" << old_diff << " < " << new_diff << ")"<< endl; */
/* break; */
/* } */
iint++; /* if (new_diff<old_diff){ */
/* cout << "best difference at iteration "<< iint << " (" << new_diff << " < " << old_diff << ")"<< "Best: "<< best_diff << " RMS: "<< sqrt(tot_eta) << endl; */
/* ; */
/* } else { */
// break;
// }
old_diff=new_diff;
iint++;
cout << "Iteration "<< iint << " Chi2: " << new_diff << endl; //" Best: "<< best_diff << " RMS: "<< rms<< endl;
} }
delete [] newhhx; delete [] newhhx;
delete [] newhhy; delete [] newhhy;
if (hhx!=besthhx) if (hhx!=besthhx)
delete [] hhx; delete [] hhx;
@ -248,34 +235,47 @@ class etaInterpolationAdaptiveBins : public etaInterpolationPosXY {
hhy=besthhy; hhy=besthhy;
cout << "Iteration "<< iint << " Chi2: " << best_diff << endl; //" Best: "<< best_diff << " RMS: "<< rms<< endl;
#ifdef SAVE_ALL #ifdef SAVE_ALL
debugSaveAll(iint);
for (int ii=0; ii<etabins*etabins; ii++) {
etah[ii]=hhx[ii];
}
sprintf(tit,"/scratch/eta_hhx_%d.tiff",id);
WriteToTiff(etah, tit, etabins, etabins);
for (int ii=0; ii<etabins*etabins; ii++) {
etah[ii]=hhy[ii];
}
sprintf(tit,"/scratch/eta_hhy_%d.tiff",id);
WriteToTiff(etah, tit, etabins, etabins);
for (int ii=0; ii<etabins*etabins; ii++) {
etah[ii]=heta[ii];
}
sprintf(tit,"/scratch/eta_%d.tiff",id);
WriteToTiff(etah, tit, etabins, etabins);
delete [] etah;
#endif #endif
return ; return ;
} }
};
class eta2InterpolationAdaptiveBins : public virtual eta2InterpolationBase, public virtual etaInterpolationAdaptiveBins {
public:
eta2InterpolationAdaptiveBins(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),etaInterpolationAdaptiveBins(nx,ny,ns, nb, emin,emax){
cout << "NSUBPIX is " << ns << " " << nSubPixels << endl;
// cout << "e2pxy " << nb << " " << emin << " " << emax << endl;
};
eta2InterpolationAdaptiveBins(eta2InterpolationAdaptiveBins *orig): etaInterpolationBase(orig), etaInterpolationAdaptiveBins(orig) {};
virtual eta2InterpolationAdaptiveBins* Clone() { return new eta2InterpolationAdaptiveBins(this);};
};
class eta3InterpolationAdaptiveBins : public virtual eta3InterpolationBase, public virtual etaInterpolationAdaptiveBins {
public:
eta3InterpolationAdaptiveBins(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), etaInterpolationAdaptiveBins(nx,ny,ns, nb, emin,emax){
cout << "e3pxy " << nbeta << " " << etamin << " " << etamax << " " << nSubPixels<< endl;
};
eta3InterpolationAdaptiveBins(eta3InterpolationAdaptiveBins *orig): etaInterpolationBase(orig), etaInterpolationAdaptiveBins(orig) {};
virtual eta3InterpolationAdaptiveBins* Clone() { return new eta3InterpolationAdaptiveBins(this);};
}; };
#endif #endif

281
slsDetectorCalibration/interpolations/etaInterpolationBase.h
View File

@ -23,97 +23,53 @@ class etaInterpolationBase : public slsInterpolation {
nbeta=nSubPixels*10; nbeta=nSubPixels*10;
} }
if (etamin>=etamax) { if (etamin>=etamax) {
// cout << "aaa:" <<endl;
etamin=-1; etamin=-1;
etamax=2; etamax=2;
} }
etastep=(etamax-etamin)/nbeta; etastep=(etamax-etamin)/nbeta;
#ifdef MYROOT1
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
heta=new int[nbeta*nbeta]; heta=new int[nbeta*nbeta];
hhx=new float[nbeta*nbeta]; hhx=new float[nbeta*nbeta];
hhy=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];
#endif
//cout << nbeta << " " << etamin << " " << etamax << endl;
}; };
etaInterpolationBase(etaInterpolationBase *orig): slsInterpolation(orig){ etaInterpolationBase(etaInterpolationBase *orig): slsInterpolation(orig){
nbeta=orig->nbeta; nbeta=orig->nbeta;
etamin=orig->etamin; etamin=orig->etamin;
etamax=orig->etamax; etamax=orig->etamax;
rangeMin=orig->rangeMin;
rangeMax=orig->rangeMax;
etastep=(etamax-etamin)/nbeta; etastep=(etamax-etamin)/nbeta;
#ifdef MYROOT1
heta=(TH2D*)(orig->heta)->Clone("heta");
hhx=(TH2D*)(orig->hhx)->Clone("hhx");
hhy=(TH2D*)(orig->hhy)->Clone("hhy");
#endif
#ifndef MYROOT1
heta=new int[nbeta*nbeta]; heta=new int[nbeta*nbeta];
memcpy(heta,orig->heta,nbeta*nbeta*sizeof(int)); memcpy(heta,orig->heta,nbeta*nbeta*sizeof(int));
hhx=new float[nbeta*nbeta]; hhx=new float[nbeta*nbeta];
memcpy(hhx,orig->hhx,nbeta*nbeta*sizeof(float)); memcpy(hhx,orig->hhx,nbeta*nbeta*sizeof(float));
hhy=new float[nbeta*nbeta]; hhy=new float[nbeta*nbeta];
memcpy(hhy,orig->hhy,nbeta*nbeta*sizeof(float)); memcpy(hhy,orig->hhy,nbeta*nbeta*sizeof(float));
hintcorr=new int [nSubPixels*nSubPixels*nPixelsX*nPixelsY];
#endif
}; };
/* virtual etaInterpolationBase* Clone()=0; {
return new etaInterpolationBase(this);
};*/
virtual void resetFlatField() { virtual void resetFlatField() {
#ifndef MYROOT1
for (int ibx=0; ibx<nbeta*nbeta; ibx++) { for (int ibx=0; ibx<nbeta*nbeta; ibx++) {
heta[ibx]=0; heta[ibx]=0;
hhx[ibx]=0; hhx[ibx]=0;
hhy[ibx]=0; hhy[ibx]=0;
} }
#endif
#ifdef MYROOT1
heta->Reset();
hhx->Reset();
hhy->Reset();
#endif
}; };
#ifdef MYROOT1
TH2D *setEta(TH2D *h, int nb=-1, double emin=1, double emax=0)
{
if (h) { heta=h;
nbeta=heta->GetNbinsX();
etamin=heta->GetXaxis()->GetXmin();
etamax=heta->GetXaxis()->GetXmax();
etastep=(etamax-etamin)/nbeta;
}
return heta;
};
TH2D *setFlatField(TH2D *h, int nb=-1, double emin=1, double emax=0)
{
return setEta(h, nb, emin, emax);
};
TH2D *getFlatField(){return setEta(NULL);};
#endif
#ifndef MYROOT1
int *setEta(int *h, int nb=-1, double emin=1, double emax=0) int *setEta(int *h, int nb=-1, double emin=1, double emax=0)
{ {
if (h) { if (h) {
@ -127,6 +83,8 @@ class etaInterpolationBase : public slsInterpolation {
etamin=-1; etamin=-1;
etamax=2; etamax=2;
} }
rangeMin=etamin;
rangeMax=etamax;
etastep=(etamax-etamin)/nbeta; etastep=(etamax-etamin)/nbeta;
} }
return heta; return heta;
@ -143,7 +101,6 @@ class etaInterpolationBase : public slsInterpolation {
int *getFlatField(int &nb, double &emin, double &emax){ int *getFlatField(int &nb, double &emin, double &emax){
nb=nbeta; nb=nbeta;
//cout << "igff* ff has " << nb << " bins " << endl;
emin=etamin; emin=etamin;
emax=etamax; emax=etamax;
return getFlatField(); return getFlatField();
@ -208,28 +165,6 @@ class etaInterpolationBase : public slsInterpolation {
#endif
/* ////////////////////////////////////////////////////////////////////////////// */
#ifdef MYROOT1
TH2D *gethhx()
{
hhx->Scale((double)nSubPixels);
return hhx;
};
TH2D *gethhy()
{
hhy->Scale((double)nSubPixels);
return hhy;
};
#endif
#ifndef MYROOT1
float *gethhx() float *gethhx()
{ {
// hhx->Scale((double)nSubPixels); // hhx->Scale((double)nSubPixels);
@ -241,53 +176,193 @@ float *gethhx()
// hhy->Scale((double)nSubPixels); // hhy->Scale((double)nSubPixels);
return hhy; return hhy;
}; };
#endif
//////////////////////////////////////////////////////////////////////////////
//////////// /*It return position hit for the event in input */ //////////////
/* virtual void getInterpolatedPosition(int x, int y, int *data, double &int_x, double &int_y)=0; */
/* virtual void getInterpolatedPosition(int x, int y, double *data, double &int_x, double &int_y)=0; */
/* virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cl,double &int_x, double &int_y)=0; */
/* virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,int *cl,double &int_x, double &int_y)=0; */
/* virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int corner, double &int_x, double &int_y)=0; */
/* virtual int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay)=0; */
/* virtual int addToFlatField(double totquad,int quad,double *cl,double &etax, double &etay)=0; */
/* virtual int addToFlatField(double *cluster, double &etax, double &etay)=0; */
/* virtual int addToFlatField(int *cluster, double &etax, double &etay)=0; */
virtual int addToFlatField(double etax, double etay){ virtual int addToFlatField(double etax, double etay){
#ifdef MYROOT1
heta->Fill(etax,etay);
#endif
#ifndef MYROOT1
int ex,ey; int ex,ey;
ex=(etax-etamin)/etastep; ex=(etax-etamin)/etastep;
ey=(etay-etamin)/etastep; ey=(etay-etamin)/etastep;
if (ey<nbeta && ex<nbeta && ex>=0 && ey>=0) if (ey<nbeta && ex<nbeta && ex>=0 && ey>=0)
heta[ey*nbeta+ex]++; heta[ey*nbeta+ex]++;
#endif
return 0; return 0;
}; };
// virtual void prepareInterpolation(int &ok)=0; // virtual void prepareInterpolation(int &ok)=0;
void debugSaveAll(int ind=0) {
int ib, 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.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: protected:
#ifdef MYROOT1
TH2D *heta; double calcDiff(double avg, float *hx, float *hy) {
TH2D *hhx; //double p_tot=0;
TH2D *hhy; double diff=0, d;
#endif double bsize=1./nSubPixels;
#ifndef MYROOT1 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);
}
int *heta; int *heta;
float *hhx; float *hhx;
float *hhy; float *hhy;
#endif
int nbeta; int nbeta;
double etamin, etamax, etastep; double etamin, etamax, etastep;
double rangeMin, rangeMax;
double *flat;
int *hintcorr;
}; };

263
slsDetectorCalibration/interpolations/etaInterpolationCleverAdaptiveBins.h Normal file
View File

@ -0,0 +1,263 @@
#ifndef ETA_INTERPOLATION_CLEVER_ADAPTIVEBINS_H
#define ETA_INTERPOLATION_CLEVER_ADAPTIVEBINS_H
#include <cmath>
#include "tiffIO.h"
//#include "etaInterpolationBase.h"
#include "etaInterpolationAdaptiveBins.h"
//#define HSIZE 1
class etaInterpolationCleverAdaptiveBins : public etaInterpolationAdaptiveBins {
private:
// double *gradientX, *gradientY, *gradientXY;
virtual void iterate(float *newhhx, float *newhhy) {
double bsize=1./nSubPixels;
/* double hy[nSubPixels*HSIZE][nbeta]; //profile y */
/* double hx[nSubPixels*HSIZE][nbeta]; //profile x */
// double hix[nSubPixels*HSIZE][nbeta]; //integral of projection x
// double hiy[nSubPixels*HSIZE][nbeta]; //integral of projection y
int ipy, ipx, ippx, ippy;
// double tot_eta_x[nSubPixels*HSIZE];
//double tot_eta_y[nSubPixels*HSIZE];
double mean=0;
double maxflat=0, minflat=0, maxgradX=0, mingradX=0, maxgradY=0, mingradY=0, maxgr=0, mingr=0;
int ix_maxflat, iy_maxflat, ix_minflat, iy_minflat, ix_maxgrX, iy_maxgrX, ix_mingrX, iy_mingrX,ix_maxgrY, iy_maxgrY, ix_mingrY, iy_mingrY, ix_mingr, iy_mingr, ix_maxgr, iy_maxgr;
int maskMin[nSubPixels*nSubPixels], maskMax[nSubPixels*nSubPixels];
//for (int ipy=0; ipy<nSubPixels; ipy++) {
for (ipy=0; ipy<nSubPixels; ipy++) {
for (ipx=0; ipx<nSubPixels; ipx++) {
// cout << ipx << " " << ipy << endl;
mean+=flat[ipx+nSubPixels*ipy]/((double)(nSubPixels*nSubPixels));
}
}
// cout << "Mean is " << mean << endl;
/*** Find local minima and maxima within the staistical uncertainty **/
for (ipy=0; ipy<nSubPixels; ipy++) {
for (ipx=0; ipx<nSubPixels; ipx++) {
if (flat[ipx+nSubPixels*ipy]<mean-3.*sqrt(mean))maskMin[ipx+nSubPixels*ipy]=1; else maskMin[ipx+nSubPixels*ipy]=0;
if (flat[ipx+nSubPixels*ipy]>mean+3.*sqrt(mean)) maskMax[ipx+nSubPixels*ipy]=1; else maskMax[ipx+nSubPixels*ipy]=0;
if (ipx>0 && ipy>0) {
if (flat[ipx+nSubPixels*ipy]<flat[ipx-1+nSubPixels*(ipy-1)]) maskMax[ipx+nSubPixels*ipy]=0;
if (flat[ipx+nSubPixels*ipy]>flat[ipx-1+nSubPixels*(ipy-1)]) maskMin[ipx+nSubPixels*ipy]=0;
}
if (ipx>0 && ipy<nSubPixels-1) {
if (flat[ipx+nSubPixels*ipy]<flat[ipx-1+nSubPixels*(ipy+1)]) maskMax[ipx+nSubPixels*ipy]=0;
if (flat[ipx+nSubPixels*ipy]>flat[ipx-1+nSubPixels*(ipy+1)]) maskMin[ipx+nSubPixels*ipy]=0;
}
if (ipy>0 && ipx<nSubPixels-1) {
if (flat[ipx+nSubPixels*ipy]<flat[ipx+1+nSubPixels*(ipy-1)]) maskMax[ipx+nSubPixels*ipy]=0;
if (flat[ipx+nSubPixels*ipy]>flat[ipx+1+nSubPixels*(ipy-1)]) maskMin[ipx+nSubPixels*ipy]=0;
}
if (ipy<nSubPixels-1 && ipx<nSubPixels-1) {
if (flat[ipx+nSubPixels*ipy]<flat[ipx+1+nSubPixels*(ipy+1)]) maskMax[ipx+nSubPixels*ipy]=0;
if (flat[ipx+nSubPixels*ipy]>flat[ipx+1+nSubPixels*(ipy+1)]) maskMin[ipx+nSubPixels*ipy]=0;
}
if (ipy<nSubPixels-1 ) {
if (flat[ipx+nSubPixels*ipy]<flat[ipx+nSubPixels*(ipy+1)]) maskMax[ipx+nSubPixels*ipy]=0;
if (flat[ipx+nSubPixels*ipy]>flat[ipx+nSubPixels*(ipy+1)]) maskMin[ipx+nSubPixels*ipy]=0;
}
if (ipx<nSubPixels-1) {
if (flat[ipx+nSubPixels*ipy]<flat[ipx+1+nSubPixels*(ipy)]) maskMax[ipx+nSubPixels*ipy]=0;
if (flat[ipx+nSubPixels*ipy]>flat[ipx+1+nSubPixels*(ipy)]) maskMin[ipx+nSubPixels*ipy]=0;
}
if (ipy>0 ) {
if (flat[ipx+nSubPixels*ipy]<flat[ipx+nSubPixels*(ipy-1)]) maskMax[ipx+nSubPixels*ipy]=0;
if (flat[ipx+nSubPixels*ipy]>flat[ipx+nSubPixels*(ipy-1)]) maskMin[ipx+nSubPixels*ipy]=0;
}
if (ipx>0 ) {
if (flat[ipx+nSubPixels*ipy]<flat[ipx-1+nSubPixels*(ipy)]) maskMax[ipx+nSubPixels*ipy]=0;
if (flat[ipx+nSubPixels*ipy]>flat[ipx-1+nSubPixels*(ipy)]) maskMin[ipx+nSubPixels*ipy]=0;
}
// if (maskMin[ipx+nSubPixels*ipy]) cout << ipx << " " << ipy << " is a local minimum " << flat[ipx+nSubPixels*ipy] << endl;
// if (maskMax[ipx+nSubPixels*ipy]) cout << ipx << " " << ipy << " is a local maximum "<< flat[ipx+nSubPixels*ipy] << endl;
}
}
int is_a_border=0;
//initialize the new partition to the previous one
// int ibx_p, iby_p, ibx_n, iby_n;
int ibbx, ibby;
memcpy(newhhx,hhx,nbeta*nbeta*sizeof(float));
memcpy(newhhy,hhy,nbeta*nbeta*sizeof(float));
for (int ibx=0; ibx<nbeta; ibx++) {
for (int iby=0; iby<nbeta; iby++) {
ippy=hhy[ibx+iby*nbeta]*nSubPixels;
ippx=hhx[ibx+iby*nbeta]*nSubPixels;
is_a_border=0;
if (maskMin[ippx+nSubPixels*ippy] || maskMax[ippx+nSubPixels*ippy]) {
for (int ix=-1; ix<2; ix++) {
ibbx=ibx+ix;
if (ibbx<0) ibbx=0;
if (ibbx>nbeta-1) ibbx=nbeta-1;
for (int iy=-1; iy<2; iy++) {
ibby=iby+iy;
if (ibby<0) ibby=0;
if (ibby>nbeta-1) ibby=nbeta-1;
ipy=hhy[ibbx+ibby*nbeta]*nSubPixels;
ipx=hhx[ibbx+ibby*nbeta]*nSubPixels;
if (ipx!=ippx || ipy!=ippy) {
is_a_border=1;
if (maskMin[ippx+nSubPixels*ippy]) {
//increase the region
newhhx[ibbx+ibby*nbeta]=((double)ippx+0.5)/((double)nSubPixels);
newhhy[ibbx+ibby*nbeta]=((double)ippy+0.5)/((double)nSubPixels);
}
if (maskMax[ippx+nSubPixels*ippy]) {
//reduce the region
newhhx[ibx+iby*nbeta]=((double)ipx+0.5)/((double)nSubPixels);
newhhy[ibx+iby*nbeta]=((double)ipy+0.5)/((double)nSubPixels);
}
// cout << ippx << " " << ippy << " " << ibx << " " << iby << " * " << ipx << " " << ipy << " " << ibbx << " " << ibby << endl;
}
}
}
}
}
}
//Check that the resulting histograms are monotonic and they don't have holes!
for (int ibx=0; ibx<nbeta-1; ibx++) {
for (int iby=0; iby<nbeta-1; iby++) {
ippy=newhhy[ibx+iby*nbeta]*nSubPixels;
ippx=newhhx[ibx+iby*nbeta]*nSubPixels;
ipy=newhhy[ibx+(iby+1)*nbeta]*nSubPixels;
ipx=newhhx[ibx+1+iby*nbeta]*nSubPixels;
if ( ippx>ipx)
newhhx[ibx+1+iby*nbeta]=newhhx[ibx+iby*nbeta];
else if (ipx >ippx+1)
newhhx[ibx+1+iby*nbeta]=((double)(ippx+1+0.5))/((double)nSubPixels);
if ( ippy>ipy)
newhhy[ibx+(iby+1)*nbeta]=newhhy[ibx+iby*nbeta];
else if (ipy >ippy+1)
newhhy[ibx+(iby+1)*nbeta]=((double)(ippy+1+0.5))/((double)nSubPixels);
}
}
}
public:
etaInterpolationCleverAdaptiveBins(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationAdaptiveBins(nx,ny,ns, nb, emin,emax){
};
etaInterpolationCleverAdaptiveBins(etaInterpolationCleverAdaptiveBins *orig): etaInterpolationAdaptiveBins(orig){};
virtual etaInterpolationCleverAdaptiveBins* Clone()=0;
/* return new etaInterpolationCleverAdaptiveBins(this); */
/* }; */
};
class eta2InterpolationCleverAdaptiveBins : public virtual eta2InterpolationBase, public virtual etaInterpolationCleverAdaptiveBins {
public:
eta2InterpolationCleverAdaptiveBins(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),etaInterpolationCleverAdaptiveBins(nx,ny,ns, nb, emin,emax){
};
eta2InterpolationCleverAdaptiveBins(eta2InterpolationCleverAdaptiveBins *orig): etaInterpolationBase(orig), etaInterpolationCleverAdaptiveBins(orig) {};
virtual eta2InterpolationCleverAdaptiveBins* Clone() { return new eta2InterpolationCleverAdaptiveBins(this);};
// virtual int *getInterpolatedImage(){return eta2InterpolationBase::getInterpolatedImage();};
/* virtual int *getInterpolatedImage(){ */
/* int ipx, ipy; */
/* cout << "ff" << endl; */
/* calcDiff(1, hhx, hhy); //get flat */
/* double avg=0; */
/* for (ipx=0; ipx<nSubPixels; ipx++) */
/* for (ipy=0; ipy<nSubPixels; ipy++) */
/* avg+=flat[ipx+ipy*nSubPixels]; */
/* avg/=nSubPixels*nSubPixels; */
/* for (int ibx=0 ; ibx<nSubPixels*nPixelsX; ibx++) { */
/* ipx=ibx%nSubPixels-nSubPixels; */
/* if (ipx<0) ipx=nSubPixels+ipx; */
/* for (int iby=0 ; iby<nSubPixels*nPixelsY; iby++) { */
/* ipy=iby%nSubPixels-nSubPixels; */
/* if (ipy<0) ipy=nSubPixels+ipy; */
/* if (flat[ipx+ipy*nSubPixels]>0) */
/* hintcorr[ibx+iby*nSubPixels*nPixelsX]=hint[ibx+iby*nSubPixels*nPixelsX]*(avg/flat[ipx+ipy*nSubPixels]); */
/* else */
/* hintcorr[ibx+iby*nSubPixels*nPixelsX]=hint[ibx+iby*nSubPixels*nPixelsX]; */
/* } */
/* } */
/* return hintcorr; */
/* }; */
};
class eta3InterpolationCleverAdaptiveBins : public virtual eta3InterpolationBase, public virtual etaInterpolationCleverAdaptiveBins {
public:
eta3InterpolationCleverAdaptiveBins(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), etaInterpolationCleverAdaptiveBins(nx,ny,ns, nb, emin,emax){
};
eta3InterpolationCleverAdaptiveBins(eta3InterpolationCleverAdaptiveBins *orig): etaInterpolationBase(orig), etaInterpolationCleverAdaptiveBins(orig) {};
virtual eta3InterpolationCleverAdaptiveBins* Clone() { return new eta3InterpolationCleverAdaptiveBins(this);};
};
#endif

125
slsDetectorCalibration/interpolations/etaInterpolationPosXY.h
View File

@ -51,16 +51,26 @@ class etaInterpolationPosXY : public virtual etaInterpolationBase{
double hix[nbeta]; //integral of projection x double hix[nbeta]; //integral of projection x
double hiy[nbeta]; //integral of projection y double hiy[nbeta]; //integral of projection y
int ii=0; int ii=0;
double etax, etay;
for (int ib=0; ib<nbeta; ib++) { for (int ib=0; ib<nbeta; ib++) {
tot_eta_x=0; tot_eta_x=0;
tot_eta_y=0; tot_eta_y=0;
for (int iby=0; iby<nbeta; iby++) { for (int iby=0; iby<nbeta; iby++) {
hx[iby]=heta[iby+ib*nbeta]; etax=etamin+iby*etastep;
tot_eta_x+=hx[iby]; //cout << etax << endl;
hy[iby]=heta[ib+iby*nbeta]; if (etax>=0 && etax<=1)
tot_eta_y+=hy[iby]; hx[iby]=heta[iby+ib*nbeta];
else {
hx[iby]=0;
}
// tot_eta_x+=hx[iby];
if (etax>=0 && etax<=1)
hy[iby]=heta[ib+iby*nbeta];
else
hy[iby]=0;
// tot_eta_y+=hy[iby];
} }
hix[0]=hx[0]; hix[0]=hx[0];
@ -72,22 +82,17 @@ class etaInterpolationPosXY : public virtual etaInterpolationBase{
} }
ii=0; ii=0;
tot_eta_x=hix[nbeta-1]+1;
tot_eta_y=hiy[nbeta-1]+1;
for (int ibx=0; ibx<nbeta; ibx++) { for (int ibx=0; ibx<nbeta; ibx++) {
if (tot_eta_x==0) { if (tot_eta_x<=0) {
hhx[ibx+ib*nbeta]=((float)ibx)/((float)nbeta); hhx[ibx+ib*nbeta]=-1;
ii=(ibx)/nbeta; //ii=(ibx)/nbeta;
} else //if (hix[ibx]>(ii+1)*tot_eta_x*bsize) } else //if (hix[ibx]>(ii+1)*tot_eta_x*bsize)
{ {
//ii++; //if (hix[ibx]>tot_eta_x*(ii+1)/nSubPixels) ii++;
// cout << ib << " x " << ibx << " " << tot_eta_x << " " << (ii)*tot_eta_x*bsize << " " << ii << endl; hhx[ibx+ib*nbeta]=hix[ibx]/tot_eta_x;
// }
#ifdef MYROOT1
hhx->SetBinContent(ibx+1,ib+1,ii);
#endif
#ifndef MYROOT1
hhx[ibx+ib*nbeta]=hix[ibx]/((float)tot_eta_x);//ii;
#endif
} }
} }
/* if (ii!=(nSubPixels-1)) */ /* if (ii!=(nSubPixels-1)) */
@ -96,53 +101,55 @@ class etaInterpolationPosXY : public virtual etaInterpolationBase{
ii=0; ii=0;
for (int ibx=0; ibx<nbeta; ibx++) { for (int ibx=0; ibx<nbeta; ibx++) {
if (tot_eta_y==0) { if (tot_eta_y<=0) {
hhx[ibx+ib*nbeta]=((float)ibx)/((float)nbeta); hhy[ib+ibx*nbeta]=-1;
ii=(ibx*nSubPixels)/nbeta; //ii=(ibx*nSubPixels)/nbeta;
} else //if (hiy[ibx]>(ii+1)*tot_eta_y*bsize) } else {
{ //if (hiy[ibx]>tot_eta_y*(ii+1)/nSubPixels) ii++;
//ii++; hhy[ib+ibx*nbeta]=hiy[ibx]/tot_eta_y;
//cout << ib << " y " << ibx << " " << tot_eta_y << " "<< (ii)*tot_eta_y*bsize << " " << ii << endl; }
//}
#ifdef MYROOT1
hhy->SetBinContent(ib+1,ibx+1,ii);
#endif
#ifndef MYROOT1
hhy[ib+ibx*nbeta]=hiy[ibx]/((float)tot_eta_y);//ii;
#endif
}
} }
/* if (ii!=(nSubPixels-1)) */
/* cout << ib << " y " << tot_eta_y << " " << (ii+1)*tot_eta_y*bsize << " " << ii << " " << hiy[nbeta-1]<< endl; */
// cout << "y " << nbeta << " " << (ii+1)*tot_eta_x*bsize << " " << ii << endl;
} }
int ibx, iby, ib;
iby=0;
while (hhx[iby*nbeta+nbeta/2]<0) iby++;
for (ib=0; ib<iby;ib++) {
for (ibx=0; ibx<nbeta;ibx++)
hhx[ibx+nbeta*ib]=hhx[ibx+nbeta*iby];
}
iby=nbeta-1;
while (hhx[iby*nbeta+nbeta/2]<0) iby--;
for (ib=iby+1; ib<nbeta;ib++) {
for (ibx=0; ibx<nbeta;ibx++)
hhx[ibx+nbeta*ib]=hhx[ibx+nbeta*iby];
}
iby=0;
while (hhy[nbeta/2*nbeta+iby]<0) iby++;
for (ib=0; ib<iby;ib++) {
for (ibx=0; ibx<nbeta;ibx++)
hhy[ib+nbeta*ibx]=hhy[iby+nbeta*ibx];
}
iby=nbeta-1;
while (hhy[nbeta/2*nbeta+iby]<0) iby--;
for (ib=iby+1; ib<nbeta;ib++) {
for (ibx=0; ibx<nbeta;ibx++)
hhy[ib+nbeta*ibx]=hhy[iby+nbeta*ibx];
}
#ifdef SAVE_ALL #ifdef SAVE_ALL
char tit[10000]; debugSaveAll();
float *etah=new float[nbeta*nbeta];
int etabins=nbeta;
for (int ii=0; ii<etabins*etabins; ii++) {
etah[ii]=hhx[ii];
}
sprintf(tit,"/scratch/eta_hhx_%d.tiff",id);
WriteToTiff(etah, tit, etabins, etabins);
for (int ii=0; ii<etabins*etabins; ii++) {
etah[ii]=hhy[ii];
}
sprintf(tit,"/scratch/eta_hhy_%d.tiff",id);
WriteToTiff(etah, tit, etabins, etabins);
for (int ii=0; ii<etabins*etabins; ii++) {
etah[ii]=heta[ii];
}
sprintf(tit,"/scratch/eta_%d.tiff",id);
WriteToTiff(etah, tit, etabins, etabins);
delete [] etah;
#endif #endif
return ; return ;
} }

233
slsDetectorCalibration/interpolations/slsInterpolation.h
View File

@ -1,12 +1,6 @@
#ifndef SLS_INTERPOLATION_H #ifndef SLS_INTERPOLATION_H
#define SLS_INTERPOLATION_H #define SLS_INTERPOLATION_H
#ifdef MYROOT1
#include <TObject.h>
#include <TTree.h>
#include <TH2F.h>
#endif
#include <cstdlib> #include <cstdlib>
#ifndef MY_TIFF_IO_H #ifndef MY_TIFF_IO_H
#include "tiffIO.h" #include "tiffIO.h"
@ -37,13 +31,7 @@ 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) :nPixelsX(nx), nPixelsY(ny), nSubPixels(ns), id(0) {
#ifdef MYROOT1
hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
#endif
#ifndef MYROOT1
hint=new int[ns*nx*ns*ny]; hint=new int[ns*nx*ns*ny];
#endif
}; };
@ -51,14 +39,9 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
nPixelsX=orig->nPixelsX; nPixelsX=orig->nPixelsX;
nPixelsY=orig->nPixelsY; nPixelsY=orig->nPixelsY;
nSubPixels=orig->nSubPixels; nSubPixels=orig->nSubPixels;
#ifdef MYROOT1
hint=(TH2F*)(orig->hint)->Clone("hint");
#endif
#ifndef MYROOT1
hint=new int[nSubPixels*nPixelsX*nSubPixels*nPixelsY]; hint=new int[nSubPixels*nPixelsX*nSubPixels*nPixelsY];
memcpy(hint, orig->hint,nSubPixels*nPixelsX*nSubPixels*nPixelsY*sizeof(int)); memcpy(hint, orig->hint,nSubPixels*nPixelsX*nSubPixels*nPixelsY*sizeof(int));
#endif
}; };
@ -94,11 +77,6 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
//create interpolated image //create interpolated image
//returns interpolated image //returns interpolated image
#ifdef MYROOT1
virtual TH2F *getInterpolatedImage(){return hint;};
#endif
#ifndef MYROOT1
virtual int *getInterpolatedImage(){ virtual int *getInterpolatedImage(){
// cout << "return interpolated image " << endl; // cout << "return interpolated image " << endl;
/* for (int i=0; i<nSubPixels* nSubPixels* nPixelsX*nPixelsY; i++) { */ /* for (int i=0; i<nSubPixels* nSubPixels* nPixelsX*nPixelsY; i++) { */
@ -106,7 +84,6 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
/* } */ /* } */
return hint; return hint;
}; };
#endif
@ -114,11 +91,12 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
void *writeInterpolatedImage(const char * imgname) { void *writeInterpolatedImage(const char * imgname) {
//cout << "!" <<endl; //cout << "!" <<endl;
float *gm=NULL; float *gm=NULL;
int *dummy=getInterpolatedImage();
gm=new float[ nSubPixels* nSubPixels* nPixelsX*nPixelsY]; gm=new float[ nSubPixels* nSubPixels* nPixelsX*nPixelsY];
if (gm) { if (gm) {
for (int ix=0; ix<nPixelsX*nSubPixels; ix++) { for (int ix=0; ix<nPixelsX*nSubPixels; ix++) {
for (int iy=0; iy<nPixelsY*nSubPixels; iy++) { for (int iy=0; iy<nPixelsY*nSubPixels; iy++) {
gm[iy*nPixelsX*nSubPixels+ix]=hint[iy*nPixelsX*nSubPixels+ix]; gm[iy*nPixelsX*nSubPixels+ix]=dummy[iy*nPixelsX*nSubPixels+ix];
} }
} }
WriteToTiff(gm, imgname,nSubPixels* nPixelsX ,nSubPixels* nPixelsY); WriteToTiff(gm, imgname,nSubPixels* nPixelsX ,nSubPixels* nPixelsY);
@ -142,16 +120,11 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
virtual void clearInterpolatedImage() { virtual void clearInterpolatedImage() {
#ifdef MYROOT1
hint->Reset();
#endif
#ifndef MYROOT1
for (int ix=0; ix<nPixelsX*nSubPixels; ix++) { for (int ix=0; ix<nPixelsX*nSubPixels; ix++) {
for (int iy=0; iy<nPixelsY*nSubPixels; iy++) { for (int iy=0; iy<nPixelsY*nSubPixels; iy++) {
hint[iy*nPixelsX*nSubPixels+ix]=0; hint[iy*nPixelsX*nSubPixels+ix]=0;
} }
} }
#endif
}; };
@ -159,11 +132,6 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
#ifdef MYROOT1
TH2F *addToImage(double int_x, double int_y){hint->Fill(int_x, int_y); return hint;};
#endif
#ifndef MYROOT1
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)nSubPixels)*int_y;
int ix=((double)nSubPixels)*int_x; int ix=((double)nSubPixels)*int_x;
@ -176,7 +144,6 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
return hint; return hint;
}; };
#endif
virtual int addToFlatField(double *cluster, double &etax, double &etay)=0; virtual int addToFlatField(double *cluster, double &etax, double &etay)=0;
@ -185,19 +152,11 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
virtual int addToFlatField(double totquad,int quad,double *cluster,double &etax, double &etay)=0; virtual int addToFlatField(double totquad,int quad,double *cluster,double &etax, double &etay)=0;
virtual int addToFlatField(double etax, double etay)=0; virtual int addToFlatField(double etax, double etay)=0;
#ifdef MYROOT1
virtual TH2D *getFlatField(){return NULL;};
virtual TH2D *setFlatField(TH2D *h, int nb=-1, double emin=-1, double emax=-1){return NULL;};
virtual TH2D *getFlatField(int &nb, double &emin, double &emax){nb=0; emin=0; emax=0; return getFlatField();};
#endif
#ifndef MYROOT1
virtual int *getFlatField(){return NULL;}; virtual int *getFlatField(){return NULL;};
virtual int *setFlatField(int *h, int nb=-1, double emin=-1, double emax=-1){return NULL;}; virtual int *setFlatField(int *h, int nb=-1, double emin=-1, double emax=-1){return NULL;};
virtual void *writeFlatField(const char * imgname){return NULL;}; virtual void *writeFlatField(const char * imgname){return NULL;};
virtual void *readFlatField(const char * imgname, int nb=-1, double emin=1, double emax=0){return NULL;}; virtual void *readFlatField(const char * imgname, int nb=-1, double emin=1, double emax=0){return NULL;};
virtual int *getFlatField(int &nb, double &emin, double &emax){nb=0; emin=0; emax=0; return getFlatField();}; virtual int *getFlatField(int &nb, double &emin, double &emax){nb=0; emin=0; emax=0; return getFlatField();};
#endif
virtual void resetFlatField()=0; virtual void resetFlatField()=0;
@ -215,10 +174,10 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
static int calcQuad(double *cl, double &sum, double &totquad, double sDum[2][2]){ static int calcQuad(double *cl, double &sum, double &totquad, double sDum[2][2]){
int corner = UNDEFINED_QUADRANT; int corner = UNDEFINED_QUADRANT;
double *cluster[3]; /* double *cluster[3]; */
cluster[0]=cl; /* cluster[0]=cl; */
cluster[1]=cl+3; /* cluster[1]=cl+3; */
cluster[2]=cl+6; /* cluster[2]=cl+6; */
sum=0; sum=0;
double sumBL=0; double sumBL=0;
@ -228,11 +187,11 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
int xoff=0, yoff=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+=cluster[iy][ix]; sum+=cl[ix+3*iy];
if (ix<=1 && iy<=1) sumBL+=cluster[iy][ix]; if (ix<=1 && iy<=1) sumBL+=cl[ix+iy*3];
if (ix<=1 && iy>=1) sumTL+=cluster[iy][ix]; if (ix<=1 && iy>=1) sumTL+=cl[ix+iy*3];
if (ix>=1 && iy<=1) sumBR+=cluster[iy][ix]; if (ix>=1 && iy<=1) sumBR+=cl[ix+iy*3];
if (ix>=1 && iy>=1) sumTR+=cluster[iy][ix]; if (ix>=1 && iy>=1) sumTR+=cl[ix+iy*3];
} }
} }
@ -240,6 +199,8 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
/* sDum[0][1] = cluster[0][1]; sDum[1][1] = cluster[1][1]; */ /* sDum[0][1] = cluster[0][1]; sDum[1][1] = cluster[1][1]; */
corner = BOTTOM_LEFT; corner = BOTTOM_LEFT;
totquad=sumBL; totquad=sumBL;
xoff=0;
yoff=0;
if(sumTL >= totquad){ if(sumTL >= totquad){
@ -274,7 +235,7 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
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++) {
sDum[iy][ix] = cluster[iy+yoff][ix+xoff]; sDum[iy][ix] = cl[ix+xoff+(iy+yoff)*3];
} }
} }
@ -396,7 +357,7 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
// int quad; // int quad;
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++) {
val=cl[iy+3*ix]; val=cl[ix+3*iy];
sum+=val; sum+=val;
if (iy==0) l+=val; if (iy==0) l+=val;
if (iy==2) r+=val; if (iy==2) r+=val;
@ -440,92 +401,92 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
} }
static int calcMyEta(double totquad, int quad, double *cl, double &etax, double &etay) { /* static int calcMyEta(double totquad, int quad, double *cl, double &etax, double &etay) { */
double l,r,t,b, sum; /* double l,r,t,b, sum; */
int yoff; /* int yoff; */
switch (quad) { /* switch (quad) { */
case BOTTOM_LEFT: /* case BOTTOM_LEFT: */
case BOTTOM_RIGHT: /* case BOTTOM_RIGHT: */
yoff=0; /* yoff=0; */
break; /* break; */
case TOP_LEFT: /* case TOP_LEFT: */
case TOP_RIGHT: /* case TOP_RIGHT: */
yoff=1; /* yoff=1; */
break; /* break; */
default: /* default: */
; /* ; */
} /* } */
l=cl[0+yoff*3]+cl[0+yoff*3+3]; /* l=cl[0+yoff*3]+cl[0+yoff*3+3]; */
r=cl[2+yoff*3]+cl[2+yoff*3+3]; /* r=cl[2+yoff*3]+cl[2+yoff*3+3]; */
b=cl[0+yoff*3]+cl[1+yoff*3]*cl[2+yoff*3]; /* b=cl[0+yoff*3]+cl[1+yoff*3]*cl[2+yoff*3]; */
t=cl[0+yoff*3+3]+cl[1+yoff*3+3]*cl[0+yoff*3+3]; /* t=cl[0+yoff*3+3]+cl[1+yoff*3+3]*cl[0+yoff*3+3]; */
sum=t+b; /* sum=t+b; */
if (sum>0) { /* if (sum>0) { */
etax=(-l+r)/sum; /* etax=(-l+r)/sum; */
etay=(+t)/sum; /* etay=(+t)/sum; */
} /* } */
return -1; /* return -1; */
} /* } */
static int calcMyEta(double totquad, int quad, int *cl, double &etax, double &etay) { /* static int calcMyEta(double totquad, int quad, int *cl, double &etax, double &etay) { */
double l,r,t,b, sum; /* double l,r,t,b, sum; */
int yoff; /* int yoff; */
switch (quad) { /* switch (quad) { */
case BOTTOM_LEFT: /* case BOTTOM_LEFT: */
case BOTTOM_RIGHT: /* case BOTTOM_RIGHT: */
yoff=0; /* yoff=0; */
break; /* break; */
case TOP_LEFT: /* case TOP_LEFT: */
case TOP_RIGHT: /* case TOP_RIGHT: */
yoff=1; /* yoff=1; */
break; /* break; */
default: /* default: */
; /* ; */
} /* } */
l=cl[0+yoff*3]+cl[0+yoff*3+3]; /* l=cl[0+yoff*3]+cl[0+yoff*3+3]; */
r=cl[2+yoff*3]+cl[2+yoff*3+3]; /* r=cl[2+yoff*3]+cl[2+yoff*3+3]; */
b=cl[0+yoff*3]+cl[1+yoff*3]*cl[2+yoff*3]; /* b=cl[0+yoff*3]+cl[1+yoff*3]*cl[2+yoff*3]; */
t=cl[0+yoff*3+3]+cl[1+yoff*3+3]*cl[0+yoff*3+3]; /* t=cl[0+yoff*3+3]+cl[1+yoff*3+3]*cl[0+yoff*3+3]; */
sum=t+b; /* sum=t+b; */
if (sum>0) { /* if (sum>0) { */
etax=(-l+r)/sum; /* etax=(-l+r)/sum; */
etay=(+t)/sum; /* etay=(+t)/sum; */
} /* } */
return -1; /* return -1; */
} /* } */
static int calcEta3X(double *cl, double &etax, double &etay, double &sum) { /* static int calcEta3X(double *cl, double &etax, double &etay, double &sum) { */
double l,r,t,b; /* double l,r,t,b; */
sum=cl[0]+cl[1]+cl[2]+cl[3]+cl[4]+cl[5]+cl[6]+cl[7]+cl[8]; /* sum=cl[0]+cl[1]+cl[2]+cl[3]+cl[4]+cl[5]+cl[6]+cl[7]+cl[8]; */
if (sum>0) { /* if (sum>0) { */
l=cl[3]; /* l=cl[3]; */
r=cl[5]; /* r=cl[5]; */
b=cl[1]; /* b=cl[1]; */
t=cl[7]; /* t=cl[7]; */
etax=(-l+r)/sum; /* etax=(-l+r)/sum; */
etay=(-b+t)/sum; /* etay=(-b+t)/sum; */
} /* } */
return -1; /* return -1; */
} /* } */
static int calcEta3X(int *cl, double &etax, double &etay, double &sum) { /* static int calcEta3X(int *cl, double &etax, double &etay, double &sum) { */
double l,r,t,b; /* double l,r,t,b; */
sum=cl[0]+cl[1]+cl[2]+cl[3]+cl[4]+cl[5]+cl[6]+cl[7]+cl[8]; /* sum=cl[0]+cl[1]+cl[2]+cl[3]+cl[4]+cl[5]+cl[6]+cl[7]+cl[8]; */
if (sum>0) { /* if (sum>0) { */
l=cl[3]; /* l=cl[3]; */
r=cl[5]; /* r=cl[5]; */
b=cl[1]; /* b=cl[1]; */
t=cl[7]; /* t=cl[7]; */
etax=(-l+r)/sum; /* etax=(-l+r)/sum; */
etay=(-b+t)/sum; /* etay=(-b+t)/sum; */
} /* } */
return -1; /* return -1; */
} /* } */
@ -535,14 +496,8 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
protected: protected:
int nPixelsX, nPixelsY; int nPixelsX, nPixelsY;
int nSubPixels; int nSubPixels;
#ifdef MYROOT1
TH2F *hint;
#endif
#ifndef MYROOT1
int *hint;
#endif
int id; int id;
int *hint;
}; };
#endif #endif

41
slsDetectorCalibration/moenchExecutables/moench03Interpolation.cpp
View File

@ -14,14 +14,14 @@
#include "single_photon_hit.h" #include "single_photon_hit.h"
#endif #endif
#include "etaInterpolationPosXY.h" //#include "etaInterpolationPosXY.h"
#include "noInterpolation.h" #include "noInterpolation.h"
//#include "etaInterpolationAdaptiveBins.h" #include "etaInterpolationCleverAdaptiveBins.h"
//#include "etaInterpolationRandomBins.h" //#include "etaInterpolationRandomBins.h"
using namespace std; using namespace std;
#define NC 400 #define NC 400
#define NR 400 #define NR 400
#define MAX_ITERATIONS (nSubPixels*100)
#define XTALK #define XTALK
@ -65,9 +65,10 @@ int main(int argc, char *argv[]) {
float cmax=atof(argv[iarg++]); float cmax=atof(argv[iarg++]);
cout << "Energy min: " << cmin << endl; cout << "Energy min: " << cmin << endl;
cout << "Energy max: " << cmax << endl; cout << "Energy max: " << cmax << endl;
//int etabins=500;
int etabins=1000;//nsubpix*2*100; int etabins=1000;//nsubpix*2*100;
double etamin=-1, etamax=2; double etamin=-1, etamax=2;
//double etamin=-0.1, etamax=1.1;
double eta3min=-2, eta3max=2; double eta3min=-2, eta3max=2;
int quad; int quad;
double sum, totquad; double sum, totquad;
@ -90,15 +91,18 @@ int main(int argc, char *argv[]) {
single_photon_hit cl(3,3); single_photon_hit cl(3,3);
#endif #endif
int nSubPixels=nsubpix;
//eta2InterpolationPosXY *interp=new eta2InterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax);
eta2InterpolationCleverAdaptiveBins *interp=new eta2InterpolationCleverAdaptiveBins(NC, NR, nsubpix, etabins, etamin, etamax);
eta2InterpolationPosXY *interp=new eta2InterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax);
#ifndef FF #ifndef FF
cout << "read ff " << argv[2] << endl; cout << "read ff " << argv[2] << endl;
sprintf(fname,"%s",argv[2]); sprintf(fname,"%s",argv[2]);
interp->readFlatField(fname); interp->readFlatField(fname);
interp->prepareInterpolation(ok); interp->prepareInterpolation(ok, MAX_ITERATIONS);
// return 0;
#endif #endif
#ifdef FF #ifdef FF
cout << "Will write eta file " << argv[2] << endl; cout << "Will write eta file " << argv[2] << endl;
@ -147,22 +151,27 @@ int main(int argc, char *argv[]) {
nph++; nph++;
// if (sum>200 && sum<580) { // if (sum>200 && sum<580) {
// interp->getInterpolatedPosition(cl.x,cl.y, totquad,quad,cl.get_cluster(),int_x, int_y); // interp->getInterpolatedPosition(cl.x,cl.y, totquad,quad,cl.get_cluster(),int_x, int_y);
#ifdef SOLEIL // #ifdef SOLEIL
if (cl.x>210 && cl.x<240 && cl.y>210 && cl.y<240) { // if (cl.x>210 && cl.x<240 && cl.y>210 && cl.y<240) {
#endif // #endif
#ifndef FF #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->addToImage(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);
#endif #endif
#ifdef FF #ifdef FF
interp->addToFlatField(cl.get_cluster(), etax, etay); interp->addToFlatField(cl.get_cluster(), etax, etay);
#endif #endif
#ifdef SOLEIL // #ifdef SOLEIL
} // }
#endif // #endif
if (nph%1000000==0) cout << nph << endl; if (nph%1000000==0) cout << nph << endl;
if (nph%100000000==0) { if (nph%10000000==0) {
#ifndef FF #ifndef FF
interp->writeInterpolatedImage(outfname); interp->writeInterpolatedImage(outfname);
#endif #endif

15
slsDetectorCalibration/single_photon_hit.h
View File

@ -182,6 +182,21 @@ class single_photon_hit {
return 0; return 0;
}; };
void print() {
int ix, iy;
for (int iy=0; iy<dy; iy++) {
for (int ix=0; ix<dx; ix++) {
printf("%d \t",data[ix+iy*dx]);
}
printf("\n");
}
}
/** /**
assign the value to the element of the cluster matrix, with relative coordinates where the center of the cluster is (0,0) assign the value to the element of the cluster matrix, with relative coordinates where the center of the cluster is (0,0)
\param v value to be set \param v value to be set

2
slsDetectorCalibration/tiffIO.cpp
View File

@ -18,7 +18,7 @@ void *WriteToTiff(float * imgData, const char * imgname, int nrow, int ncol){
int sampleperpixel=1; int sampleperpixel=1;
// unsigned char * buff=NULL; // unsigned char * buff=NULL;
tsize_t linebytes; tsize_t linebytes;
cout << "--" <<endl; // cout << "--" <<endl;
TIFF * tif = TIFFOpen(imgname,"w"); TIFF * tif = TIFFOpen(imgname,"w");
if (tif) { if (tif) {
TIFFSetField(tif,TIFFTAG_IMAGEWIDTH,ncol); TIFFSetField(tif,TIFFTAG_IMAGEWIDTH,ncol);