Interpolation changed for rectangular pixels

This commit is contained in:
2019-12-09 13:53:11 +01:00
parent a2c26664ee
commit 9584ea626d
7 changed files with 331 additions and 291 deletions

View File

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