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

merging refactor (replacing)

Browse Source
This commit is contained in:
maliakal_d
2019-04-12 10:53:09 +02:00
parent 0bb800cc8a
commit 89a06f099c
1176 changed files with 82698 additions and 159058 deletions
Download Patch File Download Diff File Expand all files Collapse all files

BIN
slsDetectorCalibration/interpolations/.__afs3A48
View File

Binary file not shown.

403
slsDetectorCalibration/interpolations/eta2InterpolationBase.h
View File

@ -1,403 +0,0 @@
#ifndef ETA2_INTERPOLATION_BASE_H
#define ETA2_INTERPOLATION_BASE_H
#ifdef MYROOT1
#include <TObject.h>
#include <TTree.h>
#include <TH2D.h>
#include <TH2F.h>
#endif
#include "etaInterpolationBase.h"
class eta2InterpolationBase : public virtual etaInterpolationBase {
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) {
/* if (etamin>=etamax) { */
/* etamin=-1; */
/* etamax=2; */
/* // cout << ":" <<endl; */
/* } */
/* etastep=(etamax-etamin)/nbeta; */
};
eta2InterpolationBase(eta2InterpolationBase *orig): etaInterpolationBase(orig){ };
//////////////////////////////////////////////////////////////////////////////
//////////// /*It return position hit for the event in input */ //////////////
virtual void getInterpolatedPosition(int x, int y, int *data, double &int_x, double &int_y)
{
double sDum[2][2];
double tot, totquad;
double etax,etay;
int corner;
corner=calcQuad(data, tot, totquad, sDum);
if (nSubPixels>2)
calcEta(totquad, sDum, etax, etay);
getInterpolatedPosition(x,y,etax,etay,corner,int_x,int_y);
return;
};
virtual void getInterpolatedPosition(int x, int y, double *data, double &int_x, double &int_y)
{
double sDum[2][2];
double tot, totquad;
double etax,etay;
int corner;
corner=calcQuad(data, tot, totquad, sDum);
if (nSubPixels>2)
calcEta(totquad, sDum, etax, etay);
getInterpolatedPosition(x,y,etax,etay,corner,int_x,int_y);
return;
};
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cl,double &int_x, double &int_y) {
double cc[2][2];
int xoff, yoff;
switch (quad) {
case BOTTOM_LEFT:
xoff=0;
yoff=0;
break;
case BOTTOM_RIGHT:
xoff=1;
yoff=0;
break;
case TOP_LEFT:
xoff=0;
yoff=1;
break;
case TOP_RIGHT:
xoff=1;
yoff=1;
break;
default:
;
}
double etax, etay;
if (nSubPixels>2) {
cc[0][0]=cl[xoff+3*yoff];
cc[1][0]=cl[xoff+3*(yoff+1)];
cc[0][1]=cl[xoff+1+3*yoff];
cc[1][1]=cl[xoff+1+3*(yoff+1)];
calcEta(totquad,cc,etax,etay);
}
return getInterpolatedPosition(x,y,etax, etay,quad,int_x,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];
int xoff, yoff;
switch (quad) {
case BOTTOM_LEFT:
xoff=0;
yoff=0;
break;
case BOTTOM_RIGHT:
xoff=1;
yoff=0;
break;
case TOP_LEFT:
xoff=0;
yoff=1;
break;
case TOP_RIGHT:
xoff=1;
yoff=1;
break;
default:
;
}
double etax, etay;
if (nSubPixels>2) {
cc[0][0]=cl[xoff+3*yoff];
cc[1][0]=cl[xoff+3*(yoff+1)];
cc[0][1]=cl[xoff+1+3*yoff];
cc[1][1]=cl[xoff+1+3*(xoff+1)];
calcEta(totquad,cc,etax,etay);
}
return getInterpolatedPosition(x,y,etax, etay,quad,int_x,int_y);
}
virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int corner, double &int_x, double &int_y)
{
double xpos_eta=0,ypos_eta=0;
double dX,dY;
int ex,ey;
switch (corner)
{
case TOP_LEFT:
dX=-1.;
dY=0;
break;
case TOP_RIGHT:
;
dX=0;
dY=0;
break;
case BOTTOM_LEFT:
dX=-1.;
dY=-1.;
break;
case BOTTOM_RIGHT:
dX=0;
dY=-1.;
break;
default:
cout << "bad quadrant" << endl;
dX=0.;
dY=0.;
}
if (nSubPixels>2) {
ex=(etax-etamin)/etastep;
ey=(etay-etamin)/etastep;
if (ex<0) {
cout << "x*"<< ex << endl;
ex=0;
}
if (ex>=nbeta) {
cout << "x?"<< ex << endl;
ex=nbeta-1;
}
if (ey<0) {
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);
ypos_eta=(((double)hhy[(ey*nbeta+ex)]))+dY ;///((double)nSubPixels);
} else {
xpos_eta=0.5*dX+0.25;
ypos_eta=0.5*dY+0.25;
}
int_x=((double)x) + xpos_eta+0.5;
int_y=((double)y) + ypos_eta+0.5;
}
virtual int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay) {
double cc[2][2];
int xoff, yoff;
switch (quad) {
case BOTTOM_LEFT:
xoff=0;
yoff=0;
break;
case BOTTOM_RIGHT:
xoff=1;
yoff=0;
break;
case TOP_LEFT:
xoff=0;
yoff=1;
break;
case TOP_RIGHT:
xoff=1;
yoff=1;
break;
default:
;
}
cc[0][0]=cl[xoff+3*yoff];
cc[1][0]=cl[xoff+3*(yoff+1)];
cc[0][1]=cl[xoff+1+3*yoff];
cc[1][1]=cl[xoff+1+3*(yoff+1)];
//calcMyEta(totquad,quad,cl,etax, etay);
calcEta(totquad, cc,etax, etay);
// cout <<"******"<< etax << " " << etay << endl;
return addToFlatField(etax,etay);
}
virtual int addToFlatField(double totquad,int quad,double *cl,double &etax, double &etay) {
double cc[2][2];
int xoff, yoff;
switch (quad) {
case BOTTOM_LEFT:
xoff=0;
yoff=0;
break;
case BOTTOM_RIGHT:
xoff=1;
yoff=0;
break;
case TOP_LEFT:
xoff=0;
yoff=1;
break;
case TOP_RIGHT:
xoff=1;
yoff=1;
break;
default:
;
}
cc[0][0]=cl[xoff+3*yoff];
cc[1][0]=cl[(yoff+1)*3+xoff];
cc[0][1]=cl[yoff*3+xoff+1];
cc[1][1]=cl[(yoff+1)*3+xoff+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);
calcEta(totquad, cc,etax, etay);
// cout <<"******"<< etax << " " << etay << endl;
return addToFlatField(etax,etay);
}
//////////////////////////////////////////////////////////////////////////////////////
virtual int addToFlatField(double *cluster, double &etax, double &etay){
double sDum[2][2];
double tot, totquad;
int corner;
corner=calcQuad(cluster, tot, totquad, sDum);
double xpos_eta,ypos_eta;
double dX,dY;
calcEta(totquad, sDum, etax, etay);
return addToFlatField(etax,etay);
};
virtual int addToFlatField(int *cluster, double &etax, double &etay){
double sDum[2][2];
double tot, totquad;
int corner;
corner=calcQuad(cluster, tot, totquad, sDum);
double xpos_eta,ypos_eta;
double dX,dY;
calcEta(totquad, sDum, etax, etay);
return addToFlatField(etax,etay);
};
virtual int addToFlatField(double etax, double etay){
#ifdef MYROOT1
heta->Fill(etax,etay);
#endif
#ifndef MYROOT1
int ex,ey;
ex=(etax-etamin)/etastep;
ey=(etay-etamin)/etastep;
if (ey<nbeta && ex<nbeta && ex>=0 && ey>=0)
heta[ey*nbeta+ex]++;
#endif
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: */
/* #ifdef MYROOT1 */
/* TH2D *heta; */
/* TH2D *hhx; */
/* TH2D *hhy; */
/* #endif */
/* #ifndef MYROOT1 */
/* int *heta; */
/* float *hhx; */
/* float *hhy; */
/* #endif */
/* int nbeta; */
/* double etamin, etamax, etastep; */
};
#endif

294
slsDetectorCalibration/interpolations/eta3InterpolationBase.h
View File

@ -1,294 +0,0 @@
#ifndef ETA3_INTERPOLATION_BASE_H
#define ETA3_INTERPOLATION_BASE_H
#ifdef MYROOT1
#include <TObject.h>
#include <TTree.h>
#include <TH2D.h>
#include <TH2F.h>
#endif
#include "etaInterpolationBase.h"
class eta3InterpolationBase : public virtual etaInterpolationBase {
public:
eta3InterpolationBase(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 << "e3ib " << nb << " " << emin << " " << emax << endl;
/* if (nbeta<=0) { */
/* nbeta=nSubPixels*10; */
/* } */
if (etamin>=etamax) {
etamin=-1;
etamax=1;
}
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;
};
eta3InterpolationBase(eta3InterpolationBase *orig): etaInterpolationBase(orig){ };
/* virtual eta3InterpolationBase* Clone()=0; */
// virtual void prepareInterpolation(int &ok){};
//////////////////////////////////////////////////////////////////////////////
//////////// /*It return position hit for the event in input */ //////////////
virtual void getInterpolatedPosition(int x, int y, int *data, double &int_x, double &int_y)
{
double tot, totquad;
double etax,etay;
int corner=calcEta3(data,etax,etay, totquad);
getInterpolatedPosition(x,y,etax,etay,corner,int_x,int_y);
return;
};
virtual void getInterpolatedPosition(int x, int y, double *data, double &int_x, double &int_y)
{
double sDum[2][2];
double tot, totquad;
double etax,etay;
int corner=calcEta3(data,etax,etay, totquad);
getInterpolatedPosition(x,y,etax,etay,corner,int_x,int_y);
return;
};
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cl,double &int_x, double &int_y) {
double etax, etay;
if (nSubPixels>2) {
calcEta3(cl,etax,etay, totquad);
}
return getInterpolatedPosition(x,y,etax, etay,quad,int_x,int_y);
}
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,int *cl,double &int_x, double &int_y) {
double etax, etay;
if (nSubPixels>2) {
calcEta3(cl,etax,etay, totquad);
}
return getInterpolatedPosition(x,y,etax, etay,quad,int_x,int_y);
}
virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int corner, double &int_x, double &int_y)
{
double xpos_eta=0,ypos_eta=0;
int ex,ey;
if (nSubPixels>2) {
#ifdef MYROOT1
xpos_eta=(hhx->GetBinContent(hhx->GetXaxis()->FindBin(etax),hhy->GetYaxis()->FindBin(etay)))/((double)nSubPixels);
ypos_eta=(hhy->GetBinContent(hhx->GetXaxis()->FindBin(etax),hhy->GetYaxis()->FindBin(etay)))/((double)nSubPixels);
#endif
#ifndef MYROOT1
ex=(etax-etamin)/etastep;
ey=(etay-etamin)/etastep;
if (ex<0) {
/* cout << etax << " " << etamin << " "; */
/* cout << "3x*"<< ex << endl; */
ex=0;
}
if (ex>=nbeta) {
/* cout << etax << " " << etamin << " "; */
/* cout << "3x?"<< ex << endl; */
ex=nbeta-1;
}
if (ey<0) {
/* cout << etay << " " << etamin << " "; */
/* cout << "3y*"<< ey << endl; */
ey=0;
}
if (ey>=nbeta) {
/* cout << etay << " " << etamin << " "; */
/* cout << "3y?"<< ey << endl; */
ey=nbeta-1;
}
xpos_eta=(((double)hhx[(ey*nbeta+ex)]));///((double)nSubPixels);
ypos_eta=(((double)hhy[(ey*nbeta+ex)]));///((double)nSubPixels);
#endif
} else {
switch (corner) {
case BOTTOM_LEFT:
xpos_eta=-0.25;
ypos_eta=-0.25;
break;
case BOTTOM_RIGHT:
xpos_eta=0.25;
ypos_eta=-0.25;
break;
case TOP_LEFT:
xpos_eta=-0.25;
ypos_eta=0.25;
break;
case TOP_RIGHT:
xpos_eta=0.25;
ypos_eta=0.25;
break;
default:
xpos_eta=0;
ypos_eta=0;
}
}
int_x=((double)x) + xpos_eta;
int_y=((double)y) + ypos_eta;
// int_x=5. + xpos_eta;
// int_y=5. + ypos_eta;
}
/* ///////////////////////////////////////////////////////////////////////////////////////////////// */
/* virtual void getPositionETA3(int x, int y, double *data, double &int_x, double &int_y) */
/* { */
/* double sDum[2][2]; */
/* double tot, totquad; */
/* double eta3x,eta3y; */
/* double ex,ey; */
/* calcQuad(data, tot, totquad, sDum); */
/* calcEta3(data,eta3x, eta3y,tot); */
/* double xpos_eta,ypos_eta; */
/* #ifdef MYROOT1 */
/* xpos_eta=((hhx->GetBinContent(hhx->GetXaxis()->FindBin(eta3x),hhy->GetYaxis()->FindBin(eta3y))))/((double)nSubPixels); */
/* ypos_eta=((hhy->GetBinContent(hhx->GetXaxis()->FindBin(eta3x),hhy->GetYaxis()->FindBin(eta3y))))/((double)nSubPixels); */
/* #endif */
/* #ifndef MYROOT1 */
/* ex=(eta3x-etamin)/etastep; */
/* ey=(eta3y-etamin)/etastep; */
/* if (ex<0) ex=0; */
/* if (ex>=nbeta) ex=nbeta-1; */
/* if (ey<0) ey=0; */
/* if (ey>=nbeta) ey=nbeta-1; */
/* xpos_eta=(((double)hhx[(int)(ey*nbeta+ex)]))/((double)nSubPixels); */
/* ypos_eta=(((double)hhy[(int)(ey*nbeta+ex)]))/((double)nSubPixels); */
/* #endif */
/* int_x=((double)x) + xpos_eta; */
/* int_y=((double)y) + ypos_eta; */
/* return; */
/* }; */
virtual int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay) {
calcEta3(cl, etax, etay, totquad);
return addToFlatField(etax,etay);
}
virtual int addToFlatField(double totquad,int quad,double *cl,double &etax, double &etay) {
calcEta3(cl, etax, etay, totquad);
return addToFlatField(etax,etay);
}
//////////////////////////////////////////////////////////////////////////////////////
virtual int addToFlatField(double *cluster, double &etax, double &etay){
double totquad;
calcEta3(cluster, etax, etay, totquad);
return addToFlatField(etax,etay);
};
virtual int addToFlatField(int *cluster, double &etax, double &etay){
double totquad;
calcEta3(cluster, etax, etay, totquad);
return addToFlatField(etax,etay);
};
virtual int addToFlatField(double etax, double etay){
#ifdef MYROOT1
heta->Fill(etax,etay);
#endif
#ifndef MYROOT1
int ex,ey;
ex=(etax-etamin)/etastep;
ey=(etay-etamin)/etastep;
if (ey<nbeta && ex<nbeta && ex>=0 && ey>=0)
heta[ey*nbeta+ex]++;
#endif
return 0;
};
/* protected: */
/* #ifdef MYROOT1 */
/* TH2D *heta; */
/* TH2D *hhx; */
/* TH2D *hhy; */
/* #endif */
/* #ifndef MYROOT1 */
/* int *heta; */
/* float *hhx; */
/* float *hhy; */
/* #endif */
/* int nbeta; */
/* double etamin, etamax, etastep; */
};
#endif

285
slsDetectorCalibration/interpolations/etaInterpolationAdaptiveBins.h
View File

