detectors/slsDetectorPackage
25
0
Fork 0
mirror of https://github.com/slsdetectorgroup/slsDetectorPackage.git synced 2025-06-23 18:17:59 +02:00
Code Issues Actions Packages Releases Activity

some fixes in the cluster finder

Browse Source
This commit is contained in:
bergamaschi
2018-01-30 10:24:56 +01:00
parent f264691dc1
commit d1ef6bf94b
10 changed files with 323 additions and 212 deletions
Download Patch File Download Diff File Expand all files Collapse all files

203
slsDetectorCalibration/interpolations/etaInterpolationAdaptiveBins.h
View File

@ -3,17 +3,19 @@
#include "tiffIO.h"
#include "etaInterpolationBase.h"
//#include "etaInterpolationBase.h"
#include "etaInterpolationPosXY.h"
class etaInterpolationAdaptiveBins : public etaInterpolationPosXY {
class etaInterpolationAdaptiveBins : public etaInterpolationBase{
private:
double calcDiff(double avg, double *hx, double *hy) {
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++) {
for (int ipy=0; ipy=nSubPixels; ipy++) {
for (int ipx=0; ipx<nSubPixels; ipx++) {
for (int ipy=0; ipy<nSubPixels; ipy++) {
p_tot=0;
for (int ibx=0; ibx<nbeta; ibx++) {
for (int iby=0; iby<nbeta; iby++) {
@ -23,15 +25,19 @@ class etaInterpolationAdaptiveBins : public etaInterpolationBase{
}
}
cout << p_tot << " \t ";
diff+=(p_tot-avg)*(p_tot-avg);
}
cout << "\n";
}
return diff;
}
void iterate(double *newhhx, double *newhhy) {
void iterate(float *newhhx, float *newhhy) {
double bsize=1./nSubPixels;
@ -42,55 +48,71 @@ class etaInterpolationAdaptiveBins : public etaInterpolationBase{
double tot_eta_x=0;
double tot_eta_y=0;
for (int ipy=0; ipy=nSubPixels; ipy++) {
for (int ipy=0; ipy<nSubPixels; ipy++) {
for (int ibx=0; ibx<nbeta; ibx++) {
for (int iby=0; iby<nbeta; iby++) {
if ( hhy[ibx+iby*nbeta]>=((ipy)*bsize) && hhy[ibx+iby*nbeta]<((ipy+1)*bsize)) {
hx[ibx]+=heta[ibx+iby*nbeta];
}
if ( hhx[ibx+iby*nbeta]>=((ipy)*bsize) && hhx[ibx+iby*nbeta]<((ipy+1)*bsize)) {
hy[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];
for (int ipy=0; ipy=nSubPixels; ipy++) {
for (int ibx=0; ibx<nbeta; ibx++) {
for (int iby=0; iby<nbeta; iby++) {
if ( hhy[ibx+iby*nbeta]>=((ipy)*bsize) && hhy[ibx+iby*nbeta]<((ipy+1)*bsize)) {
newhhx[ibx+iby*nbeta]=hix[ibx]/((double)tot_eta_x);
}
if ( hhx[ibx+iby*nbeta]>=((ipy)*bsize) && hhx[ibx+iby*nbeta]<((ipy+1)*bsize)) {
newhhy[ibx+iby*nbeta]=hiy[iby]/((double)tot_eta_y);
}
}
hx[ibx]=0;
hy[ibx]=0;
}
tot_eta_x=0;
tot_eta_y=0;
// cout << ipy << " " << ((ipy)*bsize) << " " << ((ipy+1)*bsize) << endl;
for (int ibx=0; ibx<nbeta; ibx++) {
for (int iby=0; iby<nbeta; iby++) {
if (hhy[ibx+iby*nbeta]>=((ipy)*bsize) && hhy[ibx+iby*nbeta]<=((ipy+1)*bsize)) {
hx[ibx]+=heta[ibx+iby*nbeta];
tot_eta_x+=heta[ibx+iby*nbeta];
}
if (hhx[ibx+iby*nbeta]>=((ipy)*bsize) && hhx[ibx+iby*nbeta]<=((ipy+1)*bsize)) {
hy[iby]+=heta[ibx+iby*nbeta];
tot_eta_y+=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 (int ibx=0; ibx<nbeta; ibx++) {
for (int iby=0; iby<nbeta; iby++) {
if ( hhy[ibx+iby*nbeta]>=((ipy)*bsize) && hhy[ibx+iby*nbeta]<=((ipy+1)*bsize)) {
newhhx[ibx+iby*nbeta]=hix[ibx]/((double)tot_eta_x);
if (newhhx[ibx+iby*nbeta]>1) 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 (hhx[ibx+iby*nbeta]>=((ipy)*bsize) && hhx[ibx+iby*nbeta]<=((ipy+1)*bsize)) {
newhhy[ibx+iby*nbeta]=hiy[iby]/((double)tot_eta_y);
if (newhhy[ibx+iby*nbeta]>1) cout << "***"<< ibx << " " << iby << newhhy[ibx+iby*nbeta] << endl;
// if (ipy==3 && iby==10) cout << newhhy[ibx+iby*nbeta] << " " << hiy[iby] << " " << ibx+iby*nbeta << endl;
}
}
}
}
}
public:
etaInterpolationAdaptiveBins(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){};
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(etaInterpolationAdaptiveBins *orig): etaInterpolationBase(orig){};
etaInterpolationAdaptiveBins(etaInterpolationAdaptiveBins *orig): etaInterpolationPosXY(orig){};
virtual etaInterpolationAdaptiveBins* Clone() {
@ -103,7 +125,7 @@ class etaInterpolationAdaptiveBins : public etaInterpolationBase{
virtual void prepareInterpolation(int &ok)
{
ok=1;
cout << "Adaptive bins" << endl;
///*Eta Distribution Rebinning*///
double bsize=1./nSubPixels; //precision
@ -113,7 +135,6 @@ class etaInterpolationAdaptiveBins : public etaInterpolationBase{
double tot_eta_y=0;
for (int ip=0; ip<nbeta*nbeta; ip++)
tot_eta+=heta[ip];
cout << "total eta entries is :"<< tot_eta << endl;
if (tot_eta<=0) {ok=0; return;};
@ -125,41 +146,82 @@ class etaInterpolationAdaptiveBins : public etaInterpolationBase{
/** initialize distribution to linear interpolation */
for (int ibx=0; ibx<nbeta; ibx++) {
for (int ib=0; ib<nbeta; ib++) {
hhx[ibx+ib*nbeta]=((float)ibx)/((float)nbeta);
hhy[ibx+ib*nbeta]=((float)ibx)/((float)nbeta);
}
}
// for (int ibx=0; ibx<nbeta; ibx++) {
// for (int ib=0; ib<nbeta; ib++) {
// hhx[ibx+ib*nbeta]=((float)ibx)/((float)nbeta);
// hhy[ibx+ib*nbeta]=((float)ib)/((float)nbeta);
// }
// }
int nint=10;
etaInterpolationPosXY::prepareInterpolation(ok);
double thr=1./((double)/nSubPixels);
double avg=tot_eta/((double)(nSubPixels*nSubPixels));
#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 avg=tot_eta/((double)(nSubPixels*nSubPixels));
cout << "total eta entries is :"<< tot_eta << " avg: "<< avg << endl;
cout << "Start " << endl;
double old_diff=calcDiff(avg, hhx, hhy), new_diff=old_diff+1;
cout << " diff= " << new_diff << endl;
int iint=0;
double newhhx=new double[nbeta*nbeta]; //profile x
double newhhy=new double[nbeta*nbeta]; //profile y
float *newhhx=new float[nbeta*nbeta]; //profile x
float *newhhy=new float[nbeta*nbeta]; //profile y
while (iint<nint) {
cout << "Iteration " << iint << endl;
iterate(newhhx,newhhy);
new_diff=calcDiff(avg, newhhx, newhhy);
if (new_diff<old_diff) {
cout << " diff= " << new_diff << 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<old_diff) {
delete [] hhx;
delete [] hhy;
hhx=newhhx;
hhy=newhhy;
newhhx=new double[nbeta*nbeta]; //profile x
newhhy=new double[nbeta*nbeta]; //profile y
newhhx=new float[nbeta*nbeta]; //profile x
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;
}
// } /* else { */
/* cout << "Difference not decreasing after "<< iint << " iterations (" << old_diff << " < " << new_diff << ")"<< endl; */
/* break; */
/* } */
iint++;
}
@ -170,10 +232,7 @@ class etaInterpolationAdaptiveBins : public etaInterpolationBase{
#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];
@ -193,7 +252,7 @@ class etaInterpolationAdaptiveBins : public etaInterpolationBase{
sprintf(tit,"/scratch/eta_%d.tiff",id);
WriteToTiff(etah, tit, etabins, etabins);
delete [] etah;
#endif
#endif
return ;
}