formatted slsDetectorCalibration

This commit is contained in:
Erik Frojdh
2022-01-24 11:26:56 +01:00
parent 623b1de8a0
commit c554bbb2d3
77 changed files with 15998 additions and 15890 deletions

View File

@ -3,285 +3,316 @@
#ifndef ETA_INTERPOLATION_ADAPTIVEBINS_H
#define ETA_INTERPOLATION_ADAPTIVEBINS_H
#include <cmath>
#include "tiffIO.h"
#include <cmath>
//#include "etaInterpolationBase.h"
#include "etaInterpolationPosXY.h"
class etaInterpolationAdaptiveBins : public etaInterpolationPosXY {
// protected:
// protected:
private:
virtual void iterate(float *newhhx, float *newhhy) {
private:
double bsize = 1. / nSubPixels;
virtual void iterate(float *newhhx, float *newhhy) {
double hy[nSubPixels][nbeta]; // profile y
double hx[nSubPixels][nbeta]; // profile x
double hix[nSubPixels][nbeta]; // integral of projection x
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 bsize=1./nSubPixels;
double hy[nSubPixels][nbeta]; //profile y
double hx[nSubPixels][nbeta]; //profile x
double hix[nSubPixels][nbeta]; //integral of projection x
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++) {
for (ipy=0; ipy<nSubPixels; ipy++) {
for (int ibx=0; ibx<nbeta; ibx++) {
hx[ipy][ibx]=0;
hy[ipy][ibx]=0;
}
}
for (ipy = 0; ipy < nSubPixels; ipy++) {
for (int ibx = 0; ibx < nbeta; ibx++) {
hx[ipy][ibx] = 0;
hy[ipy][ibx] = 0;
}
}
// cout << ipy << " " << ((ipy)*bsize) << " " << ((ipy+1)*bsize) << endl;
for (int ibx=0; ibx<nbeta; ibx++) {
for (int iby=0; iby<nbeta; iby++) {
ipy=hhy[ibx+iby*nbeta]*nSubPixels;
if (ipy<0) ipy=0;
if (ipy>=nSubPixels) ipy=nSubPixels-1;
hx[ipy][ibx]+=heta[ibx+iby*nbeta];
// cout << ipy << " " << ((ipy)*bsize) << " " << ((ipy+1)*bsize) <<
// endl;
for (int ibx = 0; ibx < nbeta; ibx++) {
for (int iby = 0; iby < nbeta; iby++) {
ipy = hhy[ibx + iby * nbeta] * nSubPixels;
if (ipy < 0)
ipy = 0;
if (ipy >= nSubPixels)
ipy = nSubPixels - 1;
hx[ipy][ibx] += heta[ibx + iby * nbeta];
ipx = hhx[ibx + iby * nbeta] * nSubPixels;
if (ipx < 0)
ipx = 0;
if (ipx >= nSubPixels)
ipx = nSubPixels - 1;
hy[ipx][iby] += heta[ibx + iby * nbeta];
}
}
ipx=hhx[ibx+iby*nbeta]*nSubPixels;
if (ipx<0) ipx=0;
if (ipx>=nSubPixels) ipx=nSubPixels-1;
hy[ipx][iby]+=heta[ibx+iby*nbeta];
}
}
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 iby=0; iby<nbeta; iby++) {
// if ( hhy[ibx+iby*nbeta]>=((ipy)*bsize) && hhy[ibx+iby*nbeta]<=((ipy+1)*bsize)) {
ipy=hhy[ibx+iby*nbeta]*nSubPixels;
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++) {
if (ipy<0) ipy=0;
if (ipy>=nSubPixels) ipy=nSubPixels-1;
for (int ibx = 0; ibx < nbeta; ibx++) {
for (int iby = 0; iby < nbeta; iby++) {
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;
// }
ipy=hhx[ibx+iby*nbeta]*nSubPixels;
//if (hhx[ibx+iby*nbeta]>=((ipy)*bsize) && hhx[ibx+iby*nbeta]<=((ipy+1)*bsize)) {
if (ipy<0) ipy=0;
if (ipy>=nSubPixels) ipy=nSubPixels-1;
// if ( hhy[ibx+iby*nbeta]>=((ipy)*bsize) &&
// hhy[ibx+iby*nbeta]<=((ipy+1)*bsize)) {
ipy = hhy[ibx + iby * nbeta] * nSubPixels;
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 < 0)
ipy = 0;
if (ipy >= nSubPixels)
ipy = nSubPixels - 1;
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){
// flat=new double[nSubPixels*nSubPixels]; flat_x=new double[nSubPixels]; flat_y=new double[nSubPixels];
// flat=new double[nSubPixels*nSubPixels];
};
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;
// }
ipy = hhx[ibx + iby * nbeta] * nSubPixels;
// if (hhx[ibx+iby*nbeta]>=((ipy)*bsize) &&
// hhx[ibx+iby*nbeta]<=((ipy+1)*bsize)) {
if (ipy < 0)
ipy = 0;
if (ipy >= nSubPixels)
ipy = nSubPixels - 1;
etaInterpolationAdaptiveBins(etaInterpolationAdaptiveBins *orig): etaInterpolationPosXY(orig){hintcorr=new int[nPixelsX*nPixelsY*nSubPixels];};
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;
// }
}
}
// }
}
virtual etaInterpolationAdaptiveBins* Clone()=0;
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){
// flat=new double[nSubPixels*nSubPixels]; flat_x=new
// double[nSubPixels]; flat_y=new double[nSubPixels];
// flat=new double[nSubPixels*nSubPixels];
};
/* return new etaInterpolationAdaptiveBins(this); */
etaInterpolationAdaptiveBins(etaInterpolationAdaptiveBins *orig)
: etaInterpolationPosXY(orig) {
hintcorr = new int[nPixelsX * nPixelsY * nSubPixels];
};
/* }; */
virtual etaInterpolationAdaptiveBins *Clone() = 0;
/* return new etaInterpolationAdaptiveBins(this); */
/* }; */
virtual void prepareInterpolation(int &ok) {
prepareInterpolation(ok, 1000);
}
virtual void prepareInterpolation(int &ok) {
prepareInterpolation(ok, 1000);
}
virtual void prepareInterpolation(int &ok, int nint)
{
ok=1;
virtual void prepareInterpolation(int &ok, int nint) {
ok = 1;
///*Eta Distribution Rebinning*///
double bsize=1./nSubPixels; //precision
// cout<<"nPixelsX = "<<nPixelsX<<" nPixelsY = "<<nPixelsY<<" nSubPixels = "<<nSubPixels<<endl;
double tot_eta=0;
double tot_eta_x=0;
double tot_eta_y=0;
for (int ip=0; ip<nbeta*nbeta; ip++)
tot_eta+=heta[ip];
if (tot_eta<=0) {ok=0; return;};
///*Eta Distribution Rebinning*///
double bsize = 1. / nSubPixels; // precision
// cout<<"nPixelsX = "<<nPixelsX<<" nPixelsY = "<<nPixelsY<<" nSubPixels
// = "<<nSubPixels<<endl;
double tot_eta = 0;
double tot_eta_x = 0;
double tot_eta_y = 0;
for (int ip = 0; ip < nbeta * nbeta; ip++)
tot_eta += heta[ip];
if (tot_eta <= 0) {
ok = 0;
return;
};
double hx[nbeta]; // profile x
double hy[nbeta]; // profile y
double hix[nbeta]; // integral of projection x
double hiy[nbeta]; // integral of projection y
int ii = 0;
double hx[nbeta]; //profile x
double hy[nbeta]; //profile y
double hix[nbeta]; //integral of projection x
double hiy[nbeta]; //integral of projection y
int ii=0;
/** 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)ib)/((float)nbeta);
// }
// }
etaInterpolationPosXY::prepareInterpolation(ok);
/** 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)ib)/((float)nbeta);
// }
// }
etaInterpolationPosXY::prepareInterpolation(ok);
double thr = 1. / ((double)nSubPixels);
double avg = tot_eta / ((double)(nSubPixels * nSubPixels));
double rms = sqrt(tot_eta);
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;
// cout << " chi2= " << old_diff << " (rms= " << sqrt(tot_eta) << ")" <<
// endl;
cout << endl;
cout << endl;
debugSaveAll(0);
int iint = 0;
float *newhhx = new float[nbeta * nbeta]; // profile x
float *newhhy = new float[nbeta * nbeta]; // profile y
float *besthhx = hhx; // profile x
float *besthhy = hhy; // profile y
cout << "Iteration " << iint << " Chi2: " << old_diff
<< endl; //" Best: "<< best_diff << " RMS: "<< rms<< endl;
while (iint < nint && best_diff > rms) {
double thr=1./((double)nSubPixels);
double avg=tot_eta/((double)(nSubPixels*nSubPixels));
double rms=sqrt(tot_eta);
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;
// cout << " chi2= " << old_diff << " (rms= " << sqrt(tot_eta) << ")" << endl;
cout << endl;
cout << endl;
debugSaveAll(0);
int iint=0;
float *newhhx=new float[nbeta*nbeta]; //profile x
float *newhhy=new float[nbeta*nbeta]; //profile y
float *besthhx=hhx; //profile x
float *besthhy=hhy; //profile y
/* #ifdef SAVE_ALL */
/* if (iint%10==0) */
/* debugSaveAll(iint); */
/* #endif */
// cout << "Iteration " << iint << endl;
iterate(newhhx, newhhy);
new_diff = calcDiff(avg, newhhx, newhhy);
// cout << " chi2= " << new_diff << " (rms= " << sqrt(tot_eta) <<
// ")"<<endl;
cout << "Iteration "<< iint << " Chi2: " << old_diff << endl; //" Best: "<< best_diff << " RMS: "<< rms<< endl;
while (iint<nint && best_diff > rms) {
if (new_diff < best_diff) {
best_diff = new_diff;
besthhx = newhhx;
besthhy = newhhy;
}
/* #ifdef SAVE_ALL */
/* if (iint%10==0) */
/* debugSaveAll(iint); */
/* #endif */
// cout << "Iteration " << iint << endl;
iterate(newhhx,newhhy);
new_diff=calcDiff(avg, newhhx, newhhy);
// cout << " chi2= " << new_diff << " (rms= " << sqrt(tot_eta) << ")"<<endl;
if (hhx != besthhx)
delete[] hhx;
if (hhy != besthhy)
delete[] hhy;
if (new_diff<best_diff) {
best_diff=new_diff;
besthhx=newhhx;
besthhy=newhhy;
}
if (hhx!=besthhx)
delete [] hhx;
if (hhy!=besthhy)
delete [] hhy;
hhx=newhhx;
hhy=newhhy;
hhx = newhhx;
hhy = newhhy;
#ifdef SAVE_ALL
if (new_diff<=best_diff) {
debugSaveAll(iint);
}
if (new_diff <= best_diff) {
debugSaveAll(iint);
}
#endif
newhhx = new float[nbeta * nbeta]; // profile x
newhhy = new float[nbeta * nbeta]; // profile y
newhhx=new float[nbeta*nbeta]; //profile x
newhhy=new float[nbeta*nbeta]; //profile y
/* 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;
// }
/* 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;
old_diff = new_diff;
iint++;
cout << "Iteration "<< iint << " Chi2: " << new_diff << endl; //" Best: "<< best_diff << " RMS: "<< rms<< endl;
}
delete [] newhhx;
delete [] newhhy;
if (hhx!=besthhx)
delete [] hhx;
if (hhy!=besthhy)
delete [] hhy;
hhx=besthhx;
hhy=besthhy;
iint++;
cout << "Iteration " << iint << " Chi2: " << new_diff
<< endl; //" Best: "<< best_diff << " RMS: "<< rms<< endl;
}
delete[] newhhx;
delete[] newhhy;
cout << "Iteration "<< iint << " Chi2: " << best_diff << endl; //" Best: "<< best_diff << " RMS: "<< rms<< endl;
if (hhx != besthhx)
delete[] hhx;
if (hhy != besthhy)
delete[] hhy;
hhx = besthhx;
hhy = besthhy;
cout << "Iteration " << iint << " Chi2: " << best_diff
<< endl; //" Best: "<< best_diff << " RMS: "<< rms<< endl;
#ifdef SAVE_ALL
debugSaveAll(iint);
debugSaveAll(iint);
#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 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){};
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);};
virtual eta3InterpolationAdaptiveBins *Clone() {
return new eta3InterpolationAdaptiveBins(this);
};
};
#endif