@ -1,285 +0,0 @@
#ifndef ETA_INTERPOLATION_ADAPTIVEBINS_H
#define ETA_INTERPOLATION_ADAPTIVEBINS_H
#include <cmath>
#include "tiffIO.h"
//#include "etaInterpolationBase.h"
#include "etaInterpolationPosXY.h"
class etaInterpolationAdaptiveBins : public etaInterpolationPosXY {
// protected:
private:
virtual void iterate(float *newhhx, float *newhhy) {
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;
}
}
// 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];
}
}
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;
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;
// }
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 (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;
// }
}
}
// }
}
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];
};
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, 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;};
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);
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) {
/* #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 (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;
#ifdef SAVE_ALL
if (new_diff<=best_diff) {
debugSaveAll(iint);
}
#endif
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;
// }
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;
cout << "Iteration "<< iint << " Chi2: " << best_diff << endl; //" Best: "<< best_diff << " RMS: "<< rms<< endl;
#ifdef SAVE_ALL
debugSaveAll(iint);
#endif
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

369
slsDetectorCalibration/interpolations/etaInterpolationBase.h
View File

@ -1,369 +0,0 @@
#ifndef ETA_INTERPOLATION_BASE_H
#define ETA_INTERPOLATION_BASE_H
#ifdef MYROOT1
#include <TObject.h>
#include <TTree.h>
#include <TH2D.h>
#include <TH2F.h>
#endif
#include "slsInterpolation.h"
#include "tiffIO.h"
class etaInterpolationBase : public slsInterpolation {
public:
etaInterpolationBase(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : slsInterpolation(nx,ny,ns), hhx(NULL), hhy(NULL), heta(NULL), nbeta(nb), etamin(emin), etamax(emax) {
// cout << "eb " << nb << " " << emin << " " << emax << endl;
// cout << nb << " " << etamin << " " << etamax << endl;
if (nbeta<=0) {
//cout << "aaa:" <<endl;
nbeta=nSubPixels*10;
}
if (etamin>=etamax) {
etamin=-1;
etamax=2;
}
etastep=(etamax-etamin)/nbeta;
heta=new int[nbeta*nbeta];
hhx=new float[nbeta*nbeta];
hhy=new float[nbeta*nbeta];
rangeMin=etamin;
rangeMax=etamax;
flat= new double[nSubPixels*nSubPixels];
hintcorr=new int [nSubPixels*nSubPixels*nPixelsX*nPixelsY];
};
etaInterpolationBase(etaInterpolationBase *orig): slsInterpolation(orig){
nbeta=orig->nbeta;
etamin=orig->etamin;
etamax=orig->etamax;
rangeMin=orig->rangeMin;
rangeMax=orig->rangeMax;
etastep=(etamax-etamin)/nbeta;
heta=new int[nbeta*nbeta];
memcpy(heta,orig->heta,nbeta*nbeta*sizeof(int));
hhx=new float[nbeta*nbeta];
memcpy(hhx,orig->hhx,nbeta*nbeta*sizeof(float));
hhy=new float[nbeta*nbeta];
memcpy(hhy,orig->hhy,nbeta*nbeta*sizeof(float));
hintcorr=new int [nSubPixels*nSubPixels*nPixelsX*nPixelsY];
};
virtual void resetFlatField() {
for (int ibx=0; ibx<nbeta*nbeta; ibx++) {
heta[ibx]=0;
hhx[ibx]=0;
hhy[ibx]=0;
}
};
int *setEta(int *h, int nb=-1, double emin=1, double emax=0)
{
if (h) {
if (heta) delete [] heta;
heta=h;
nbeta=nb;
if (nb<=0) nbeta=nSubPixels*10;
etamin=emin;
etamax=emax;
if (etamin>=etamax) {
etamin=-1;
etamax=2;
}
rangeMin=etamin;
rangeMax=etamax;
etastep=(etamax-etamin)/nbeta;
}
return heta;
};
int *setFlatField(int *h, int nb=-1, double emin=1, double emax=0)
{
return setEta(h, nb, emin, emax);
};
int *getFlatField(){return setEta(NULL);};
int *getFlatField(int &nb, double &emin, double &emax){
nb=nbeta;
emin=etamin;
emax=etamax;
return getFlatField();
};
void *writeFlatField(const char * imgname) {
float *gm=NULL;
gm=new float[nbeta*nbeta];
for (int ix=0; ix<nbeta; ix++) {
for (int iy=0; iy<nbeta; iy++) {
gm[iy*nbeta+ix]=heta[iy*nbeta+ix];
}
}
WriteToTiff(gm, imgname, nbeta, nbeta);
delete [] gm;
return NULL;
};
int readFlatField(const char * imgname, double emin=1, double emax=0) {
if (emax>=1) etamax=emax;
if (emin<=0) etamin=emin;
if (etamin>=etamax) {
etamin=-1;
etamax=2;
}
etastep=(etamax-etamin)/nbeta;
uint32 nnx;
uint32 nny;
float *gm=ReadFromTiff(imgname, nnx, nny);
if (nnx!=nny) {
cout << "different number of bins in x " << nnx << " and y " << nny<< " !"<< endl;
cout << "Aborting read"<< endl;
return 0;
}
nbeta=nnx;
if (gm) {
if (heta) {
delete [] heta;
delete [] hhx;
delete [] hhy;
}
heta=new int[nbeta*nbeta];
hhx=new float[nbeta*nbeta];
hhy=new float[nbeta*nbeta];
for (int ix=0; ix<nbeta; ix++) {
for (int iy=0; iy<nbeta; iy++) {
heta[iy*nbeta+ix]=gm[iy*nbeta+ix];
}
}
delete [] gm;
return 1;
}
return 0;
};
float *gethhx()
{
// hhx->Scale((double)nSubPixels);
return hhx;
};
float *gethhy()
{
// hhy->Scale((double)nSubPixels);
return hhy;
};
virtual int addToFlatField(double etax, double etay){
int ex,ey;
ex=(etax-etamin)/etastep;
ey=(etay-etamin)/etastep;
if (ey<nbeta && ex<nbeta && ex>=0 && ey>=0)
heta[ey*nbeta+ex]++;
return 0;
};
// virtual void prepareInterpolation(int &ok)=0;
void debugSaveAll(int ind=0) {
int 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:
double calcDiff(double avg, float *hx, float *hy) {
//double p_tot=0;
double diff=0, d;
double bsize=1./nSubPixels;
int nbad=0;
double p_tot_x[nSubPixels], p_tot_y[nSubPixels], p_tot[nSubPixels*nSubPixels];
double maxdiff=0, mindiff=avg*nSubPixels*nSubPixels;
int ipx, ipy;
for (ipy=0; ipy<nSubPixels; ipy++) {
for (ipx=0; ipx<nSubPixels; ipx++) {
p_tot[ipx+ipy*nSubPixels]=0;
}
p_tot_y[ipy]=0;
p_tot_x[ipy]=0;
}
for (int ibx=0; ibx<nbeta; ibx++) {
for (int iby=0; iby<nbeta; iby++) {
ipx=hx[ibx+iby*nbeta]*nSubPixels;
if (ipx<0) ipx=0;
if (ipx>=nSubPixels) ipx=nSubPixels-1;
ipy=hy[ibx+iby*nbeta]*nSubPixels;
if (ipy<0) ipy=0;
if (ipy>=nSubPixels) ipy=nSubPixels-1;
p_tot[ipx+ipy*nSubPixels]+=heta[ibx+iby*nbeta];
p_tot_y[ipy]+=heta[ibx+iby*nbeta];
p_tot_x[ipx]+=heta[ibx+iby*nbeta];
}
}
// cout << endl << endl;
for (ipy=0; ipy<nSubPixels; ipy++) {
cout.width(5);
//flat_y[ipy]=p_tot_y[ipy];//avg/nSubPixels;
for (ipx=0; ipx<nSubPixels; ipx++) {
// flat_x[ipx]=p_tot_x[ipx];///avg/nSubPixels;
flat[ipx+nSubPixels*ipy]=p_tot[ipx+nSubPixels*ipy];///avg;
d=p_tot[ipx+nSubPixels*ipy]-avg;
if (d<0) d*=-1.;
if (d>5*sqrt(avg) )
nbad++;
diff+=d*d;
if (d<mindiff) mindiff=d;
if (d>maxdiff) maxdiff=d;
// cout << setprecision(4) << p_tot[ipx+nSubPixels*ipy] << " ";
}
/* cout << "** " << setprecision(4) << flat_y[ipy]; */
//cout << "\n";
}
/* cout << "**" << endl; cout.width(5); */
/* for (ipx=0; ipx<nSubPixels; ipx++) { */
/* cout << setprecision(4) << flat_x[ipx] << " "; */
/* } */
//cout << "**" << endl; cout.width(5);
//cout << "Min diff: " << mindiff/sqrt(avg) << " Max diff: " << maxdiff/sqrt(avg) << " Nbad: " << nbad << endl;
// cout << "Bad pixels: " << 100.*(float)nbad/((float)(nSubPixels*nSubPixels)) << " %" << endl;
return sqrt(diff);
}
int *heta;
float *hhx;
float *hhy;
int nbeta;
double etamin, etamax, etastep;
double rangeMin, rangeMax;
double *flat;
int *hintcorr;
};
#endif

263
slsDetectorCalibration/interpolations/etaInterpolationCleverAdaptiveBins.h
View File

@ -1,263 +0,0 @@
#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

85
slsDetectorCalibration/interpolations/etaInterpolationGlobal.h
View File

@ -1,85 +0,0 @@
#ifndef ETA_INTERPOLATION_GLOBAL_H
#define ETA_INTERPOLATION_GLOBAL_H
#include "etaInterpolationBase.h"
class etaInterpolationGlobal : public etaInterpolationBase{
public:
globalEtaInterpolation(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){};
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*///
double bsize=1./nSubPixels; //precision
// cout<<"nPixelsX = "<<nPixelsX<<" nPixelsY = "<<nPixelsY<<" nSubPixels = "<<nSubPixels<<endl;
double tot_eta=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;};
double hx[nbeta]; //projection x
double hy[nbeta]; //projection y
for (int ibx=0; ibx<nbeta; ibx++) {
for (int iby=0; iby<nbeta; iby++) {
hx[ibx]=hx[ibx]+heta[ibx+iby*nbeta];
hy[iby]=hx[iby]+heta[ibx+iby*nbeta];
}
}
double hix[nbeta]; //integral of projection x
double hiy[nbeta]; //integral of projection y
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]+hx[ib];
}
int ib=0;
for (int ibx=0; ibx<nbeta; ibx++) {
if (hix[ibx]>(ib+1)*tot_eta*bsize) ib++;
for (int iby=0; iby<nbeta; iby++) {
#ifdef MYROOT1
hhx->SetBinContent(ibx+1,iby+1,ib);
#endif
#ifndef MYROOT1
hhx[ibx+iby*nbeta]=ib;
#endif
}
}
ib=0;
for (int iby=0; iby<nbeta; iby++) {
if (hiy[iby]>(ib+1)*tot_eta*bsize) ib++;
for (int ibx=0; ibx<nbeta; ibx++) {
#ifdef MYROOT1
hhy->SetBinContent(ibx+1,iby+1,ib);
#endif
#ifndef MYROOT1
hhy[ibx+iby*nbeta]=ib;
#endif
}
}
return ;
};
};
#endif

184
slsDetectorCalibration/interpolations/etaInterpolationPosXY.h
View File

@ -1,184 +0,0 @@
#ifndef ETA_INTERPOLATION_POSXY_H
#define ETA_INTERPOLATION_POSXY_H
//#include "tiffIO.h"
#include "etaInterpolationBase.h"
#include "eta2InterpolationBase.h"
#include "eta3InterpolationBase.h"
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){
// cout << "epxy " << nb << " " << emin << " " << emax << endl; cout << nbeta << " " << etamin << " " << etamax << endl;
};
etaInterpolationPosXY(etaInterpolationPosXY *orig): etaInterpolationBase(orig) {};
virtual etaInterpolationPosXY* Clone()=0; /* { */
/* return new etaInterpolationPosXY(this); */
/* }; */
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*///
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];
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
int ii=0;
double etax, etay;
for (int ib=0; ib<nbeta; ib++) {
tot_eta_x=0;
tot_eta_y=0;
for (int iby=0; iby<nbeta; iby++) {
etax=etamin+iby*etastep;
//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];
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];
hiy[iby]=hiy[iby-1]+hy[iby];
}
ii=0;
tot_eta_x=hix[nbeta-1]+1;
tot_eta_y=hiy[nbeta-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++) {
if (tot_eta_y<=0) {
hhy[ib+ibx*nbeta]=-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;
}
}
}
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
debugSaveAll();
#endif
return ;
}
};
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){
// cout << "e2pxy " << nb << " " << emin << " " << emax << endl;
};
eta2InterpolationPosXY(eta2InterpolationPosXY *orig): etaInterpolationBase(orig), etaInterpolationPosXY(orig) {};
virtual eta2InterpolationPosXY* Clone() { return new eta2InterpolationPosXY(this);};
};
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(eta3InterpolationPosXY *orig): etaInterpolationBase(orig), etaInterpolationPosXY(orig) {};
virtual eta3InterpolationPosXY* Clone() { return new eta3InterpolationPosXY(this);};
};
#endif

417
slsDetectorCalibration/interpolations/etaInterpolationRandomBins.h
View File

@ -1,417 +0,0 @@
#ifndef ETA_INTERPOLATION_RANDOMBINS_H
#define ETA_INTERPOLATION_RANDOMBINS_H
#include "tiffIO.h"
//#include "etaInterpolationBase.h"
#include "etaInterpolationPosXY.h"
#include <cstdlib>
#include <algorithm>
//#include <math>
#include <cmath> // std::abs
#define PI 3.14159265
#define TWOPI 2.*PI
using namespace std;
class etaInterpolationRandomBins : public etaInterpolationPosXY {
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++) {
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;
}
double iterate(float *newhhx, float *newhhy, double avg) {
double bsize=1./nSubPixels;
double hy[nbeta]; //profile y
double hx[nbeta]; //profile x
double hix[nbeta]; //integral of projection x
double hiy[nbeta]; //integral of projection y
double tot_eta_x=0;
double tot_eta_y=0;
int p0;
int vx[(nSubPixels+1)*(nSubPixels+1)], vy[(nSubPixels+1)*(nSubPixels+1)];
int arrx[nSubPixels+1], arry[nSubPixels+1];
int bad=1;
int isby, isbx;
int ii=0;
// using default comparison (operator <):
// std::sort (myvector.begin(), myvector.begin()+4); //(12 32 45 71)26 80 53 33
for (isby=0; isby<(nSubPixels+1)/2+1; isby++) {
for (isbx=0; isbx<(nSubPixels+1)/2+1; isbx++) {
p0=isby*(nSubPixels+1)+isbx;
// for (int iv=0; iv<(nSubPixels+1)*(nSubPixels+1); iv++) {
if (isbx==0) {
vy[p0]=isby*nbeta/nSubPixels;
vx[p0]=0;
} else if ( isby==0 ) {
vy[p0]=0;
vx[p0]=isbx*nbeta/nSubPixels;
}
else {
vy[p0]=rand()%(nbeta/2);
vx[p0]=rand()%(nbeta/2);
if (nSubPixels%2==0 && isbx==nSubPixels/2)
vx[p0]=nbeta/2;
if (nSubPixels%2==0 && isby==nSubPixels/2 )
vy[p0]=nbeta/2;
}
// cout << "(" << vx[p0] << " , " << vy[p0] << " ) \t" ;
// }
}
//cout << endl;
}
// cout << "rand" << endl;
while (bad) {
for (isby=0; isby<(nSubPixels+1)/2+1; isby++) {
for (isbx=0; isbx<(nSubPixels+1)/2+1; isbx++) {
arrx[isbx]=vx[isby*(nSubPixels+1)+isbx];
arry[isbx]=vy[isbx*(nSubPixels+1)+isby];
//cout << isbx << " " << arrx[isbx] << " " << isby << " " << arry[isbx] << endl;
}
sort(arrx,arrx+(nSubPixels+1)/2+1);
sort(arry,arry+(nSubPixels+1)/2+1);
// cout << "*****"<< endl;
// cout << endl;
for (int isbx=0; isbx<(nSubPixels+1)/2+1; isbx++) {
vx[isby*(nSubPixels+1)+isbx]=arrx[isbx];
vy[isbx*(nSubPixels+1)+isby]=arry[isbx];
vx[(nSubPixels-isby)*(nSubPixels+1)+(nSubPixels-isbx)]=nbeta-arrx[isbx];
vy[(nSubPixels-isbx)*(nSubPixels+1)+(nSubPixels-isby)]=nbeta-arry[isbx];
vx[isby*(nSubPixels+1)+(nSubPixels-isbx)]=nbeta-arrx[isbx];
vy[isbx*(nSubPixels+1)+(nSubPixels-isby)]=arry[isbx];
vx[(nSubPixels-isby)*(nSubPixels+1)+(isbx)]=arrx[isbx];
vy[(nSubPixels-isbx)*(nSubPixels+1)+(isby)]=nbeta-arry[isbx];
}
}
/* for (isby=0; isby<nSubPixels+1; isby++) { */
/* for (isbx=0; isbx<nSubPixels+1; isbx++) { */
/* cout << "("<< vx[isby*(nSubPixels+1)+isbx] << " " << vy[isby*(nSubPixels+1)+isbx] << ")\t";//<< endl; */
/* } */
/* cout << endl; */
/* } */
bad=0;
for (isby=1; isby<(nSubPixels+1)/2+1; isby++) {
for (isbx=1; isbx<(nSubPixels+1)/2+1; isbx++) {
if (heta[vx[isby*(nSubPixels+1)+isbx]+vy[isby*(nSubPixels+1)+isbx]*nbeta]<avg*(nSubPixels*nSubPixels)/(nbeta*nbeta)) {
// cout << ii << " " << isbx << " " << isby << " " << vx[isby*(nSubPixels+1)+isbx] << " " << vy[isby*(nSubPixels+1)+isbx] << " " << heta[vx[isby*(nSubPixels+1)+isbx]+vy[isby*(nSubPixels+1)+isbx]*nbeta] << endl;
if (nSubPixels%2==0 && isbx==nSubPixels/2)
;
else
vx[isby*(nSubPixels+1)+isbx]=rand()%(nbeta/2);
if (nSubPixels%2==0 && isbx==nSubPixels/2)
;
else
vy[isby*(nSubPixels+1)+isbx]=rand()%(nbeta/2);
if (bad==0)
ii++;
bad=1;
// break;
}
}
//if (bad) break;
}
// cout << "sort" << endl;
}
cout << ii << " sub iteractions " << avg*(nSubPixels*nSubPixels)/(nbeta*nbeta) << endl;
double m,q;
int in_quad;
int p[4];
int p1x,p2x, p1y, p2y;
// cout << nbeta << endl;
double angle;
double dtheta,theta1,theta2;
for (int ibx=0; ibx<nbeta; ibx++) {
for (int iby=0; iby<nbeta; iby++) {
in_quad=0;
/* if (ibx==0) */
/* isbx=0; */
/* else */
/* isbx= (newhhx[ibx-1+iby*nbeta])/bsize-1; */
/* if (isbx<0) isbx=0; */
/* if (isbx>nSubPixels-1) isbx=nSubPixels-1; */
/* if (iby==0) */
/* isby=0; */
/* else */
/* isby= (newhhx[ibx+(iby-1)*nbeta])/bsize-1; */
/* if (isby<0) isbx=0; */
/* if (isby>nSubPixels-1) isby=nSubPixels-1; */
/* // cout << isbx << " " << isby << endl; */
for (isby=0; isby<nSubPixels; isby++) {
for (isbx=0; isbx<nSubPixels; isbx++) {
// cout << ibx << " " << iby << " " << isbx << " " << isby << endl;
p[0]=isby*(nSubPixels+1)+isbx;
p[1]=isby*(nSubPixels+1)+isbx+1;
p[2]=(isby+1)*(nSubPixels+1)+isbx+1;
p[3]=(isby+1)*(nSubPixels+1)+isbx;
angle=0;
for (int i=0;i<4;i++) {
p1x = vx[p[i]] - ibx;
p1y = vy[p[i]] - iby;
p2x = vx[p[(i+1)%4]] - ibx;
p2y = vy[p[(i+1)%4]] - iby;
theta1 = atan2(p1y,p1x);
theta2 = atan2(p2y,p2x);
dtheta = theta2 - theta1;
while (dtheta > PI)
dtheta -= TWOPI;
while (dtheta < -PI)
dtheta += TWOPI;
angle += dtheta;
}
if (abs((double)angle) < PI)
in_quad=0;
else
in_quad=1;
if (in_quad) {
newhhx[ibx+iby*nbeta]=bsize*((double)isbx);
newhhy[ibx+iby*nbeta]=bsize*((double)isby);
break;
}
}
if (in_quad) break;
}
}
}
// cout << "hist" << endl;
return calcDiff(avg, newhhx, newhhy);
}
public:
etaInterpolationRandomBins(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){};
etaInterpolationRandomBins(etaInterpolationRandomBins *orig): etaInterpolationPosXY(orig){};
virtual etaInterpolationRandomBins* Clone() {
return new etaInterpolationRandomBins(this);
};
virtual void prepareInterpolation(int &ok)
{
ok=1;
cout << "Adaptive bins" << endl;
///*Eta Distribution Rebinning*///
double bsize=1./nSubPixels; //precision
// cout<<"nPixelsX = "<<nPixelsX<<" nPixelsY = "<<nPixelsY<<" nSubPixels = "<<nSubPixels<<endl;
double tot_eta=0;
for (int ip=0; ip<nbeta*nbeta; ip++)
tot_eta+=heta[ip];
if (tot_eta<=0) {ok=0; return;};
int ii=0;
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=-1, new_diff=-1;
// cout << " diff= " << new_diff << endl;
etaInterpolationPosXY::prepareInterpolation(ok);
old_diff=calcDiff(avg, hhx, hhy);
cout << " diff= " << old_diff << endl;
int iint=0;
float *newhhx=new float[nbeta*nbeta]; //profile x
float *newhhy=new float[nbeta*nbeta]; //profile y
int igood=0, ibad=0;
#ifdef SAVE_ALL
int etabins=nbeta;
float *etah=new float[nbeta*nbeta];
char tit[1000];
#endif
while (iint<nint) {
cout << "Iteration " << iint << endl;
new_diff=iterate(newhhx,newhhy, avg);
//new_diff=calcDiff(avg, newhhx, newhhy);
cout << " diff= " << new_diff << " ( " << old_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/randeta_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/randeta_hhy_%d.tiff",iint); */
/* WriteToTiff(etah, tit, etabins, etabins); */
/* #endif */
if (new_diff<old_diff) {
cout << "******************** GOOD! ***********************"<< endl;
delete [] hhx;
delete [] hhy;
igood++;
hhx=newhhx;
hhy=newhhy;
newhhx=new float[nbeta*nbeta]; //profile x */
newhhy=new float[nbeta*nbeta]; //profile y */
old_diff=new_diff;
} else
ibad++;
iint++;
}
delete [] newhhx;
delete [] newhhy;
cout << "performed " << iint << " iterations of which " << igood << " positive " << endl;
/* #ifdef SAVE_ALL */
/* 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 */
return ;
}
};
#endif

