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

some changes in slsCalibration for moench

Browse Source
This commit is contained in:
bergamaschi
2018-03-15 12:26:45 +01:00
parent e95b444908
commit 10209b75df
22 changed files with 842 additions and 1050 deletions
Download Patch File Download Diff File Expand all files Collapse all files

76
slsDetectorCalibration/interpolations/etaInterpolationAdaptiveBins.h
View File

@ -181,53 +181,73 @@ class etaInterpolationAdaptiveBins : public etaInterpolationPosXY {
double avg=tot_eta/((double)(nSubPixels*nSubPixels));
cout << "total eta entries is :"<< tot_eta << " avg: "<< avg << endl;
cout << "Start " << endl;
double old_diff=calcDiff(avg, hhx, hhy), new_diff=old_diff+1;
cout << " diff= " << new_diff << endl;
double old_diff=calcDiff(avg, hhx, hhy), new_diff=old_diff+1, best_diff=old_diff;
cout << " diff= " << old_diff << endl;
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
while (iint<nint) {
cout << "Iteration " << iint << endl;
iterate(newhhx,newhhy);
new_diff=calcDiff(avg, newhhx, newhhy);
cout << " diff= " << new_diff << endl;
#ifdef SAVE_ALL
for (int ii=0; ii<etabins*etabins; ii++) {
etah[ii]=newhhx[ii];
if (etah[ii]>1 || etah[ii]<0 ) cout << "***"<< ii << etah[ii] << endl;
/* #ifdef SAVE_ALL */
/* for (int ii=0; ii<etabins*etabins; ii++) { */
/* etah[ii]=newhhx[ii]; */
/* if (etah[ii]>1 || etah[ii]<0 ) cout << "***"<< ii << etah[ii] << endl; */
}
sprintf(tit,"/scratch/neweta_hhx_%d.tiff",iint);
WriteToTiff(etah, tit, etabins, etabins);
/* } */
/* sprintf(tit,"/scratch/neweta_hhx_%d.tiff",iint); */
/* WriteToTiff(etah, tit, etabins, etabins); */
for (int ii=0; ii<etabins*etabins; ii++) {
etah[ii]=newhhy[ii];
if (etah[ii]>1 || etah[ii]<0 ) cout << "***"<< ii << etah[ii] << endl;
/* for (int ii=0; ii<etabins*etabins; ii++) { */
/* etah[ii]=newhhy[ii]; */
/* if (etah[ii]>1 || etah[ii]<0 ) cout << "***"<< ii << etah[ii] << endl; */
/* } */
/* sprintf(tit,"/scratch/neweta_hhy_%d.tiff",iint); */
/* WriteToTiff(etah, tit, etabins, etabins); */
/* #endif */
if (new_diff<best_diff) {
best_diff=new_diff;
besthhx=newhhx;
besthhy=newhhy;
}
sprintf(tit,"/scratch/neweta_hhy_%d.tiff",iint);
WriteToTiff(etah, tit, etabins, etabins);
#endif
// if (new_diff<old_diff) {
delete [] hhx;
delete [] hhy;
hhx=newhhx;
hhy=newhhy;
newhhx=new float[nbeta*nbeta]; //profile x
newhhy=new float[nbeta*nbeta]; //profile y
old_diff=new_diff;
// } /* else { */
if (hhx!=besthhx)
delete [] hhx;
if (hhy!=besthhy)
delete [] hhy;
hhx=newhhx;
hhy=newhhy;
newhhx=new float[nbeta*nbeta]; //profile x
newhhy=new float[nbeta*nbeta]; //profile y
old_diff=new_diff;
//} /* else { */
/* cout << "Difference not decreasing after "<< iint << " iterations (" << old_diff << " < " << new_diff << ")"<< endl; */
/* break; */
/* } */
iint++;
iint++;
}
delete [] newhhx;
delete [] newhhy;
delete [] newhhx;
delete [] newhhy;
if (hhx!=besthhx)
delete [] hhx;
if (hhy!=besthhy)
delete [] hhy;
hhx=besthhx;
hhy=besthhy;

291
slsDetectorCalibration/interpolations/etaInterpolationBase.h
View File

@ -14,10 +14,16 @@
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) {
if (nb<=0)
nbeta=nSubPixels*10;
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) {
cout << "aaa:" <<endl;
etamin=-1;
etamax=2;
}
@ -34,6 +40,7 @@ class etaInterpolationBase : public slsInterpolation {
#endif
//cout << nbeta << " " << etamin << " " << etamax << endl;
};
etaInterpolationBase(etaInterpolationBase *orig): slsInterpolation(orig){
@ -60,11 +67,11 @@ class etaInterpolationBase : public slsInterpolation {
};
virtual etaInterpolationBase* Clone() {
virtual etaInterpolationBase* Clone()=0;/*{
return new etaInterpolationBase(this);
};
*/
@ -179,8 +186,6 @@ class etaInterpolationBase : public slsInterpolation {
#endif
virtual void prepareInterpolation(int &ok){};
/* ////////////////////////////////////////////////////////////////////////////// */
@ -214,277 +219,35 @@ float *gethhx()
//////////////////////////////////////////////////////////////////////////////
//////////// /*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);
/* virtual void getInterpolatedPosition(int x, int y, int *data, double &int_x, double &int_y)=0; */
/* virtual void getInterpolatedPosition(int x, int y, double *data, double &int_x, double &int_y)=0; */
/* virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cl,double &int_x, double &int_y)=0; */
/* virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,int *cl,double &int_x, double &int_y)=0; */
/* virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int corner, double &int_x, double &int_y)=0; */
return;
};
/* virtual int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay)=0; */
/* virtual int addToFlatField(double totquad,int quad,double *cl,double &etax, double &etay)=0; */
/* virtual int addToFlatField(double *cluster, double &etax, double &etay)=0; */
/* virtual int addToFlatField(int *cluster, double &etax, double &etay)=0; */
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cl,double &int_x, double &int_y) {
double cc[2][2];
double *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);
}
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) {
#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) {
ex=0;
cout << "x*"<< ex << endl;
}
if (ex>=nbeta) {
ex=nbeta-1;
cout << "x?"<< ex << endl;
}
if (ey<0) {
ey=0;
cout << "y*"<< ey << endl;
}
if (ey>=nbeta) {
ey=nbeta-1;
cout << "y?"<< ey << endl;
}
xpos_eta=(((double)hhx[(ey*nbeta+ex)]))+dX ;///((double)nSubPixels);
ypos_eta=(((double)hhy[(ey*nbeta+ex)]))+dY ;///((double)nSubPixels);
//else
//return 0;
#endif
} 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 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,double *cl,double &etax, double &etay) {
double cc[2][2];
double *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:
;
}
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];
/* 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(double etax, double etay){
int ex,ey;
#ifdef MYROOT1
heta->Fill(etax,etay);
#endif
#ifndef MYROOT1
int ex,ey;
ex=(etax-etamin)/etastep;
ey=(etay-etamin)/etastep;
// cout << etax << " " << ex << " " << etay << " " << ey << " " << ey*nbeta+ex << endl;
if (ey<nbeta && ex<nbeta && ex>=0 && ey>=0)
heta[ey*nbeta+ex]++;
// cout << "*"<< etax << " " << etay << endl;
/* cout << etax << " " << etay << " " << ex << " " << ey << " " << ey*nbeta+ex << endl; */
/* cout <<"********"<< endl << endl ; */
#endif
return 0;
};
};
// virtual void prepareInterpolation(int &ok)=0;
protected:
#ifdef MYROOT1