678
slsDetectorCalibration/interpolations/etaVEL/EVELAlg.C
View File

@ -1,678 +0,0 @@
#include <TH1D.h>
#include <TH2D.h>
#include <TPad.h>
#include <TDirectory.h>
#include <TEntryList.h>
#include <TFile.h>
#include <TMath.h>
#include <TTree.h>
#include <TChain.h>
#include <THStack.h>
#include <TCanvas.h>
#include <TF1.h>
#include <TLegend.h>
#include <stdio.h>
#include <iostream>
#include <deque>
#include <list>
#include <queue>
#include <fstream>
#include "EtaVEL.h"
#include "EtaVEL.cpp"
/*
Zum erstellen der correction map ist createGainAndEtaFile(...) in EVELAlg.C der entry point.
Zum erstellen des HR images ist createImage(...) der entry point.
*/
int etabins = 25;
int nEtas = 25;
Double_t dum[3][3];
Int_t x,y,f,q;
int counter[5];
int remoteCounter[5];
//TH2D *sum = new TH2D("sum","sum",3,-0.1,2.1,3,-0.1,2.1);
//TH2F *subPos = new TH2F("subPos","subPos", 100, -1.,1. ,100, -1.,1.);
TH2D *subPosAEta = new TH2D("subPosAEta","subPosAEta", 50, -.5,1.5 ,50, -.5,1.5);
TH2D *subPosBEta = new TH2D("subPosBEta","subPosBEta", 50, -.5,1.5 ,50, -.5,1.5);
TH1D *cE = new TH1D("clusterEnergy","clusterEnergy",400, 0.,4000.);
//TH1D *cES = new TH1D("clusterEnergyS","clusterEnergyS",400, 0.,4000.);
TH2D *cES3vs2 = new TH2D("clusterEnergy3vs2","clusterEnergy3vs2",800, 0.,8000.,600,0.,6000.);
TH2D *cES3vs2S = new TH2D("clusterEnergy3vs2S","clusterEnergy3vs2S",800, 0.,8000.,600,0.,6000.);
double th = 0.99;
double sigmas = 1.0;
TH2D *imgRLR = new TH2D("imgRLR","imgRLR",160,0.0,160.0 ,160 ,0.0,160.0);
TH2D *imgLR = new TH2D("imgLR","imgLR",160*2,0.0,160.0 ,160*2 ,0.0,160.0);
TH2D *clusHist= new TH2D("clusHist","clusHist",3,-0.5,2.5,3,-0.5,2.5);
TH2D *clusHistC= new TH2D("clusHistC","clusHistC",3,-0.5,2.5,3,-0.5,2.5);
int **imgArray;
int findShape(Double_t cluster[3][3], double sDum[2][2]){
int corner = -1;
double sum = cluster[0][0] + cluster[1][0] + cluster[2][0] + cluster[0][1] + cluster[1][1] + cluster[2][1] + cluster[0][2] + cluster[1][2] + cluster[2][2];
double sumTL = cluster[0][0] + cluster[1][0] + cluster[0][1] + cluster[1][1]; //2 ->BL
double sumTR = cluster[1][0] + cluster[2][0] + cluster[2][1] + cluster[1][1]; //0 ->TL
double sumBL = cluster[0][1] + cluster[0][2] + cluster[1][2] + cluster[1][1]; //3 ->BR
double sumBR = cluster[1][2] + cluster[2][1] + cluster[2][2] + cluster[1][1]; //1 ->TR
double sumMax = 0;
//double **sDum = subCluster;
Double_t ssDum[2][2];
// if(sumTL >= sumMax){
sDum[0][0] = cluster[0][0]; sDum[1][0] = cluster[1][0];
sDum[0][1] = cluster[0][1]; sDum[1][1] = cluster[1][1];
ssDum[0][0] = cluster[0][0]; ssDum[1][0] = cluster[0][1];
ssDum[0][1] = cluster[1][0]; ssDum[1][1] = cluster[1][1];
corner = 2;
sumMax=sumTL;
// }
if(sumTR >= sumMax){
sDum[0][0] = cluster[1][0]; sDum[1][0] = cluster[2][0];
sDum[0][1] = cluster[1][1]; sDum[1][1] = cluster[2][1];
ssDum[0][0] = cluster[2][0]; ssDum[1][0] = cluster[2][1];
ssDum[0][1] = cluster[1][0]; ssDum[1][1] = cluster[1][1];
corner = 0;
sumMax=sumTR;
}
if(sumBL >= sumMax){
sDum[0][0] = cluster[0][1]; sDum[1][0] = cluster[1][1];
sDum[0][1] = cluster[0][2]; sDum[1][1] = cluster[1][2];
ssDum[0][0] = cluster[0][2]; ssDum[1][0] = cluster[0][1];
ssDum[0][1] = cluster[1][2]; ssDum[1][1] = cluster[1][1];
corner = 3;
sumMax=sumBL;
}
if(sumBR >= sumMax){
sDum[0][0] = cluster[1][1]; sDum[1][0] = cluster[2][1];
sDum[0][1] = cluster[1][2]; sDum[1][1] = cluster[2][2];
ssDum[0][0] = cluster[2][2]; ssDum[1][0] = cluster[2][1];
ssDum[0][1] = cluster[1][2]; ssDum[1][1] = cluster[1][1];
corner = 1;
sumMax=sumBR;
}
switch(corner){
case 0:
cES3vs2->Fill(sum,sumTR); break;
case 1:
cES3vs2->Fill(sum,sumBR); break;
case 2:
cES3vs2->Fill(sum,sumTL); break;
case 3:
cES3vs2->Fill(sum,sumBL); break;
}
counter[corner]++;
remoteCounter[q]++;
// cout << "local corner is: " << corner << " remote corner is: " << q << endl;
return corner;
}
int placePhoton( TH2D *img, double subCluster[2][2], int cX, int cY, int corner, double *sX, double *sY, double *scX, double *scY){
double tot = subCluster[0][0] + subCluster[0][1] + subCluster[1][0] + subCluster[1][1];
double t = subCluster[1][0] + subCluster[1][1];
double r = subCluster[0][1] + subCluster[1][1];
double xHitC = r/tot;
double yHitC = t/tot;
imgRLR->Fill(cX,cY);
cE->Fill(tot);
double dX, dY;
//before looking at annas code
/* if(corner == 0){ dX=-1.; dY=-1.; }
if(corner == 1){ dX=-1.; dY=+1.; }
if(corner == 2){ dX=+1.; dY=-1.; }
if(corner == 3){ dX=+1.; dY=+1.; }*/
if(corner == 0){ dX=-1.; dY=+1.; } //top left
if(corner == 1){ dX=+1.; dY=+1.; } //top right
if(corner == 2){ dX=-1.; dY=-1.; } //bottom left
if(corner == 3){ dX=+1.; dY=-1.; } //bottom right
imgLR->Fill(cX+0.25*dX,cY+0.25*dY);
double posX = ((double)cX) + 0.5*dX + xHitC;
double posY = ((double)cY) + 0.5*dY + yHitC;
subPosBEta->Fill(xHitC ,yHitC);
if(img){
img->Fill(posX,posY);
}
if(xHitC < 0.02&& yHitC < 0.02){
cES3vs2S->Fill(dum[0][0]+dum[0][1]+dum[0][2]+dum[1][0]+dum[1][1]+dum[1][2]+dum[2][0]+dum[2][1]+dum[2][2],subCluster[0][0]+subCluster[0][1]+subCluster[1][0]+subCluster[1][1]);
}
if(sX && sY && scX && scY){
*sX = xHitC; //0.5 + 0.5*dX + xHitC;
*sY = yHitC; //0.5 + 0.5*dY + yHitC;
*scX = ((double)cX) + 0.5*dX;
*scY = ((double)cY) + 0.5*dY;
}
return 1;
}
void placePhotonCorr(TH2D *img, EtaVEL *e,double sX, double sY, double scX, double scY){
int bin = e->findBin(sX,sY);
if(bin <= 0) return;
double subX = ((double)(e->getXBin(bin))+.5)/((double)e->getNPixels());
double subY = ((double)(e->getYBin(bin))+.5)/((double)e->getNPixels());
if(img!=NULL){
img->Fill(scX+ subX , scY+ subY);
}
subPosAEta->Fill(subX,subY);
int iscx = scX;
int iscy = scY;
if(iscx >=nx || iscx<0 || iscy >=ny || iscy<0) return;
//cout << iscx*e->getNPixels()+e->getXBin(bin) << " " << iscy*e->getNPixels()+e->getXBin(bin) << endl;
if(img==NULL) return;
imgArray[iscx*e->getNPixels()+e->getXBin(bin)][iscy*e->getNPixels()+e->getYBin(bin)]++;
}
void gainCorrection(Double_t corrected[3][3], TH2D *gainMap){
for(int xx = 0; xx < 3; xx++)
for(int yy = 0; yy < 3; yy++){
if(gainMap && gainMap->GetBinContent(x+xx+2,y+yy+2) != 0){
corrected[xx][yy] = dum[xx][yy] / gainMap->GetBinContent(x+xx+2,y+yy+2);
clusHistC->Fill(xx,yy,corrected[xx][yy]);
}
else
corrected[xx][yy] = dum[xx][yy];
clusHist->Fill(xx,yy,dum[xx][yy]);
}
}
EtaVEL *plotEtaDensity(TChain* tree2, TEntryList *el, EtaVEL *oldEta = NULL, TH2D **img = NULL, TH2D *gainMap=NULL, int nPixels=25) {
EtaVEL *newEta = new EtaVEL(25,-0.02,1.02);
Long64_t listEntries=el->GetN();
Long64_t treeEntry;
Long64_t chainEntry;
Int_t treenum=0;
tree2->SetEntryList(el);
double gainCorrC[3][3];
double subCluster[2][2];
double sX, sY, scX, scY;
cout << "Events: " << listEntries << endl;
if(oldEta == NULL){ cout << "Old Eta is NULL " << endl; }
for(int i = 0; i<4; i++){ counter[i] = 0; remoteCounter[i] = 0; }
for (Long64_t il =0; il<listEntries;il++) {
treeEntry = el->GetEntryAndTree(il,treenum);
chainEntry = treeEntry+tree2->GetTreeOffset()[treenum];
if (tree2->GetEntry(chainEntry)) {
gainCorrection(gainCorrC,gainMap);
//cout << gainCorrC[1][1] << endl;
//finds corner
int corner = findShape(gainCorrC,subCluster);
int validEvent;
if(img){
validEvent = placePhoton(img[0],subCluster,x,y, corner, &sX, &sY, &scX, &scY);
}else{
//calc etaX, etaY
validEvent = placePhoton(NULL,subCluster,x,y, corner, &sX, &sY, &scX, &scY);
}
//fill etavel
newEta->fill(sX,sY);
if(oldEta && img && img[1]){
placePhotonCorr(img[1],oldEta, sX,sY, scX, scY);
}else{
placePhotonCorr(NULL,newEta,sX,sY,scX,scY);
}
}
//cout << il << endl;
int ssize = 500000;
if(il % ssize == 0 && il != 0 && oldEta==NULL){
cout << " -------------- "<< endl;
newEta->updatePixelPos();
//newEta->resolveSelfIntersect();
char tit[1000];
/* TFile *ff = new TFile("/scratch/Spider.root","UPDATE");
sprintf(tit,"subPosAEta%i",newEta->getIt()); subPosAEta->SetName(tit);
subPosAEta->Write(); subPosAEta->Reset();
sprintf(tit,"subPosBEta%i",newEta->getIt()); subPosBEta->SetName(tit);
subPosBEta->Write(); subPosBEta->Reset();
sprintf(tit,"Eta%i",newEta->getIt()); newEta->Write(tit);
ff->Close(); */
//il = 0;
}
if(il % ssize == ssize-1){
double prog = (double)il/(double)listEntries*100.;
cout << prog << "%" << endl;
//if(prog > 19.) return newEta;
if(newEta->converged == 1){ cout << "converged ... " << endl; return newEta; }
}
}
cout << "local corners: " ;
for(int i = 0; i<4; i++) cout << i << ": " << counter[i] << " || " ;
cout << endl;
//cout << "remote corners: " ;
//for(int i = 0; i<4; i++) cout << i << ": " << remoteCounter[i] << " || " ;
//cout << endl;
return newEta;
}
TChain *openTree(char *tname, char *fname,double lEc, double hEc, double rms=5., char *chainName=">>thischan"){
TChain *tree2;
// TH1D **etaDI;
char cut[1000];
tree2=new TChain(tname);
tree2->Add(fname);
tree2->Print();
//sprintf(cut,"(x<=40) && (data[%d][%d]>%f*rms) && Sum$(data)<%f && Sum$(data)>%f",1,1,rms, hEc, lEc);
// sprintf(cut,"(x<=40) && (data[%d][%d]>%f*rms)",1,1,rms);// && Sum$(data)<%f && Sum$(data)>%f",1,1,rms, hEc, lEc);
sprintf(cut,"(x<=40) && Sum$(data)<%f && Sum$(data)>%f", hEc, lEc);
// sprintf(cut,"");
cout << cut << endl;
tree2->Draw(chainName, cut, "entrylist");
tree2->SetBranchAddress("iFrame",&f);
tree2->SetBranchAddress("x",&x);
tree2->SetBranchAddress("y",&y);
tree2->SetBranchAddress("data",dum);
//tree2->SetBranchAddress("q",&q);
cout << "openTree : end" << endl;
return tree2;
}
EtaVEL *etaDensity(char *tname, char *fname, double lEc = 1000, double hEc=3000, TH2D *gainMap=NULL, int nPixels=25) {
/** open tree and make selection */
TChain *tree2 = openTree(tname,fname,lEc,hEc);
TEntryList *elist = (TEntryList*)gDirectory->Get("thischan");
if(elist == NULL) { cout << "could not open tree " << endl; return NULL; }
EtaVEL *etaDen = plotEtaDensity(tree2,elist,NULL,NULL,gainMap,nPixels);
//etaDen->Draw("colz");
cout << "done" << endl;
return etaDen;
}
void interpolate(char *tname, char *fname, EtaVEL *etaDI, double lEc = 1000, double hEc=3000, TH2D *gainMap=NULL) {
TChain *tree2 = openTree(tname,fname,lEc,hEc,5.,">>intChain");
TEntryList *elist = (TEntryList*)gDirectory->Get("intChain");
if(elist == NULL) { cout << "could not open tree " << endl; return; }
double nPixels = (double)etaDI->getNPixels();
TH2D **img = new TH2D*[3];
img[0] = new TH2D("img","img",nPixels*160,0.0,160.0 ,nPixels*160 ,0.0,160.0);
img[1] = new TH2D("imgE","imgE",nPixels*160,0.0,160.0 ,nPixels*160 ,0.0,160.0);
int inPixels = etaDI->getNPixels();
imgArray = new int*[inPixels*160];
for(int i = 0; i < inPixels*160; i++){
imgArray[i] = new int[inPixels*160];
for(int j = 0; j < inPixels*160; j++){
imgArray[i][j] = 0;
}
}
cout << "starting" << endl;
plotEtaDensity(tree2,elist, etaDI,img,gainMap);
//img->Draw("colz");
}
TH2D *createGainMap(char *tname, char *fname, double lEc = 0,double hEc=10000){
char name[100];
TH1D *avgSpec3 = new TH1D("avgSpec3", "avgSpec3",hEc/20,0,hEc);
TH1D ***specs3 = new TH1D**[160];
TH1D ***specs1 = new TH1D**[160];
for(int xx = 0; xx < 160; xx++){
specs3[xx] = new TH1D*[160];
specs1[xx] = new TH1D*[160];
for(int yy = 0; yy < 160; yy++){
sprintf(name,"S3x%iy%i",xx,yy);
specs3[xx][yy] = new TH1D(name,name,hEc/20,0,hEc);
sprintf(name,"S1x%iy%i",xx,yy);
specs1[xx][yy] = new TH1D(name,name,hEc/20,0,hEc);
}
}
TChain *tree2 = openTree(tname,fname,0,hEc,5.,">>gainChan");
TEntryList *elist = (TEntryList*)gDirectory->Get("gainChan");
if(elist == NULL) { cout << "could not open tree " << endl; return NULL; }
Long64_t listEntries=elist->GetN();
Long64_t treeEntry;
Long64_t chainEntry;
Int_t treenum=0;
tree2->SetEntryList(elist);
cout << "Events: " << listEntries << endl;
for(int i = 0; i<4; i++) counter[i] = 0;
for (Long64_t il =0; il<listEntries;il++) {
treeEntry = elist->GetEntryAndTree(il,treenum);
chainEntry = treeEntry+tree2->GetTreeOffset()[treenum];
if (tree2->GetEntry(chainEntry)) {
double sum = 0;
for(int xx = 0; xx < 3; xx++)
for(int yy = 0; yy < 3; yy++)
sum += dum[xx][yy];
specs3[x][y]->Fill(sum);
specs1[x][y]->Fill(dum[1][1]);
avgSpec3->Fill(sum);
}
}
TH2D *gainMap3 = new TH2D("gainMap3","gainMap3",160,-0.5,160.-0.5,160,-.5,160.-.5);
TH2D *gainMap1 = new TH2D("gainMap1","gainMap1",160,-0.5,160.-0.5,160,-.5,160.-.5);
for(int xx = 0; xx < 160; xx++){
for(int yy = 0; yy < 160; yy++){
TF1 *gf3 = new TF1("gf3","gaus", lEc, hEc);
specs3[xx][yy]->Fit(gf3,"Q");
double e3 = gf3->GetParameter(1);
gainMap3->Fill(xx,yy,e3);
TF1 *gf1 = new TF1("gf1","gaus", lEc, hEc);
specs1[xx][yy]->Fit(gf1,"Q");
double e1 = gf1->GetParameter(1);
gainMap1->Fill(xx,yy,e1);
}
}
return gainMap3;
}
void writeMatlab2DHisto(int xx, int yy,char *outFileName){
ofstream outFile;
outFile.open (outFileName);
cout << "create matlab file with " << xx << " xbins and " << yy << " ybins" << endl;
for(int y = 0; y < yy; y++){
for(int x = 0; x < xx; x++){
outFile << imgArray[x][y] << "\t";
}
outFile << endl;
}
outFile.close();
}
//COMPLETE STUFF
void createImage(char *tdir, char *tname, char *ftname, char *ifname = NULL, int useGM=0, double lEth=-1., double hEth=-1.){
imgRLR->Reset();
imgLR->Reset();
char fname[1000];
char inFName[1000];
char outFName[1000];
char moutFName[1000];
if(ifname == NULL){
sprintf(fname,"%s/%s_*.root",tdir,tname);
}else{
sprintf(fname,"%s",ifname);
}
if(useGM) sprintf(inFName,"%s/%s-PlotsWGMVEL.root",tdir,ftname);
else sprintf(inFName,"%s/%s-PlotsVEL.root",tdir,ftname);
sprintf(outFName,"%s/%s-ImgVEL.root",tdir,tname);
sprintf(moutFName,"%s/%s-ImgVEL.mf",tdir,tname);
TFile *inFile = new TFile(inFName,"READ");
cout << "Image Tree File Name: " << fname << endl;
cout << "Eta File Name: " << inFName << endl;
cout << "Out File Name: " << outFName << endl;
cout << "Matlab Out File Name: " << moutFName << endl;
TH2D *gm = NULL;
if(useGM){
cout << "Load gain map" << endl;
gm = (TH2D *)gDirectory->Get("gainMap");
if(gm == NULL){ cout << "can not find gainMap in file" << endl; return; }
}
cout << "Load eta" << endl;
EtaVEL *ee = (EtaVEL *)gDirectory->Get("etaDist");
cout << "Select Energy BW" << endl;
TH1D *spec = (TH1D *)gDirectory->Get("avgSpec3");
if(spec == NULL){ cout << "can not find avgSpec3" << endl; return; }
TF1 *gf3 = new TF1("gf3","gaus", 0, 10000);
spec->Fit(gf3,"Q");
double avgE = gf3->GetParameter(1);
double sigE = gf3->GetParameter(2);
cout << "avgE: " << avgE << " sigE: " << sigE << endl;
cout << endl;
if(lEth == -1.) lEth = avgE-5.*sigE;
if(hEth == -1.) hEth = avgE+5.*sigE;
cout << lEth << " < E < " << hEth << " (eV)" << endl;
cout << "start with interpolation" << endl;
interpolate( tname, fname, ee,lEth,hEth ,gm);
TH2D *img = (TH2D *)gDirectory->Get("img");
if(img == NULL){ cout << "could not find 2d-histogram: img " << endl; return; }
TH2D *imgE = (TH2D *)gDirectory->Get("imgE");
if(imgE == NULL){ cout << "could not find 2d-histogram: imgE " << endl; return; }
//TH2D *imgEOM = (TH2D *)gDirectory->Get("imgEOM");
//if(imgEOM == NULL){ cout << "could not find 2d-histogram: imgEOM " << endl; return; }
TFile *outFile = new TFile(outFName,"UPDATE");
imgLR->Write();
imgRLR->Write();
imgE->Write();
//imgEOM->Write();
img->Write();
outFile->Close();
inFile->Close();
cout << "writing matlab file: " << moutFName << endl;
writeMatlab2DHisto(160*ee->getNPixels(),160*ee->getNPixels(),moutFName);
cout << "Done : " << outFName << endl;
}
/**
\par tdir input tree directory
\par tname input tree name
\par ifname input file name if different than tdir/tname_*.root
\par useGM use gain map
\par maxExpEinEv spectrum maximum
\par nPixels sub-pixels bins
\par lEth low threshold
\par hEth high threshold
*/
EtaVEL *createGainAndEtaFile(char *tdir, char *tname, char *ifname=NULL, int useGM=0, double maxExpEinEv=25000., int nPixels =25, double lEth=-1., double hEth=-1.){
char fname[1000];
char outFName[1000];
if(ifname == NULL){
sprintf(fname,"%s/%s_*.root",tdir,tname);
}else{
sprintf(fname,"%s",ifname);
}
if(useGM) sprintf(outFName,"%s/%s-PlotsWGVEL.root",tdir,tname);
else sprintf(outFName,"%s/%s-PlotsVEL.root",tdir,tname);
cout << "Tree File Name: " << fname << endl;
cout << "Output File Name: " << outFName << endl;
/** creates gain map and 3x3 spectrum */
cout << "Creating gain map: " << endl;
TH2D *gm = createGainMap(tname,fname,0,maxExpEinEv/10.);
gm->SetName("gainMap");
/** gets average 3x3 spectrum and fits it with a gaus */
TH1D *spec = (TH1D *)gDirectory->Get("avgSpec3");
if(spec == NULL){ cout << "can not find avgSpec3" << endl; return NULL; }
TF1 *gf3 = new TF1("gf3","gaus", 0, maxExpEinEv/10.);
spec->Fit(gf3,"Q");
double avgE = gf3->GetParameter(1);
double sigE = gf3->GetParameter(2);
cout << "avgE: " << avgE << " sigE: " << sigE << endl;
cout << endl;
/** sets high and low threshold if not given by the user */
if(lEth == -1.) lEth = avgE-5.*sigE;
if(hEth == -1.) hEth = avgE+5.*sigE;
cout << lEth << " < E < " << hEth << " (eV)" << endl;
cout << "calculating eta stuff" << endl;
EtaVEL *newEta;
if(useGM) newEta = etaDensity(tname,fname,lEth,hEth,gm,nPixels);
else newEta = etaDensity(tname,fname,lEth,hEth,NULL,nPixels);
cout << "writing to file " << outFName << endl;
TFile *outFile = new TFile(outFName,"UPDATE");
newEta->Write("etaDist");
gm->Write();
spec->Write();
subPosAEta->Write();
cES3vs2->Write();
outFile->Close();
cout << "Done : " << outFName << endl;
return newEta;
}
void exportSpec(char *tdir, char *tname){
char tfname[1000];
char ofname[1000];
char cleanName[1000];
for(int p = 0; p < strlen(tname);p++){
cleanName[p+1] = '\0';
cleanName[p] = tname[p];
if(tname[p] == '-') cleanName[p] = '_';
}
sprintf(tfname,"%s/%s-PlotsVEL.root",tdir,tname);
sprintf(ofname,"%s/%s_SpecVEL.m",tdir,cleanName);
TFile *tf = new TFile(tfname);
TH1D *spec = (TH1D *)gDirectory->Get("avgSpec3");
ofstream outFile;
outFile.open (ofname);
if(outFile.fail()){
cout << "Could not open file : " << ofname << endl;
return;
}
cout << "create matlab file with with spec " << ofname << endl;
outFile << cleanName << " = [ " << endl;
for(int i = 0; i < spec->GetNbinsX(); i++){
outFile << i << " " << spec->GetBinCenter(i) << " " << spec->GetBinContent(i) << " ; " << endl;
}
outFile << " ] ; " << endl;
outFile.close();
}

679
slsDetectorCalibration/interpolations/etaVEL/EtaVEL.cpp
View File

@ -1,679 +0,0 @@
#include "EtaVEL.h"
#include <iomanip>
// ClassImp(EtaVEL);
// double Median(const TH1D * histo1) {
// int numBins = histo1->GetXaxis()->GetNbins();
// Double_t *x = new Double_t[numBins];
// Double_t* y = new Double_t[numBins];
// for (int i = 0; i < numBins; i++) {
// x[i] = histo1->GetBinCenter(i);
// y[i] = histo1->GetBinContent(i);
// }
// return TMath::Median(numBins, x, y);
// }
double *EtaVEL::getPixelCorners(int x, int y){
double tlX,tlY,trX,trY,blX,blY,brX,brY;
tlX = xPPos[getCorner(x,y+1)];
tlY = yPPos[getCorner(x,y+1)];
trX = xPPos[getCorner(x+1,y+1)];
trY = yPPos[getCorner(x+1,y+1)];
blX = xPPos[getCorner(x,y)];
blY = yPPos[getCorner(x,y)];
brX = xPPos[getCorner(x+1,y)];
brY = yPPos[getCorner(x+1,y)];
//cout << "gPC: TL: " << getCorner(x,y+1) << " TR: " << getCorner(x+1,y+1) << " BL " << getCorner(x,y) << " BR " << getCorner(x+1,y) << endl;
double *c = new double[8];
c[0] = tlX; c[1] = trX; c[2] = brX; c[3] = blX;
c[4] = tlY; c[5] = trY; c[6] = brY; c[7] = blY;
return c;
}
int EtaVEL::findBin(double xx, double yy){
double tlX,tlY,trX,trY,blX,blY,brX,brY;
/********Added by anna ******/
// if (xx<min) xx=min+1E-6;
// if (xx>max) xx=max-1E-6;
// if (yy<min) yy=min+1E-6;
// if (yy>max) yy=max-1E-6;
/**************/
int bin = -1;
for(int x = 0; x < nPixels; x++){
for(int y = 0; y < nPixels; y++){
double *c = getPixelCorners(x,y);
tlX = c[0]; trX = c[1]; brX = c[2]; blX = c[3];
tlY = c[4]; trY = c[5]; brY = c[6]; blY = c[7];
///if(y == 0){
// cout << "x: " << x << " blY " << blY << " brY " << brY << endl;
//}
int out = 0;
double tb = 0;
double bb = 0;
double lb = 0;
double rb = 0;
if((trX-tlX)>0.)
tb = (trY - tlY)/(trX-tlX);
if((brX-blX)>0.)
bb = (brY - blY)/(brX-blX);
if((tlY-blY)>0.)
lb = (tlX - blX)/(tlY-blY);
if((trY-brY)>0.)
rb = (trX - brX)/(trY-brY);
double ty = tlY + tb * (xx - tlX);
double by = blY + bb * (xx - blX);
double lx = blX + lb * (yy - blY);
double rx = brX + rb * (yy - brY);
if(yy >= ty) out++;
if(yy < by) out++;
if(xx < lx) out++;
if(xx >= rx) out++;
//cout << "ty " << ty << endl;
//cout << "by " << by << endl;
//cout << "lx " << lx << endl;
//cout << "rx " << rx << endl;
//double dist = (xx - xPPos[getBin(x,y)]) * (xx - xPPos[getBin(x,y)]) + (yy - yPPos[getBin(x,y)]) * (yy - yPPos[getBin(x,y)]);
//cout << "x " << x << " y " << y << " out " << out << " ty " << ty << endl;
//cout << "tl " << tlX << "/" << tlY << " tr " << trX << "/" << trY << endl;
//cout << "bl " << blX << "/" << blY << " br " << brX << "/" << brY << endl;
//cout << " tb " << tb << endl;
delete[] c;
if(out == 0){ return getBin(x,y); }
}
}
return -1;
}
void EtaVEL::createLogEntry(){
if(it >= nIterations){
cerr << "log full" << endl;
}
log[it].itN = it;
log[it].xPos = new double[nPixels*nPixels+1];
log[it].yPos = new double[nPixels*nPixels+1];
log[it].binCont = new double[nPixels*nPixels+1];
for(int x = 0; x < nPixels; x++)
for(int y = 0; y < nPixels; y++){
log[it].xPos[getBin(x,y)] = xPPos[getBin(x,y)];
log[it].yPos[getBin(x,y)] = yPPos[getBin(x,y)];
log[it].binCont[getBin(x,y)] = binCont[getBin(x,y)];
}
it++;
}
void EtaVEL::updatePixelCorner(){
double w = 20;
int rows = (nPixels+1)*(nPixels+1) + 4 + 4 * 4;//(4*(nPixels+1))-4;
int cols = (nPixels+1)*(nPixels+1);
double *rVx = new double[rows];
double *rVy = new double[rows];
double *posMat = new double[rows*cols];
for(int i = 0 ; i < rows*cols; i++) posMat[i] = 0;
int boundaryPoint = 0;
cout << "linear sys stuff" << endl;
double minELength = 100000000000000; int minX=-1, minY=-1;
for(int y = 0; y < nPixels+1; y++){
for(int x = 0; x < nPixels+1; x++){
double bx = 0, by = 0;
//boundary conditions
if((x == 0 && y % 5 == 0) ||
(x == nPixels && y % 5 == 0) ||
(y == 0 && x % 5 == 0) ||
(y == nPixels && x % 5 == 0)){
bx = xPPos[getCorner(x,y)];
//cout << "bP " << boundaryPoint << " bx " << bx << endl;
by = yPPos[getCorner(x,y)];
rVx[(nPixels+1)*(nPixels+1) + boundaryPoint] = bx*w;
rVy[(nPixels+1)*(nPixels+1) + boundaryPoint] = by*w;
posMat[(nPixels+1)*(nPixels+1)*cols + boundaryPoint * cols + getCorner(x,y)-1] = w;
boundaryPoint++;
}
double tot = 4 - (x == 0) - (y == 0) - (x == nPixels) - (y == nPixels);
//cout << "totW: " << tot << endl;
//tot = 4.;
double eLength = 0;
if(x != 0) eLength += edgeL[getEdgeX(x-1,y)];
if(y != 0) eLength += edgeL[getEdgeY(x,y-1)];
if(x != nPixels) eLength += edgeL[getEdgeX(x,y)];
if(y != nPixels) eLength += edgeL[getEdgeY(x,y)];
/*cout << "Corner X:" <<x << " Y: " << y ;
cout << " C# " << getCorner(x,y);
cout << " eXl " << getEdgeX(x-1,y) << "(C# " << getCorner(x-1,y) << " ) ";
cout << " eXr " << getEdgeX(x,y) << "(C# " << getCorner(x+1,y) << " ) ";
cout << " eYb " << getEdgeY(x,y-1) << "(C# " << getCorner(x,y-1) << " ) ";
cout << " eYt " << getEdgeY(x,y) << "(C# " << getCorner(x,y+1) << " ) " << endl; */
//" totW: " << tot << " totE: " << eLength << endl;
if(eLength < minELength & tot == 4){
minELength = eLength;
minX = x; minY = y;
}
//matrixes updated
if(x != 0) posMat[y*(nPixels+1)*cols+x*cols+getCorner(x-1,y)-1] = -edgeL[getEdgeX(x-1,y)]/eLength;
if(y != 0) posMat[y*(nPixels+1)*cols+x*cols+getCorner(x,y-1)-1] = -edgeL[getEdgeY(x,y-1)]/eLength;;
if(x != nPixels) posMat[y*(nPixels+1)*cols+x*cols+getCorner(x+1,y)-1] = -edgeL[getEdgeX(x,y)]/eLength;;
if(y != nPixels) posMat[y*(nPixels+1)*cols+x*cols+getCorner(x,y+1)-1] = -edgeL[getEdgeY(x,y)]/eLength;;
posMat[y*(nPixels+1)*cols+x*cols+getCorner(x,y)-1] = 1.;
rVx[getCorner(x,y)-1] = 0.;
rVy[getCorner(x,y)-1] = 0.;
}
}
cout << "Min Corner X: " <<minX << " Y: " << minY << " C# " << getCorner(minX,minY) << " length " << minELength << endl;
TMatrixD *k = new TMatrixD(rows,cols);
TVectorD *fx = new TVectorD(rows,rVx);
TVectorD *fy = new TVectorD(rows,rVy);
// f->Print();
k->SetMatrixArray(posMat);
// k->Print();
//solve linear system
Bool_t ok;
TDecompSVD *s = new TDecompSVD(*k);
s->Solve(*fx);
s->Solve(*fy);
double *fxA = fx->GetMatrixArray();
double *fyA = fy->GetMatrixArray();
for(int y = 0; y < nPixels+1; y++){
for(int x = 0; x < nPixels+1; x++){
//do not update boundaries
if(!(x == 0 ||
x == nPixels||
y == 0 ||
y == nPixels)){
xPPos[getCorner(x,y)] = fxA[getCorner(x,y)-1];
yPPos[getCorner(x,y)] = fyA[getCorner(x,y)-1];
}
}
}
}
void EtaVEL::updatePixelPos(){
double xMov, yMov, d1Mov, d2Mov;
createLogEntry();
double *chMap = getChangeMap();
int ch =0;
cout << "update edge lengths" << endl;
for(int x = 0; x < nPixels; x++)
for(int y = 0; y < nPixels; y++){
/*cout << "Pixel X:" <<x << " Y: " << y << " P# " << getBin(x,y) << " eXb " << getEdgeX(x,y);
cout << " eXt " << getEdgeX(x,y+1) << " eYl " << getEdgeY(x,y) << " eYr " << getEdgeY(x+1,y) << endl;
*/
edgeL[getEdgeX(x,y)] *= chMap[getBin(x,y)];
edgeL[getEdgeX(x,y+1)] *= chMap[getBin(x,y)];
edgeL[getEdgeY(x,y)] *= chMap[getBin(x,y)];
edgeL[getEdgeY(x+1,y)] *= chMap[getBin(x,y)];
//cout << "Pixel x: " << x << " y: " << y << " Ch: " << chMap[getBin(x,y)] << " counts: " << binCont[getBin(x,y)] << endl;
//cout << "BE " << getEdgeX(x,y) << endl;
//cout << "TE " << getEdgeX(x,y+1) << endl;
//cout << "LE " << getEdgeY(x,y) << endl;
//cout << "RE " << getEdgeY(x+1,y) << endl;
binCont[getBin(x,y)] = 0;
}
updatePixelCorner();
//double *pSize = getSizeMap();
double totEdgeLength = 0;
for(int e = 1; e < 2*nPixels*(nPixels+1)+1; e++){
totEdgeLength += edgeL[e];
}
cout << "tot edge Length: " << totEdgeLength << endl;
totCont = 0.;
}
double *EtaVEL::getSizeMap(){
double tlX,tlY,trX,trY,blX,blY,brX,brY;
double *szMap = new double[nPixels*nPixels+1];
for(int x = 1; x < nPixels-1; x++)
for(int y = 1; y < nPixels-1; y++){
double *c = getPixelCorners(x,y);
tlX = c[0]; trX = c[1]; brX = c[2]; blX = c[3];
tlY = c[4]; trY = c[5]; brY = c[6]; blY = c[7];
//double area = dtl * dtr / 2. + dtr * dbr / 2. + dbr * dbl / 2. + dbl * dtl / 2.;
//http://en.wikipedia.org/wiki/Shoelace_formula
double sl1 = tlX * trY + trX * brY + brX * blY + blX * tlY;
double sl2 = tlY * trX + trY * brX + brY * blX + blY * tlX;
double area = 1./2. * (- sl1 + sl2);
if(area < 0.){
cout << "negative area: X " << x << " Y " << y << " area " << endl;
edgeL[getEdgeX(x,y)] *= 2.;
edgeL[getEdgeX(x,y+1)] *= 2.;
edgeL[getEdgeY(x,y)] *= 2.;
edgeL[getEdgeY(x+1,y)] *= 2.;
}
szMap[getBin(x,y)] = area / (max - min) / (max - min) * nPixels * nPixels;
delete[] c;
}
return szMap;
}
double *EtaVEL::getChangeMap(){
double *chMap = new double[nPixels*nPixels+1];
double avg = totCont/(double)(nPixels*nPixels);
// TH1D *hmed=getCounts();
// double med = Median(hmed);
// delete hmed;
double acc = TMath::Sqrt(avg);
cout << "totC: " << totCont << " avg " << avg << " acc: " << acc << endl;//<< " med " << med
double totOffAcc = 0.;
int totInRange03s = 0;
int totInRange07s = 0;
int totInRange12s = 0;
int totInRange20s = 0;
int totInRange25s = 0;
double dd;
int totInBins = 0;
//double
chi_sq=0;
int maxC = 0, maxX=-1, maxY=-1;
double minC = 1000000000000000, minX, minY;
for(int x = 0; x < nPixels; x++){
for(int y = 0; y < nPixels; y++){
totInBins += binCont[getBin(x,y)];
double r = (double)binCont[getBin(x,y)];
if(r > 0. & totCont > 0.){
dd=sqrt(r/avg);
/**Added by Anna */
if (dd>2.) dd=1.5;
if (dd<0.5) dd=0.75;
chMap[getBin(x,y)] = dd;
/** */
//if( chMap[getBin(x,y)] < 1.){ chMap[getBin(x,y)] = 1/1.2; }
//if( chMap[getBin(x,y)] > 1.){ chMap[getBin(x,y)] = 1.2; }
//if( chMap[getBin(x,y)] < 1/1.2){ chMap[getBin(x,y)] = 1/1.2; }
//if( chMap[getBin(x,y)] > 1.2){ chMap[getBin(x,y)] = 1.2; }
}else if(totCont > 0.){
chMap[getBin(x,y)] =0.5; //1/1.2;
}else{
chMap[getBin(x,y)] = 1.;
}
//if(r < avg + 2*acc && r > avg - 2*acc){ totInRange++;}// chMap[getBin(x,y)] = 1.; }
/** Commente away by Anna
if(converged == 0 && r < med+20*acc){ chMap[getBin(x,y)] = 1.; }
if(converged == 2 && r < med+20*acc && r > med-03*acc){ chMap[getBin(x,y)] = 1.; }
if(r < med+03*acc){ totInRange03s++; }
if(r < med+07*acc){ totInRange07s++; }
if(r < med+12*acc){ totInRange12s++; }
if(r < med+20*acc){ totInRange20s++; }
if(r < med+25*acc){ totInRange25s++; }
*/
//cout << "x " << x << " y " << y << " r " << r << " ch " << chMap[getBin(x,y)] << endl;
// if(r - avg > acc){ totOffAcc += r-avg;}
//if(r - avg < -acc){ totOffAcc += avg-r;}
totOffAcc += (avg-r)*(avg-r);
chi_sq+=(avg-r)*(avg-r)/r;
//cout << " x " << x << " y " << y << " bC " << binCont[x*nPixels+y] << " r " << r << endl;
if(r > maxC){ maxC = r; maxX = x; maxY = y; }
if(r < minC){minC = r; minX = x; minY = y; }
}
}
// cout << "totInBins " << totInBins << " zero Bin " << binCont[0] << endl;
cout << "AvgOffAcc: " << sqrt(totOffAcc/(double)(nPixels*nPixels)) << endl;
cout << "***********Reduced Chi Square: " << chi_sq/((double)(nPixels*nPixels)) << endl;
// cout << "totInRange03 (<" << med+03*acc << "): " << totInRange03s << endl;
// cout << "totInRange07 (<" << med+07*acc << "): " << totInRange07s << endl;
// cout << "totInRange12 (<" << med+12*acc << "): " << totInRange12s << endl;
// cout << "totInRange20 (<" << med+20*acc << "): " << totInRange20s << endl;
// cout << "totInRange25 (<" << med+25*acc << "): " << totInRange25s << endl;
double maxSig = (maxC - avg)*(maxC - avg) / avg;//acc;
double minSig = (avg - minC)*(avg - minC) / avg;//acc;
cout << "Max Pixel X: " << maxX << " Y: " << maxY << " P# " << getBin(maxX,maxY) << " count: " << maxC << " sig: "<< maxSig << endl;
cout << "Min Pixel X: " << minX << " Y: " << minY << " P# " << getBin(minX,minY) << " count: " << minC << " sig: "<< minSig << endl;
// if(maxSig <= 25){ converged = 2; cout << "reached first converstion step!!!" << endl; }
//if(minSig <= 7 && converged == 2) { converged = 1; }
if (chi_sq<3) converged=2;
if (chi_sq<1) converged=1;
cout << "Conversion step "<< converged << endl;
return chMap;
}
TH2D *EtaVEL::getContent(int it, int changeType){
TH2D *cont = new TH2D("cont","cont",nPixels,min,max,nPixels,min,max);
double *chMap = NULL;
if(changeType ==1) chMap = getChangeMap();
double *szMap = getSizeMap();
for(int x = 0; x < nPixels; x++)
for(int y = 0; y < nPixels; y++){
if(changeType ==2 ){
cont->SetBinContent(x+1,y+1,szMap[getBin(x,y)]);
}
if(changeType ==1 ){
cont->SetBinContent(x+1,y+1,chMap[getBin(x,y)]);
}
if(changeType ==0 ){
if(it == -1){
cont->SetBinContent(x+1,y+1,binCont[getBin(x,y)]);
//cout << "x " << x << " y " << y << " cont " << binCont[getBin(x,y)] << endl;
}
else{cont->SetBinContent(x+1,y+1,log[it].binCont[getBin(x,y)]);}
}
}
return cont;
}
TH1D *EtaVEL::getCounts(){
TH1D *ch = new TH1D("ch","ch",500,0,totCont/(nPixels*nPixels)*4);
for(int x = 0; x < nPixels; x++)
for(int y = 0; y < nPixels; y++){
ch->Fill(binCont[getBin(x,y)]);
}
return ch;
}
void EtaVEL::printGrid(){
double *colSum = new double[nPixels+1];
double *rowSum = new double[nPixels+1];
for(int i = 0; i < nPixels+1; i++){
colSum[i] = 0.;
rowSum[i] = 0.;
for(int j = 0; j < nPixels; j++){
rowSum[i] += edgeL[getEdgeX(j,i)];
colSum[i] += edgeL[getEdgeY(i,j)];
}
}
cout << endl;
cout.precision(3); cout << fixed;
cout << " ";
for(int x = 0; x < nPixels+1; x++){
cout << setw(2) << x << " (" << colSum[x] << ") ";
}
cout << endl;
for(int y = 0; y < nPixels+1; y++){
cout << setw(2) << y << " ";
for(int x = 0; x < nPixels+1; x++){
cout << "(" << xPPos[getCorner(x,y)] << "/" << yPPos[getCorner(x,y)] << ") " ;
if(x < nPixels) cout << " -- " << edgeL[getEdgeX(x,y)]/rowSum[y]*(max-min) << " -- ";
}
cout << " | " << rowSum[y] << endl;
if(y < nPixels){
cout << " ";
for(int x = 0; x < nPixels+1; x++){
cout << edgeL[getEdgeY(x,y)]/colSum[x]*(max-min) << " ";
}
cout << endl;
}
}
delete[] colSum;
delete[] rowSum;
}
TMultiGraph *EtaVEL::plotPixelBorder(int plotCenters){
TMultiGraph *mg = new TMultiGraph();
double cx[5], cy[5];
for(int x = 0; x < nPixels; x++)
for(int y = 0; y < nPixels; y++){
double *c = getPixelCorners(x,y);
cx[0]=c[0]; cx[1]=c[1]; cx[2]=c[2]; cx[3]=c[3]; cx[4]=c[0];
cy[0]=c[4]; cy[1]=c[5]; cy[2]=c[6]; cy[3]=c[7]; cy[4]=c[4];
TGraph *g = new TGraph(5,cx,cy);
mg->Add(g);
if(plotCenters){
g = new TGraph(1,&(xPPos[getBin(x,y)]),&(yPPos[getBin(x,y)]));
mg->Add(g);
}
delete[] c;
}
return mg;
}
TMultiGraph *EtaVEL::plotLog(int stepSize, int maxIt){
int mIt;
TMultiGraph *mg = new TMultiGraph();
double **xposl = new double*[nPixels*nPixels+1];
double **yposl = new double*[nPixels*nPixels+1];
if(maxIt==-1){ mIt = it; } else{ mIt = maxIt; };
cout << "mIt " << mIt << " steps " << mIt/stepSize << endl;
for(int x = 0; x < nPixels; x++){
for(int y = 0; y < nPixels; y++){
xposl[getBin(x,y)] = new double[mIt/stepSize];
yposl[getBin(x,y)] = new double[mIt/stepSize];
for(int i = 0; i < mIt/stepSize; i++){
xposl[getBin(x,y)][i] = log[i*stepSize].xPos[getBin(x,y)];
yposl[getBin(x,y)][i] = log[i*stepSize].yPos[getBin(x,y)];
}
TGraph *g = new TGraph(mIt/stepSize,xposl[getBin(x,y)],yposl[getBin(x,y)]);
g->SetLineColor((x*y % 9) + 1);
if(x == 0) g->SetLineColor(2);
if(y == 0) g->SetLineColor(3);
if(x == nPixels-1) g->SetLineColor(4);
if(y == nPixels-1) g->SetLineColor(5);
mg->Add(g);
}
}
return mg;
}
void EtaVEL::serialize(ostream &o){
// b.WriteVersion(EtaVEL::IsA());
char del = '|';
o << min << del;
o << max << del;
o << ds << del;
o << nPixels << del;
o << it << del;
o << totCont << del;
for(int i = 0; i < (nPixels+1)*(nPixels+1)+1; i++){
o << xPPos[i] << del;
o << yPPos[i] << del;
}
for(int i = 0; i < nPixels*nPixels+1; i++){
o << binCont[i] << del;
}
for(int i = 0; i < it; i++){
o << log[i].itN << del;
for(int j = 0; j < (nPixels+1)*(nPixels+1)+1; j++){
o << log[i].xPos[j] << del;
o << log[i].yPos[j] << del;
}
for(int j = 0; j < nPixels*nPixels+1; j++){
o << log[i].binCont[j] << del;
}
}
}
void EtaVEL::deserialize(istream &is){
delete[] xPPos;
delete[] yPPos;
delete[] binCont;
char del;
is >> min >> del;
is >> max >> del;
is >> ds >> del;
is >> nPixels >> del;
is >> it >> del;
is >> totCont >> del;
xPPos = new double[(nPixels+1)*(nPixels+1)+1];
yPPos = new double[(nPixels+1)*(nPixels+1)+1];
binCont = new double[nPixels*nPixels+1];
cout << "d";
for(int i = 0; i < (nPixels+1)*(nPixels+1)+1; i++){
is >> xPPos[i] >> del;
is >> yPPos[i] >> del;
}
cout << "d";
for(int i = 0; i < nPixels*nPixels+1; i++){
is >> binCont[i] >> del;
}
cout << "d";
for(int i = 0; i < it; i++){
is >> log[i].itN >> del;
log[i].xPos = new double[(nPixels+1)*(nPixels+1)+1];
log[i].yPos = new double[(nPixels+1)*(nPixels+1)+1];
log[i].binCont = new double[nPixels*nPixels+1];
for(int j = 0; j < (nPixels+1)*(nPixels+1)+1; j++){
is >> log[i].xPos[j] >> del;
is >> log[i].yPos[j] >> del;
}
for(int j = 0; j < nPixels*nPixels+1; j++){
is >> log[i].binCont[j] >> del;
}
cout << "d";
}
cout << endl;
}
void EtaVEL::Streamer(TBuffer &b){
if (b.IsReading()) {
Version_t v = b.ReadVersion();
delete[] xPPos;
delete[] yPPos;
delete[] binCont;
b >> min;
b >> max;
b >> ds;
b >> nPixels;
b >> it;
b >> totCont;
xPPos = new double[(nPixels+1)*(nPixels+1)+1];
yPPos = new double[(nPixels+1)*(nPixels+1)+1];
binCont = new double[nPixels*nPixels+1];
for(int i = 0; i < (nPixels+1)*(nPixels+1)+1; i++){
b >> xPPos[i];
b >> yPPos[i];
}
for(int i = 0; i < nPixels*nPixels+1; i++){
b >> binCont[i];
}
for(int i = 0; i < it; i++){
b >> log[i].itN;
log[i].xPos = new double[(nPixels+1)*(nPixels+1)+1];
log[i].yPos = new double[(nPixels+1)*(nPixels+1)+1];
log[i].binCont = new double[nPixels*nPixels+1];
for(int j = 0; j < (nPixels+1)*(nPixels+1)+1; j++){
b >> log[i].xPos[j];
b >> log[i].yPos[j];
}
for(int j = 0; j < nPixels*nPixels+1; j++){
b >> log[i].binCont[j];
}
}
} else {
b.WriteVersion(EtaVEL::IsA());
b << min;
b << max;
b << ds;
b << nPixels;
b << it;
b << totCont;
for(int i = 0; i < (nPixels+1)*(nPixels+1)+1; i++){
b << xPPos[i];
b << yPPos[i];
}
for(int i = 0; i < nPixels*nPixels+1; i++){
b << binCont[i];
}
for(int i = 0; i < it; i++){
b << log[i].itN;
for(int j = 0; j < (nPixels+1)*(nPixels+1)+1; j++){
b << log[i].xPos[j];
b << log[i].yPos[j];
}
for(int j = 0; j < nPixels*nPixels+1; j++){
b << log[i].binCont[j];
}
}
}
}