42
slsDetectorCalibration/interpolations/etaInterpolationPosXY.h
View File

@ -2,24 +2,27 @@
#define ETA_INTERPOLATION_POSXY_H
#include "tiffIO.h"
//#include "tiffIO.h"
#include "etaInterpolationBase.h"
#include "eta2InterpolationBase.h"
#include "eta3InterpolationBase.h"
class etaInterpolationPosXY : public etaInterpolationBase{
class etaInterpolationPosXY : public virtual etaInterpolationBase{
public:
etaInterpolationPosXY(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny,ns, nb, emin,emax){};
etaInterpolationPosXY(int nx=400, int ny=400, int ns=25, int 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){};
etaInterpolationPosXY(etaInterpolationPosXY *orig): etaInterpolationBase(orig) {};
virtual etaInterpolationPosXY* Clone() {
virtual etaInterpolationPosXY* Clone()=0;/** {
return new etaInterpolationPosXY(this);
};
};*/
virtual void prepareInterpolation(int &ok)
{
cout <<"?"<< endl;
ok=1;
#ifdef MYROOT1
if (hhx) delete hhx;
@ -114,10 +117,6 @@ class etaInterpolationPosXY : public etaInterpolationBase{
// cout << "y " << nbeta << " " << (ii+1)*tot_eta_x*bsize << " " << ii << endl;
}
#ifdef SAVE_ALL
@ -150,4 +149,25 @@ class etaInterpolationPosXY : public etaInterpolationBase{
};
class eta2InterpolationPosXY : public virtual eta2InterpolationBase, public virtual etaInterpolationPosXY {
public:
eta2InterpolationPosXY(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny,ns, nb, emin,emax),eta2InterpolationBase(nx,ny,ns, nb, emin,emax),etaInterpolationPosXY(nx,ny,ns, nb, emin,emax){
// 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

36
slsDetectorCalibration/interpolations/noInterpolation.h
View File

@ -38,13 +38,31 @@ class noInterpolation : public slsInterpolation{
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, NULL, int_x, 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, NULL, etax, 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)
@ -55,13 +73,27 @@ class noInterpolation : public slsInterpolation{
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;};
protected:
;
// TRandom *eventGenerator;

217
slsDetectorCalibration/interpolations/slsInterpolation.h
View File

@ -8,7 +8,10 @@
#endif
#include <cstdlib>
#ifndef MY_TIFF_IO_H
#include "tiffIO.h"
#endif
#ifndef DEF_QUAD
#define DEF_QUAD
enum quadrant {
@ -88,7 +91,7 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
void *writeInterpolatedImage(const char * imgname) {
cout << "!" <<endl;
//cout << "!" <<endl;
float *gm=NULL;
gm=new float[ nSubPixels* nSubPixels* nPixelsX*nPixelsY];
if (gm) {
@ -107,8 +110,10 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
//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;
@ -138,7 +143,9 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
#endif
#ifndef MYROOT1
virtual int *addToImage(double int_x, double int_y){ int iy=nSubPixels*int_y; int ix=nSubPixels*int_x;
virtual int *addToImage(double int_x, double int_y){
int iy=((double)nSubPixels)*int_y;
int ix=((double)nSubPixels)*int_x;
// cout << int_x << " " << int_y << " " << " " << ix << " " << iy << " " << ix+iy*nPixelsX*nSubPixels << endl;
if (ix>=0 && ix<(nPixelsX*nSubPixels) && iy<(nSubPixels*nPixelsY) && iy>=0 )(*(hint+ix+iy*nPixelsX*nSubPixels))+=1;
return hint;
@ -147,9 +154,10 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
virtual int addToFlatField(double *cluster, double &etax, double &etay)=0;
virtual int addToFlatField(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;
#ifdef MYROOT1
virtual TH2D *getFlatField(){return NULL;};
@ -167,6 +175,14 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
//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]){
@ -175,48 +191,64 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
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+=cluster[iy][ix];
if (ix<=1 && iy<=1) sumBL+=cluster[iy][ix];
if (ix<=1 && iy>=1) sumTL+=cluster[iy][ix];
if (ix>=1 && iy<=1) sumBR+=cluster[iy][ix];
if (ix>=1 && iy>=1) sumTR+=cluster[iy][ix];
}
}
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 sumBL = cluster[0][0] + cluster[1][0] + cluster[0][1] + cluster[1][1]; //2 ->BL
double sumTL = cluster[1][0] + cluster[2][0] + cluster[2][1] + cluster[1][1]; //0 ->TL
double sumBR = cluster[0][1] + cluster[0][2] + cluster[1][2] + cluster[1][1]; //3 ->BR
double sumTR = cluster[1][2] + cluster[2][1] + cluster[2][2] + cluster[1][1]; //1 ->TR
double sumMax = 0;
double t, r;
// 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];
/* 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;
sumMax=sumBL;
// }
if(sumTL >= 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];
totquad=sumBL;
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;
sumMax=sumTL;
totquad=sumTL;
xoff=0;
yoff=1;
}
if(sumBR >= 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];
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;
sumMax=sumBR;
totquad=sumBR;
}
if(sumTR >= 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];
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;
sumMax=sumTR;
totquad=sumTR;
}
totquad=sumMax;
for (int ix=0; ix<2; ix++) {
for (int iy=0; iy<2; iy++) {
sDum[iy][ix] = cluster[iy+yoff][ix+xoff];
}
}
return corner;
@ -235,7 +267,6 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
}
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);
@ -244,6 +275,14 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
}
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) {
@ -294,6 +333,13 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
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]){
@ -306,23 +352,95 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
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,r,t,b;
sum=cl[0]+cl[1]+cl[2]+cl[3]+cl[4]+cl[5]+cl[6]+cl[7]+cl[8];
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[iy+3*ix];
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) {
l=cl[0]+cl[3]+cl[6];
r=cl[2]+cl[5]+cl[8];
b=cl[0]+cl[1]+cl[2];
t=cl[6]+cl[7]+cl[8];
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, double *cl, double &etax, double &etay) {
static int calcMyEta(double totquad, int quad, int *cl, double &etax, double &etay) {
double l,r,t,b, sum;
int yoff;
switch (quad) {
@ -367,6 +485,21 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
}
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;
}