164
slsDetectorCalibration/interpolations/etaVEL/EtaVEL.h
View File

@ -1,164 +0,0 @@
#include <iostream>
#include <TGraph.h>
#include <TAxis.h>
#include <TMultiGraph.h>
#include <TH2D.h>
#include <TMath.h>
#include <TObject.h>
#include <TBuffer.h>
#include <TMatrixD.h>
#include <TDecompSVD.h>
//#include <TDecompQRH.h>
#include <TH1.h>
#include <TMath.h>
#include <vector>
#include <ostream>
#include <istream>
using namespace std;
#ifndef ETAVPS
#define ETAVPS
typedef struct {
int itN;
double *xPos;
double *yPos;
double *binCont;
} itLog;
class EtaVEL : public TObject{
public:
EtaVEL(int numberOfPixels = 25, double minn=0., double maxx=1., int nnx=160, int nny=160) : nPixels(numberOfPixels), min(minn), max(maxx), converged(0), nx(nnx), ny(nny), chi_sq(0){
//acc = 0.02;
ds = 0.005;
init();
}
void init(){
double pOffset = (max-min)/(double)nPixels;
xPPos = new double[(nPixels+1)*(nPixels+1)+1];
yPPos = new double[(nPixels+1)*(nPixels+1)+1];
binCont = new double[nPixels*nPixels+1];
totCont = 0.;
edgeL = new double[2*nPixels*(nPixels+1)+1];
for(int ii = 0; ii < 2*nPixels*(nPixels+1)+1; ii++){
edgeL[ii] = 1.0;
//cout << "ii " << ii << endl;
}
for(int x = 0; x < nPixels+1; x++){
for(int y = 0; y < nPixels+1; y++){
xPPos[getCorner(x,y)] = min + (double)x * pOffset;
yPPos[getCorner(x,y)] = min + (double)y * pOffset;
if(x < nPixels && y < nPixels) binCont[getBin(x,y)] = 0;
}
}
// edgeL[1] = 3.0;
updatePixelCorner();
it = 0;
log = new itLog[nIterations];
}
void fill(double x, double y, double amount = 1.){
totCont+=amount;
int bin = findBin(x,y);
if(bin < 0) {
//cout << "can not find bin x: " << x << " y: " << y << endl;
totCont-=amount;
}
binCont[bin]+=amount;
}
int getBin(int x, int y){
if(x < 0 || x >= nPixels || y < 0 || y >= nPixels){
//cout << "getBin: out of bounds : x " << x << " y " << y << endl;
return 0;
}
return y*nPixels+x+1;
}
int getXBin(int bin){
return (bin-1)%nPixels;
}
int getYBin(int bin){
return (bin-1)/nPixels;
}
int getCorner(int x, int y){
return y*(nPixels+1)+x+1;
}
int getEdgeX(int x,int row){
int ret = row*nPixels+x+1;
//cout << "| edge X x " << x << " row " << row << ": "<< ret << " | ";
return ret;
}
int getEdgeY(int col, int y){
int ret = nPixels*(nPixels+1)+col*nPixels+y+1;
//cout << "| edge Y col " << col << " y " << y << ": "<< ret << " | ";
return ret;
}
int getIt(){ return it; };
int getNPixels(){ return nPixels; }
double *getXPPos(){ return xPPos; }
double *getYPPos(){ return yPPos; }
void updatePixelCorner();
double *getPixelCorners(int x, int y);
int findBin(double xx, double yy);
void createLogEntry();
void updatePixelPos();
double *getSizeMap();
double *getChangeMap();
TH2D *getContent(int it=-1, int changeType = 0);
TMultiGraph *plotPixelBorder(int plotCenters=0);
TMultiGraph *plotLog(int stepSize=1, int maxIt=-1);
void printGrid();
TH1D *getCounts();
void serialize(ostream &o);
void deserialize(istream &is);
int converged ;
double getChiSq(){return chi_sq;};
private:
itLog *log;
int it;
const static int nIterations =10000;
int nx, ny;
int nPixels;
double *xPPos;
double *yPPos;
double *binCont;
double totCont;
double *edgeL;
// double acc;
double ds;
double min,max;
double chi_sq;
ClassDefNV(EtaVEL,1);
#pragma link C++ class EtaVEL-;
};
#endif

393
slsDetectorCalibration/interpolations/etaVEL/EtaVELTr.py
View File

@ -1,393 +0,0 @@
import numpy as np
import math
maxf = 2
minf = 0.5
class EtaVELTr:
def __init__(self, numberOfPixels = 25, minn=0., maxx=1.):
self.nPixels = numberOfPixels
self.nCorners = self.nPixels + 1
self.minEta = minn
self.maxEta = maxx
#self.corners = []
self.edgesX = []
self.edgesY = []
self.edgesE = []
self.edgesF = []
self.edgesG = []
self.edgesH = []
self.counts = []
self.pOffset = (self.maxEta-self.minEta)/self.nPixels
self.cPosX = []
self.cPosY = []
self.zPosX = []
self.zPosY = []
self.sqSums = []
self.cIteration = 0
self.bSqSum = 0
self.initGrid()
#self.calculatePixelCorners()
self.update()
def initGrid(self):
dd = 1 / math.sqrt(2)
#self.corners = [ [self.minEta + x * pOffset, self.minEta + y * pOffset] for y in range(self.nPixels)] for x in range(self.nPixels)
self.cPosX = [ [self.minEta + x * self.pOffset for x in range(self.nCorners)] for y in range(self.nCorners) ]
self.cPosY = [ [self.minEta + y * self.pOffset for x in range(self.nCorners)] for y in range(self.nCorners) ]
self.counts = [ [ [0,0,0,0] for x in range(self.nPixels) ] for y in range(self.nPixels) ]
self.edgesX = [ [ 1 for x in range(self.nCorners) ] for y in range(self.nCorners + 1) ]
self.edgesY = [ [ 1 for x in range(self.nCorners+1) ] for y in range(self.nCorners) ]
self.edgesE = [ [ dd for x in range(self.nPixels) ] for y in range(self.nPixels) ]
self.edgesF = [ [ dd for x in range(self.nPixels) ] for y in range(self.nPixels) ]
self.edgesG = [ [ dd for x in range(self.nPixels) ] for y in range(self.nPixels) ]
self.edgesH = [ [ dd for x in range(self.nPixels) ] for y in range(self.nPixels) ]
def update(self):
self.normalizeEdgeLengths()
self.calculateEdgeLengths2()
self.calculatePixelCorners2()
conv = False
out = 0
outList = []
tot = self.nPixels*self.nPixels*4
sqSum = 0
avg = self.getAvgCounts()
avgPS = self.getAvgCounts() + 1*math.sqrt(self.getAvgCounts())
avgMS = self.getAvgCounts() - 1*math.sqrt(self.getAvgCounts())
for y in range(self.nPixels):
for x in range(self.nPixels):
for t in range(4):
sqSum += (avg -self.counts[y][x][t]) * (avg -self.counts[y][x][t])
if self.counts[y][x][t] > avgPS or self.counts[y][x][t] < avgMS:
out += 1
outList.append([y,x,t,self.counts[y][x][t]])
outList = sorted(outList,key=lambda t: abs(self.counts[t[0]][t[1]][t[2]]/self.getAvgCounts()))
self.counts = [ [ [0,0,0,0] for x in range(self.nPixels) ] for y in range(self.nPixels) ]
print("There are {} of {} triangles out of 1 std ({} %)".format(out,tot,out/tot*100))
print("Total Sq Err: {}".format(sqSum))
self.sqSums.append(sqSum)
if len(self.sqSums) > 2 and np.diff(self.sqSums)[-1] > 0:
print("converged after {} steps: sqSums {} diff(sqSums) {}".format(self.cIteration,self.sqSums,np.diff(self.sqSums)))
conv = True
self.bSqSum = self.sqSums[-2]
self.cIteration += 1
return [conv,outList,sqSum]
def normalizeEdgeLengths(self):
sumL = 0
sumL += sum(map(sum,zip(*self.edgesX)))
sumL += sum(map(sum,zip(*self.edgesY)))
sumL += sum(map(sum,zip(*self.edgesE)))
sumL += sum(map(sum,zip(*self.edgesF)))
sumL += sum(map(sum,zip(*self.edgesG)))
sumL += sum(map(sum,zip(*self.edgesH)))
avgL = sumL/(4*self.nPixels*self.nPixels+2*self.nCorners*(self.nCorners+1))
print("total Sum is {} avg: {}".format(sumL,avgL))
self.edgesX = [ [ x/avgL for x in y ] for y in self.edgesX ]
self.edgesY = [ [ x/avgL for x in y ] for y in self.edgesY ]
self.edgesE = [ [ x/avgL for x in y ] for y in self.edgesE ]
self.edgesF = [ [ x/avgL for x in y ] for y in self.edgesF ]
self.edgesG = [ [ x/avgL for x in y ] for y in self.edgesG ]
self.edgesH = [ [ x/avgL for x in y ] for y in self.edgesH ]
def _shapeF(self,f):
f = (f - 1) * 0.6 + 1
if f > maxf:
return maxf
if f < minf:
return minf
return f
def calculateEdgeLengths2(self):
if self.getTotalCounts() == 0:
return
avg = self.getAvgCounts()
for y in range(self.nPixels):
for x in range(self.nPixels):
for t in range(4):
pc = self.counts[y][x][t]
if pc == 0:
f = maxf
else:
f = math.sqrt(avg/pc)
if pc > avg-math.sqrt(avg) and pc < avg+math.sqrt(avg):
f = 1.
sf = self._shapeF(f)
if t == 0:
self.edgesX[y][x] = self.edgesX[y][x] / sf
self.edgesE[y][x] = self.edgesE[y][x] / sf
self.edgesF[y][x] = self.edgesF[y][x] / sf
if t == 1:
self.edgesY[y][x+1] = self.edgesY[y][x+1] / sf
self.edgesF[y][x] = self.edgesF[y][x] / sf
self.edgesH[y][x] = self.edgesH[y][x] / sf
if t == 2:
self.edgesX[y+1][x] = self.edgesX[y+1][x] / sf
self.edgesH[y][x] = self.edgesH[y][x] / sf
self.edgesG[y][x] = self.edgesG[y][x] / sf
if t == 3:
self.edgesY[y][x] = self.edgesY[y][x] / sf
self.edgesG[y][x] = self.edgesG[y][x] / sf
self.edgesE[y][x] = self.edgesE[y][x] / sf
def calculatePixelCorners2(self):
w = 20
posMat = []
CrVx = np.zeros((self.nCorners,self.nCorners))
CrVy = np.zeros((self.nCorners,self.nCorners))
ZrVx = np.zeros((self.nPixels,self.nPixels))
ZrVy = np.zeros((self.nPixels,self.nPixels))
#boundary conditions matrix/vectors
BCposMatX = []
BCposMatY = []
BCrVx = []
BCrVy = []
for y in range(self.nCorners):
for x in range(self.nCorners):
BClineX = np.zeros((self.nCorners,self.nCorners))
BClineY = np.zeros((self.nCorners,self.nCorners))
if (x == 0 and y == 0) or \
(x == 0 and y == self.nPixels) or \
(x == self.nPixels and y == 0) or \
(x == self.nPixels and y == self.nPixels):
BClineX[y][x] = w
BClineY[y][x] = w
BCrVx.append(self.getCornerPos(y,x)[0] * w)
BCrVy.append(self.getCornerPos(y,x)[1] * w)
#print("bclinex shape {} zeros shape {}".format( BClineX.reshape((self.nCorners*self.nCorners,)).shape , np.zeros((self.nPixels*self.nPixels)).shape ))
BCposMatX.append(np.hstack((BClineX.reshape((self.nCorners*self.nCorners,)),np.zeros((self.nPixels*self.nPixels,)) )) )
BCposMatY.append(np.hstack((BClineY.reshape((self.nCorners*self.nCorners,)),np.zeros((self.nPixels*self.nPixels,)) )) )
elif x == 0 or x == self.nPixels:
BClineX[y][x] = w
#BClineY[y][x] = 1
BCrVx.append(self.getCornerPos(y,x)[0] * w)
#BCrVy.append(self.getCornerPos(y,x)[1])
BCposMatX.append(np.hstack((BClineX.reshape((self.nCorners*self.nCorners,)),np.zeros((self.nPixels*self.nPixels,)) )) )
#BCposMatY.append(np.hstack((BClineY.reshape((self.nCorners*self.nCorners,)),np.zeros((self.nPixels*self.nPixels,)) )) )
elif y == 0 or y == self.nPixels:
#BClineX[y][x] = 1
BClineY[y][x] = w
#BCrVx.append(self.getCornerPos(y,x)[0])
BCrVy.append(self.getCornerPos(y,x)[1] * w)
#BCposMatX.append(np.hstack((BClineX.reshape((self.nCorners*self.nCorners,)),np.zeros((self.nPixels*self.nPixels,)) )) )
BCposMatY.append(np.hstack((BClineY.reshape((self.nCorners*self.nCorners,)),np.zeros((self.nPixels*self.nPixels,)) )) )
eLength = 0
if x != 0:
eLength += self.edgesX[y][x-1]
if y != 0:
eLength += self.edgesY[y-1][x]
if x != self.nPixels:
eLength += self.edgesX[y][x]
if y != self.nPixels:
eLength += self.edgesY[y][x]
if y != 0 and x != 0:
eLength += self.edgesH[y-1][x-1]
if y != self.nPixels and x != 0:
eLength += self.edgesF[y][x-1]
if y != 0 and x != self.nPixels:
eLength += self.edgesG[y-1][x]
if y != self.nPixels and x != self.nPixels:
eLength += self.edgesE[y][x]
line = np.zeros((self.nCorners,self.nCorners))
lineZ = np.zeros((self.nPixels,self.nPixels))
if x != 0:
line[y][x-1] = - self.edgesX[y][x-1]/eLength
if y != 0:
line[y-1][x] = - self.edgesY[y-1][x]/eLength
if x != self.nPixels:
line[y][x+1] = - self.edgesX[y][x]/eLength
if y != self.nPixels:
line[y+1][x] = - self.edgesY[y][x]/eLength
if y != 0 and x != 0:
lineZ[y-1][x-1] = -self.edgesH[y-1][x-1]/eLength
if y != self.nPixels and x != 0:
lineZ[y][x-1] = -self.edgesF[y][x-1]/eLength
if y != 0 and x != self.nPixels:
lineZ[y-1][x] = -self.edgesG[y-1][x]/eLength
if y != self.nPixels and x != self.nPixels:
lineZ[y][x] = -self.edgesE[y][x]/eLength
line[y][x] = 1
CrVx[y][x] = 0
CrVy[y][x] = 0
posMat.append( \
np.hstack(( \
line.reshape((self.nCorners*self.nCorners,)), \
lineZ.reshape((self.nPixels*self.nPixels,)) \
)) \
)
for y in range(self.nPixels):
for x in range(self.nPixels):
line = np.zeros((self.nCorners,self.nCorners))
lineZ = np.zeros((self.nPixels,self.nPixels))
eLength = self.edgesE[y][x] + self.edgesF[y][x] +self.edgesG[y][x] +self.edgesH[y][x]
line[y][x] = -self.edgesE[y][x] / eLength
line[y][x+1] = -self.edgesF[y][x] / eLength
line[y+1][x] = -self.edgesG[y][x] / eLength
line[y+1][x+1] = -self.edgesH[y][x] / eLength
lineZ[y][x] = 1
ZrVx[y][x] = 0
ZrVy[y][x] = 0
posMat.append( \
np.hstack(( \
line.reshape((self.nCorners*self.nCorners,)), \
lineZ.reshape((self.nPixels*self.nPixels,)) \
)) \
)
CrVxFlat = CrVx.reshape((self.nCorners*self.nCorners,))
CrVyFlat = CrVy.reshape((self.nCorners*self.nCorners,))
ZrVxFlat = ZrVx.reshape((self.nPixels*self.nPixels,))
ZrVyFlat = ZrVy.reshape((self.nPixels*self.nPixels,))
posMat = np.asarray(posMat)
BCrVyFlat = np.asarray(BCrVy)
BCrVxFlat = np.asarray(BCrVx)
BCposMatX = np.asarray(BCposMatX)
BCposMatY = np.asarray(BCposMatY)
print ("BCposMatY vy {} shape posMat {}".format(BCposMatY.shape,posMat.shape))
FinalrVy = np.hstack((CrVyFlat,ZrVyFlat,BCrVyFlat))
FinalrVx = np.hstack((CrVxFlat,ZrVxFlat,BCrVxFlat))
FinalposMatX = np.vstack((posMat,BCposMatX))
FinalposMatY = np.vstack((posMat,BCposMatY))
print("posMat shape {}".format(posMat.shape))
print("posMatX shape {}".format(FinalposMatX.shape))
print("rVxFlat shape {}".format(FinalrVx.shape))
#print("posMat {}".format(FinalposMat))
#print("rVxFlat {}".format(FinalrVx))
xPos = np.linalg.lstsq(FinalposMatX,FinalrVx)[0]
yPos = np.linalg.lstsq(FinalposMatY,FinalrVy)[0]
print("xPosShape {} cutXPosShape {}".format(xPos.shape,xPos[:self.nCorners][:self.nCorners].shape))
self.cPosX = xPos[:self.nCorners*self.nCorners].reshape((self.nCorners,self.nCorners))
self.cPosY = yPos[:self.nCorners*self.nCorners].reshape((self.nCorners,self.nCorners))
self.zPosX = xPos[self.nCorners*self.nCorners:].reshape((self.nPixels,self.nPixels))
self.zPosY = yPos[self.nCorners*self.nCorners:].reshape((self.nPixels,self.nPixels))
def fill(self,yy,xx,count = 1):
[y,x,t] = self.getPixel(yy,xx)
self.counts[y][x][t] += count
def getCountDist(self):
c = []
for y in range(self.nPixels):
for x in range(self.nPixels):
c.append(self.counts[y][x])
return c
def getPixel(self,yy,xx, debug = False):
for y in range(self.nPixels):
for x in range(self.nPixels):
for t in range(4):
[v1x,v1y,v2x,v2y,v3x,v3y] = self.getTriangleCorner(y,x,t)
if self.pointInTriangle([xx,yy],[v1x,v1y],[v2x,v2y],[v3x,v3y]):
return [y,x,t]
if not debug:
raise Exception("not inside a pixel")
else:
print("no pixel found")
return [0,0]
#http://stackoverflow.com/questions/2049582/how-to-determine-a-point-in-a-2d-triangle
def trSign (self, p1, p2, p3):
return (p1[0] - p3[0]) * (p2[1] - p3[1]) - (p2[0] - p3[0]) * (p1[1] - p3[1]);
def pointInTriangle (self,pt, v1, v2, v3):
b1 = self.trSign(pt, v1, v2) < 0.0
b2 = self.trSign(pt, v2, v3) < 0.0
b3 = self.trSign(pt, v3, v1) < 0.0
return ((b1 == b2) and (b2 == b3));
def getAvgCounts(self):
return self.getTotalCounts() / self.nPixels/self.nPixels/4.
def getTotalCounts(self):
tot = 0
for y in range(self.nPixels):
for x in range(self.nPixels):
for t in range(4):
tot += self.counts[y][x][t]
return tot
#tl tr bl br
def getPixelCorners(self,iy,ix):
return self.getCornerPos(iy,ix) + self.getCornerPos(iy,ix+1) + self.getCornerPos(iy+1,ix) + self.getCornerPos(iy+1,ix+1)
def getTriangleCorner(self,iy,ix,tr):
if tr == 0:
return self.getCornerPos(iy,ix) + self.getCornerPos(iy,ix+1) + [self.zPosX[iy][ix],self.zPosY[iy][ix]]
if tr == 1:
return self.getCornerPos(iy,ix+1) + self.getCornerPos(iy+1,ix+1) + [self.zPosX[iy][ix],self.zPosY[iy][ix]]
if tr == 2:
return self.getCornerPos(iy+1,ix+1) + self.getCornerPos(iy+1,ix) + [self.zPosX[iy][ix],self.zPosY[iy][ix]]
if tr == 3:
return self.getCornerPos(iy+1,ix) + self.getCornerPos(iy,ix) + [self.zPosX[iy][ix],self.zPosY[iy][ix]]
def getCornerPos(self,iy,ix):
return [self.cPosX[iy][ix],self.cPosY[iy][ix]]
def getXEdgePos(self,iy,ix):
p1 = self.getCornerPos(iy,ix)
p2 = self.getCornerPos(iy,ix+1)
return [p1[0],p1[1],p2[0],p2[1]]
def getYEdgePos(self,iy,ix):
p1 = self.getCornerPos(iy,ix)
p2 = self.getCornerPos(iy+1,ix)
return [p1[0],p1[1],p2[0],p2[1]]
def getEEdgePos(self,iy,ix):
p1 = self.getCornerPos(iy,ix)
return [p1[0],p1[1],self.zPosX[iy][ix],self.zPosY[iy][ix]]
def getFEdgePos(self,iy,ix):
p1 = self.getCornerPos(iy,ix+1)
return [p1[0],p1[1],self.zPosX[iy][ix],self.zPosY[iy][ix]]
def getGEdgePos(self,iy,ix):
p1 = self.getCornerPos(iy+1,ix)
return [p1[0],p1[1],self.zPosX[iy][ix],self.zPosY[iy][ix]]
def getHEdgePos(self,iy,ix):
p1 = self.getCornerPos(iy+1,ix+1)
return [p1[0],p1[1],self.zPosX[iy][ix],self.zPosY[iy][ix]]

134
slsDetectorCalibration/interpolations/etaVEL/etaVELInterpolation.cpp
View File

@ -1,134 +0,0 @@
#include "interpolation_EtaVEL.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TROOT.h"
//#include "EtaVEL.h"
#include "EtaVEL.cpp"
/*
Zum erstellen der correction map ist createGainAndEtaFile(...) in EVELAlg.C der entry point.
Zum erstellen des HR images ist createImage(...) der entry point.
*/
interpolation_EtaVEL::interpolation_EtaVEL(int nx, int ny, int ns, double etamin, double etamax, int p) : slsInterpolation(nx, ny, ns), newEta(NULL), heta(NULL), plot(p) {
newEta = new EtaVEL(nSubPixels,etamin,etamax,nPixelsX, nPixelsY);
heta= new TH2F("heta","heta",50*nSubPixels, etamin,etamax,50*nSubPixels, etamin,etamax);
heta->SetStats(kFALSE);
}
interpolation_EtaVEL::~interpolation_EtaVEL() {
delete newEta;
delete heta;
}
void interpolation_EtaVEL::prepareInterpolation(int &ok, int maxit) {
int nit=0;
while ((newEta->converged != 1) && nit++<maxit) {
cout << " -------------- new step "<< nit << endl;
iterate();
}
if (plot) {
Draw();
gPad->Modified();
gPad->Update();
}
if (newEta->converged==1) ok=1; else ok=0;
}
int interpolation_EtaVEL::addToFlatField(Double_t *cluster, Double_t &etax, Double_t &etay) {
Double_t sum, totquad, sDum[2][2];
int corner =calcEta(cluster, etax, etay, sum, totquad, sDum);
//check if it's OK...should redo it every time?
//or should we fill a finer histogram and afterwards re-fill the newEta?
addToFlatField(etax, etay);
return corner;
}
int interpolation_EtaVEL::addToFlatField(Double_t etax, Double_t etay) {
// newEta->fill(etaX,etaY);
heta->Fill(etax,etay);
return 0;
}
void interpolation_EtaVEL::iterate() {
cout << " -------------- newEta refilled"<< endl;
for (int ibx=0; ibx<heta->GetNbinsX(); ibx++) {
for (int iby=0; iby<heta->GetNbinsY(); iby++) {
newEta->fill(heta->GetXaxis()->GetBinCenter(ibx+1),heta->GetYaxis()->GetBinCenter(iby+1),heta->GetBinContent(ibx+1,iby+1));
}
}
newEta->updatePixelPos();
cout << " -------------- pixelPosition updated"<< endl;
}
void interpolation_EtaVEL::DrawH() {
heta->Draw("col");
(newEta->plotPixelBorder())->Draw();
}
void interpolation_EtaVEL::getInterpolatedPosition(Int_t x, Int_t y, Double_t *cluster, Double_t &int_x, Double_t &int_y) {
Double_t etax, etay, sum, totquad, sDum[2][2];
int corner =calcEta(cluster, etax, etay, sum, totquad, sDum);
int bin = newEta->findBin(etax,etay);
if (bin<=0) {
int_x=-1;
int_y=-1;
return;
}
double subX = ((double)(newEta->getXBin(bin))+.5)/((double)newEta->getNPixels());
double subY = ((double)(newEta->getYBin(bin))+.5)/((double)newEta->getNPixels());
double dX, dY;
switch (corner) {
case TOP_LEFT:
dX=-1.;
dY=+1.;
break;
case TOP_RIGHT:
dX=+1.;
dY=+1.;
break;
case BOTTOM_LEFT:
dX=-1.;
dY=-1.;
break;
case BOTTOM_RIGHT:
dX=+1.;
dY=-1.;
break;
default:
dX=0;
dY=0;
}
int_x=((double)x)+ subX+0.5*dX;
int_y=((double)y)+ subY+0.5*dY;
// cout << corner << " " << subX<< " " << subY << " " << dX << " " << dY << " " << int_x << " " << int_y << endl;
};
// void interpolation_EtaVEL::Streamer(TBuffer &b){newEta->Streamer(b);};
void interpolation_EtaVEL::getInterpolatedBin(Double_t *cluster, Int_t &int_x, Int_t &int_y) {
Double_t etax, etay, sum, totquad, sDum[2][2];
int corner =calcEta(cluster, etax, etay, sum, totquad, sDum);
int bin = newEta->findBin(etax,etay);
if (bin<0) {
int_x=-1;
int_y=-1;
return;
}
int_x=newEta->getXBin(bin);
int_y=newEta->getYBin(bin);
};

55
slsDetectorCalibration/interpolations/etaVEL/etaVELInterpolation.h
View File

@ -1,55 +0,0 @@
#ifndef INTERPOLATION_ETAVEL_H
#define INTERPOLATION_ETAVEL_H
#include <slsInterpolation.h>
#include "EtaVEL.h"
//#include "TH2F.h"
//#include "EtaVEL.cpp"
//class EtaVEL;
class etaVELInterpolation: public etaInterpolationBase {
public:
interpolation_EtaVEL(int nx=40, int ny=160, int ns=25, double etamin=-0.02, double etamax=1.02, int p=0);
~interpolation_EtaVEL();
//create eta distribution, eta rebinnining etc.
//returns flat field image
void prepareInterpolation(int &ok){prepareInterpolation(ok,10000);};
void prepareInterpolation(int &ok, int maxit);
//create interpolated image
//returns interpolated image
//return position inside the pixel for the given photon
void getInterpolatedPosition(Int_t x, Int_t y, Double_t *data, Double_t &int_x, Double_t &int_y);
void getInterpolatedBin(Double_t *cluster, Int_t &int_x, Int_t &int_y);
int addToFlatField(Double_t *cluster, Double_t &etax, Double_t &etay);
int addToFlatField(Double_t etax, Double_t etay);
int setPlot(int p=-1) {if (p>=0) plot=p; return plot;};
// int WriteH(){newEta->Write("newEta"); heta->Write("heta");};
EtaVEL *setEta(EtaVEL *ev){if (ev) {delete newEta; newEta=ev;} return newEta;};
// TH2F *setEta(TH2F *ev){if (ev) {delete heta; heta=ev;} return heta;};
void iterate();
// void DrawH();
double getChiSq(){return newEta->getChiSq();};
protected:
EtaVEL *newEta;
// TH2F *heta;
int plot;
// ClassDefNV(interpolation_EtaVEL,1);
// #pragma link C++ class interpolation_EtaVEL-;
};
#endif

234
slsDetectorCalibration/interpolations/linearInterpolation.h
View File

@ -1,234 +0,0 @@
#ifndef LINEAR_INTERPOLATION_H
#define LINEAR_INTERPOLATION_H
//#include <TObject.h>
//#include <TTree.h>
//#include <TH2F.h>
#include "slsInterpolation.h"
class linearInterpolation : public slsInterpolation{
public:
linearInterpolation(int nx=400, int ny=400, int ns=25) : slsInterpolation(nx,ny,ns) {};
linearInterpolation(linearInterpolation *orig) : slsInterpolation(orig) {};
virtual void prepareInterpolation(int &ok){ok=1;};
virtual linearInterpolation* Clone() {
return new linearInterpolation(this);
};
//////////////////////////////////////////////////////////////////////////////
//////////// /*It return position hit for the event in input */ //////////////
virtual void getInterpolatedPosition(int x, int y, double *data, double &int_x, double &int_y)
{
double sDum[2][2];
double tot, totquad;
double etax,etay;
int corner;
corner=calcQuad(data, tot, totquad, sDum);
if (nSubPixels>2) {
calcEta(totquad, sDum, etax, etay);
}
getInterpolatedPosition(x, y, etax,etay, corner, int_x, int_y);
return;
};
virtual void getInterpolatedPosition(int x, int y, int *data, double &int_x, double &int_y)
{
double sDum[2][2];
double tot, totquad;
double etax,etay;
int corner;
corner=calcQuad(data, tot, totquad, sDum);
if (nSubPixels>2)
calcEta(totquad, sDum, etax, etay);
getInterpolatedPosition(x,y,etax,etay,corner,int_x,int_y);
return;
};
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cl,double &int_x, double &int_y) {
double eta_x, eta_y;
if (nSubPixels>2) {
double cc[2][2];
double *cluster[3];
cluster[0]=cl;
cluster[1]=cl+3;
cluster[2]=cl+6;
int xoff, yoff;
switch (quad) {
case BOTTOM_LEFT:
xoff=0;
yoff=0;
break;
case BOTTOM_RIGHT:
xoff=1;
yoff=0;
break;
case TOP_LEFT:
xoff=0;
yoff=1;
break;
case TOP_RIGHT:
xoff=1;
yoff=1;
break;
default:
;
}
cc[0][0]=cluster[yoff][xoff];
cc[1][0]=cluster[yoff+1][xoff];
cc[0][1]=cluster[yoff][xoff+1];
cc[1][1]=cluster[yoff+1][xoff+1];
calcEta(totquad,cc,eta_x,eta_y);
}
// cout << x << " " << y << " " << eta_x << " " << eta_y << " " << int_x << " " << int_y << endl;
return getInterpolatedPosition(x,y,eta_x, eta_y,quad,int_x,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];
int *cluster[3];
int xoff, yoff;
cluster[0]=cl;
cluster[1]=cl+3;
cluster[2]=cl+6;
switch (quad) {
case BOTTOM_LEFT:
xoff=0;
yoff=0;
break;
case BOTTOM_RIGHT:
xoff=1;
yoff=0;
break;
case TOP_LEFT:
xoff=0;
yoff=1;
break;
case TOP_RIGHT:
xoff=1;
yoff=1;
break;
default:
;
}
double etax, etay;
if (nSubPixels>2) {
cc[0][0]=cluster[yoff][xoff];
cc[1][0]=cluster[yoff+1][xoff];
cc[0][1]=cluster[yoff][xoff+1];
cc[1][1]=cluster[yoff+1][xoff+1];
calcEta(totquad,cc,etax,etay);
}
// cout << x << " " << y << " " << etax << " " << etay << " " << int_x << " " << int_y << endl;
return getInterpolatedPosition(x,y,etax, etay,quad,int_x,int_y);
}
virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int corner, double &int_x, double &int_y)
{
double xpos_eta,ypos_eta;
double dX,dY;
switch (corner)
{
case TOP_LEFT:
dX=-1.;
dY=0;
break;
case TOP_RIGHT:
dX=0;
dY=0;
break;
case BOTTOM_LEFT:
dX=-1.;
dY=-1.;
break;
case BOTTOM_RIGHT:
dX=0;
dY=-1.;
break;
default:
cout << "bad quadrant" << endl;
dX=0.;
dY=0.;
}
if (nSubPixels>2) {
xpos_eta=(etax)+dX;
ypos_eta=(etay)+dY;
} else {
xpos_eta=0.5*dX+0.25;
ypos_eta=0.5*dY+0.25;
}
int_x=((double)x) + xpos_eta;
int_y=((double)y) + ypos_eta;
// cout <<"**"<< x << " " << y << " " << xpos_eta << " " << ypos_eta << " " << corner << endl;
return;
};
////////////////////////////////////////////////////////////////////////////////////////////////////////
virtual void getPositionETA3(int x, int y, double *data, double &int_x, double &int_y)
{
double sDum[2][2];
double tot, totquad;
double eta3x,eta3y;
calcQuad(data, tot, totquad, sDum);
calcEta3(data,eta3x, eta3y,tot);
double xpos_eta,ypos_eta;
xpos_eta=eta3x;
ypos_eta=eta3y;
int_x=((double)x) + xpos_eta;
int_y=((double)y) + ypos_eta;
return;
};
////////////////////////////////////////////////////////////////////////////////////////////////////////
virtual int addToFlatField(double *cluster, double &etax, double &etay){};
virtual int addToFlatField(int *cluster, double &etax, double &etay){};
virtual int addToFlatField(double etax, double etay){};
virtual int addToFlatField(double totquad,int quad,double *cl,double &etax, double &etay) {};
virtual int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay) {};
protected:
;
};
#endif

104
slsDetectorCalibration/interpolations/noInterpolation.h
View File

@ -1,104 +0,0 @@
#ifndef NO_INTERPOLATION_H
#define NO_INTERPOLATION_H
/* #ifdef MYROOT1 */
/* #include <TObject.h> */
/* #include <TTree.h> */
/* #include <TH2F.h> */
/* #include <TROOT.h> */
/* #include <TRandom.h> */
/* #endif */
#include <cstdlib>
#include "slsInterpolation.h"
class noInterpolation : public slsInterpolation{
public:
noInterpolation(int nx=400, int ny=400, int ns=25) : slsInterpolation(nx,ny,ns) {};// {eventGenerator=new TRandom();};
noInterpolation(noInterpolation *orig) : slsInterpolation(orig){};
virtual void prepareInterpolation(int &ok){ok=1;};
//////////////////////////////////////////////////////////////////////////////
//////////// /*It return position hit for the event in input */ //////////////
virtual noInterpolation* Clone() {
return new noInterpolation(this);
};
virtual void getInterpolatedPosition(int x, int y, double *data, double &int_x, double &int_y)
{
//Random coordinate in the Pixel reference
int_x = x + ((double)rand())/((double)RAND_MAX) -0.5;//eventGenerator->Uniform(-0.5,0.5);
int_y = y + ((double)rand())/((double)RAND_MAX) -0.5;//eventGenerator->Uniform(-0.5,0.5);
return ;
};
virtual void getInterpolatedPosition(int x, int y, int *data, double &int_x, double &int_y)
{
return getInterpolatedPosition(x, y, (double*)NULL, int_x, int_y);
}
virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int corner, double &int_x, double &int_y)
{
getInterpolatedPosition(x, y, (double*)NULL, int_x, int_y);
};
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cl,double &etax, double &etay){
getInterpolatedPosition(x, y, (double*)NULL, etax, etay);
};
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,int *cl,double &etax, double &etay){
getInterpolatedPosition(x, y, (double*)NULL, etax, etay);
};
//////////////////////////////////////////////////////////////////////////////////////
virtual void getPositionETA3(int x, int y, double *data, double &int_x, double &int_y)
{
//Random coordinate in the Pixel reference
int_x = x + ((double)rand())/((double)RAND_MAX) -0.5;//eventGenerator->Uniform(-0.5,0.5);
int_y = y + ((double)rand())/((double)RAND_MAX) -0.5;//eventGenerator->Uniform(-0.5,0.5);
return ;
};
virtual void getPositionETA3(int x, int y, int *data, double &int_x, double &int_y)
{
return getPositionETA3(x, y, (double*)NULL, int_x, int_y);
};
//////////////////////////////////////////////////////////////////////////////////////
virtual int addToFlatField(double *cluster, double &etax, double &etay){return 0;};
virtual int addToFlatField(int *cluster, double &etax, double &etay){return 0;};
virtual int addToFlatField(double etax, double etay){return 0;};
virtual int addToFlatField(double totquad,int quad,double *cl,double &etax, double &etay){return 0;};
virtual int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay){return 0;};
virtual void resetFlatField(){};
protected:
;
// TRandom *eventGenerator;
// ClassDefNV(slsInterpolation,1);
// #pragma link C++ class slsInterpolation-;
};
#endif

503
slsDetectorCalibration/interpolations/slsInterpolation.h
View File

@ -1,503 +0,0 @@
#ifndef SLS_INTERPOLATION_H
#define SLS_INTERPOLATION_H
#include <cstdlib>
#ifndef MY_TIFF_IO_H
#include "tiffIO.h"
#endif
#ifndef DEF_QUAD
#define DEF_QUAD
enum quadrant {
TOP_LEFT=0,
TOP_RIGHT=1,
BOTTOM_LEFT=2,
BOTTOM_RIGHT=3,
UNDEFINED_QUADRANT=-1
};
#endif
#include <memory.h>
#include <stdio.h>
#include <iostream>
using namespace std;
//#ifdef MYROOT1
//: public TObject
//#endif
class slsInterpolation
{
public:
slsInterpolation(int nx=400, int ny=400, int ns=25) :nPixelsX(nx), nPixelsY(ny), nSubPixels(ns), id(0) {
hint=new int[ns*nx*ns*ny];
};
slsInterpolation(slsInterpolation *orig){
nPixelsX=orig->nPixelsX;
nPixelsY=orig->nPixelsY;
nSubPixels=orig->nSubPixels;
hint=new int[nSubPixels*nPixelsX*nSubPixels*nPixelsY];
memcpy(hint, orig->hint,nSubPixels*nPixelsX*nSubPixels*nPixelsY*sizeof(int));
};
virtual int setId(int i) {id=i; return id;};
virtual slsInterpolation* Clone() =0; /*{
return new slsInterpolation(this);
}*/
int getNSubPixels() {return nSubPixels;};
int setNSubPixels(int ns) {
if (ns>0 && ns!=nSubPixels) {
delete [] hint;
nSubPixels=ns;
hint=new int[nSubPixels*nPixelsX*nSubPixels*nPixelsY];
}
return nSubPixels;
}
int getImageSize(int &nnx, int &nny, int &ns) {
nnx=nSubPixels*nPixelsX;
nny=nSubPixels*nPixelsY;
ns=nSubPixels;
return nSubPixels*nSubPixels*nPixelsX*nPixelsY;
};
//create eta distribution, eta rebinnining etc.
//returns flat field image
virtual void prepareInterpolation(int &ok)=0;
//create interpolated image
//returns interpolated image
virtual int *getInterpolatedImage(){
// cout << "return interpolated image " << endl;
/* for (int i=0; i<nSubPixels* nSubPixels* nPixelsX*nPixelsY; i++) { */
/* cout << i << " " << hint[i] << endl; */
/* } */
return hint;
};
void *writeInterpolatedImage(const char * imgname) {
//cout << "!" <<endl;
float *gm=NULL;
int *dummy=getInterpolatedImage();
gm=new float[ nSubPixels* nSubPixels* nPixelsX*nPixelsY];
if (gm) {
for (int ix=0; ix<nPixelsX*nSubPixels; ix++) {
for (int iy=0; iy<nPixelsY*nSubPixels; iy++) {
gm[iy*nPixelsX*nSubPixels+ix]=dummy[iy*nPixelsX*nSubPixels+ix];
}
}
WriteToTiff(gm, imgname,nSubPixels* nPixelsX ,nSubPixels* nPixelsY);
delete [] gm;
} else cout << "Could not allocate float image " << endl;
return NULL;
}
//return position inside the pixel for the given photon
virtual void getInterpolatedPosition(int x, int y, double *data, double &int_x, double &int_y)=0;
virtual void getInterpolatedPosition(int x, int y, int *data, double &int_x, double &int_y)=0;
//return position inside the pixel for the given photon
virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int quad, double &int_x, double &int_y)=0;
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,int *cluster,double &etax, double &etay)=0;
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cluster,double &etax, double &etay)=0;
//return position inside the pixel for the given photon
virtual void clearInterpolatedImage() {
for (int ix=0; ix<nPixelsX*nSubPixels; ix++) {
for (int iy=0; iy<nPixelsY*nSubPixels; iy++) {
hint[iy*nPixelsX*nSubPixels+ix]=0;
}
}
};
virtual int *addToImage(double int_x, double int_y){
int iy=((double)nSubPixels)*int_y;
int ix=((double)nSubPixels)*int_x;
if (ix>=0 && ix<(nPixelsX*nSubPixels) && iy<(nSubPixels*nPixelsY) && iy>=0 ){
// cout << int_x << " " << int_y << " " << " " << ix << " " << iy << " " << ix+iy*nPixelsX*nSubPixels << " " << hint[ix+iy*nPixelsX*nSubPixels];
(*(hint+ix+iy*nPixelsX*nSubPixels))+=1;
// cout << " " << hint[ix+iy*nPixelsX*nSubPixels] << endl;
}// else
// cout << "bad! "<< int_x << " " << int_y << " " << " " << ix << " " << iy << " " << ix+iy*nPixelsX*nSubPixels << endl;
return hint;
};
virtual int addToFlatField(double *cluster, double &etax, double &etay)=0;
virtual int addToFlatField(int *cluster, double &etax, double &etay)=0;
virtual int addToFlatField(double totquad,int quad,int *cl,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 *getFlatField(){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 *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 void resetFlatField()=0;
//virtual void Streamer(TBuffer &b);
static int calcQuad(int *cl, double &sum, double &totquad, double sDum[2][2]){
double cli[3*3];//=new int[3*3];
for (int i=0; i<9; i++)
cli[i]=cl[i];
return calcQuad(cli, sum, totquad, sDum);
}
static int calcQuad(double *cl, double &sum, double &totquad, double sDum[2][2]){
int corner = UNDEFINED_QUADRANT;
/* double *cluster[3]; */
/* cluster[0]=cl; */
/* cluster[1]=cl+3; */
/* cluster[2]=cl+6; */
sum=0;
double sumBL=0;
double sumTL=0;
double sumBR=0;
double sumTR=0;
int xoff=0, yoff=0;
for (int ix=0; ix<3; ix++) {
for (int iy=0; iy<3; iy++) {
sum+=cl[ix+3*iy];
if (ix<=1 && iy<=1) sumBL+=cl[ix+iy*3];
if (ix<=1 && iy>=1) sumTL+=cl[ix+iy*3];
if (ix>=1 && iy<=1) sumBR+=cl[ix+iy*3];
if (ix>=1 && iy>=1) sumTR+=cl[ix+iy*3];
}
}
/* sDum[0][0] = cluster[0][0]; sDum[1][0] = cluster[1][0]; */
/* sDum[0][1] = cluster[0][1]; sDum[1][1] = cluster[1][1]; */
corner = BOTTOM_LEFT;
totquad=sumBL;
xoff=0;
yoff=0;
if(sumTL >= totquad){
/* sDum[0][0] = cluster[1][0]; sDum[1][0] = cluster[2][0]; */
/* sDum[0][1] = cluster[1][1]; sDum[1][1] = cluster[2][1]; */
corner = TOP_LEFT;
totquad=sumTL;
xoff=0;
yoff=1;
}
if(sumBR >= totquad){
/* sDum[0][0] = cluster[0][1]; sDum[1][0] = cluster[1][1]; */
/* sDum[0][1] = cluster[0][2]; sDum[1][1] = cluster[1][2]; */
xoff=1;
yoff=0;
corner = BOTTOM_RIGHT;
totquad=sumBR;
}
if(sumTR >= totquad){
xoff=1;
yoff=1;
/* sDum[0][0] = cluster[1][1]; sDum[1][0] = cluster[2][1]; */
/* sDum[0][1] = cluster[1][2]; sDum[1][1] = cluster[2][2]; */
corner = TOP_RIGHT;
totquad=sumTR;
}
for (int ix=0; ix<2; ix++) {
for (int iy=0; iy<2; iy++) {
sDum[iy][ix] = cl[ix+xoff+(iy+yoff)*3];
}
}
return corner;
}
static int calcEta(double totquad, double sDum[2][2], double &etax, double &etay){
double t,r;
if (totquad>0) {
t = sDum[1][0] + sDum[1][1];
r = sDum[0][1] + sDum[1][1];
etax=r/totquad;
etay=t/totquad;
}
return 0;
}
static int calcEta(double *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]) {
int corner = calcQuad(cl,sum,totquad,sDum);
calcEta(totquad, sDum, etax, etay);
return corner;
}
static int calcEta(int *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]) {
int corner = calcQuad(cl,sum,totquad,sDum);
calcEta(totquad, sDum, etax, etay);
return corner;
}
static int calcEtaL(double totquad, int corner, double sDum[2][2], double &etax, double &etay){
double t,r, toth, totv;
if (totquad>0) {
switch(corner) {
case TOP_LEFT:
t = sDum[1][1];
r = sDum[0][1] ;
toth=sDum[0][1]+sDum[0][0];
totv=sDum[0][1]+sDum[1][1];
break;
case TOP_RIGHT:
t = sDum[1][0] ;
r = sDum[0][1] ;
toth=sDum[0][1]+sDum[0][0];
totv=sDum[1][0]+sDum[0][0];
break;
case BOTTOM_LEFT:
r = sDum[1][1] ;
t = sDum[1][1] ;
toth=sDum[1][0]+sDum[1][1];
totv=sDum[0][1]+sDum[1][1];
break;
case BOTTOM_RIGHT:
t = sDum[1][0] ;
r = sDum[1][1] ;
toth=sDum[1][0]+sDum[1][1];
totv=sDum[1][0]+sDum[0][0];
break;
default:
etax=-1000;
etay=-1000;
return 0;
}
//etax=r/totquad;
//etay=t/totquad;
etax=r/toth;
etay=t/totv;
}
return 0;
}
static int calcEtaL(double *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]) {
int corner = calcQuad(cl,sum,totquad,sDum);
calcEtaL(totquad, corner, sDum, etax, etay);
return corner;
}
static int calcEtaL(int *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]) {
int corner = calcQuad(cl,sum,totquad,sDum);
calcEtaL(totquad, corner, sDum, etax, etay);
return corner;
}
static int calcEtaC3(double *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]){
int corner = calcQuad(cl,sum,totquad,sDum);
calcEta(sum, sDum, etax, etay);
return corner;
}
static int calcEtaC3(int *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]){
int corner = calcQuad(cl,sum,totquad,sDum);
calcEta(sum, sDum, etax, etay);
return corner;
}
static int calcEta3(double *cl, double &etax, double &etay, double &sum) {
double l=0,r=0,t=0,b=0, val;
sum=0;
// int quad;
for (int ix=0; ix<3; ix++) {
for (int iy=0; iy<3; iy++) {
val=cl[ix+3*iy];
sum+=val;
if (iy==0) l+=val;
if (iy==2) r+=val;
if (ix==0) b+=val;
if (ix==2) t+=val;
}
}
if (sum>0) {
etax=(-l+r)/sum;
etay=(-b+t)/sum;
}
/* if (etax<-1 || etax>1 || etay<-1 || etay>1) { */
/* cout << "**********" << etax << " " << etay << endl; */
/* for (int ix=0; ix<3; ix++) { */
/* for (int iy=0; iy<3; iy++) { */
/* cout << cl[iy+3*ix] << "\t" ; */
/* } */
/* cout << endl; */
/* } */
/* cout << sum << " " << l << " " << r << " " << t << " " << b << endl; */
/* } */
if (etax>=0 && etay>=0)
return TOP_RIGHT;
if (etax<0 && etay>=0)
return TOP_LEFT;
if (etax<0 && etay<0)
return BOTTOM_LEFT;
return BOTTOM_RIGHT;
}
static int calcEta3(int *cl, double &etax, double &etay, double &sum) {
double cli[9];
for (int ix=0; ix<9; ix++) cli[ix]=cl[ix];
return calcEta3(cli, etax, etay, sum);
}
/* static int calcMyEta(double totquad, int quad, double *cl, double &etax, double &etay) { */
/* double l,r,t,b, sum; */
/* int yoff; */
/* switch (quad) { */
/* case BOTTOM_LEFT: */
/* case BOTTOM_RIGHT: */
/* yoff=0; */
/* break; */
/* case TOP_LEFT: */
/* case TOP_RIGHT: */
/* yoff=1; */
/* break; */
/* default: */
/* ; */
/* } */
/* l=cl[0+yoff*3]+cl[0+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]; */
/* t=cl[0+yoff*3+3]+cl[1+yoff*3+3]*cl[0+yoff*3+3]; */
/* sum=t+b; */
/* if (sum>0) { */
/* etax=(-l+r)/sum; */
/* etay=(+t)/sum; */
/* } */
/* return -1; */
/* } */
/* static int calcMyEta(double totquad, int quad, int *cl, double &etax, double &etay) { */
/* double l,r,t,b, sum; */
/* int yoff; */
/* switch (quad) { */
/* case BOTTOM_LEFT: */
/* case BOTTOM_RIGHT: */
/* yoff=0; */
/* break; */
/* case TOP_LEFT: */
/* case TOP_RIGHT: */
/* yoff=1; */
/* break; */
/* default: */
/* ; */
/* } */
/* l=cl[0+yoff*3]+cl[0+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]; */
/* t=cl[0+yoff*3+3]+cl[1+yoff*3+3]*cl[0+yoff*3+3]; */
/* sum=t+b; */
/* if (sum>0) { */
/* etax=(-l+r)/sum; */
/* etay=(+t)/sum; */
/* } */
/* return -1; */
/* } */
/* static int calcEta3X(double *cl, double &etax, double &etay, double &sum) { */
/* double l,r,t,b; */
/* sum=cl[0]+cl[1]+cl[2]+cl[3]+cl[4]+cl[5]+cl[6]+cl[7]+cl[8]; */
/* if (sum>0) { */
/* l=cl[3]; */
/* r=cl[5]; */
/* b=cl[1]; */
/* t=cl[7]; */
/* etax=(-l+r)/sum; */
/* etay=(-b+t)/sum; */
/* } */
/* return -1; */
/* } */
/* static int calcEta3X(int *cl, double &etax, double &etay, double &sum) { */
/* double l,r,t,b; */
/* sum=cl[0]+cl[1]+cl[2]+cl[3]+cl[4]+cl[5]+cl[6]+cl[7]+cl[8]; */
/* if (sum>0) { */
/* l=cl[3]; */
/* r=cl[5]; */
/* b=cl[1]; */
/* t=cl[7]; */
/* etax=(-l+r)/sum; */
/* etay=(-b+t)/sum; */
/* } */
/* return -1; */
/* } */
protected:
int nPixelsX, nPixelsY;
int nSubPixels;
int id;
int *hint;
};
#endif