mirror of
https://github.com/slsdetectorgroup/slsDetectorPackage.git
synced 2026-01-16 16:45:55 +01:00
Common mode added, working for gotthard double modules
This commit is contained in:
436
slsDetectorCalibration/interpolations/eta2InterpolationBase.h
Normal file
436
slsDetectorCalibration/interpolations/eta2InterpolationBase.h
Normal file
@@ -0,0 +1,436 @@
|
||||
#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) {
|
||||
// cout << "e2ib " << nb << " " << emin << " " << emax << endl;
|
||||
|
||||
if (nbeta<=0) {
|
||||
nbeta=nSubPixels*10;
|
||||
}
|
||||
if (etamin>=etamax) {
|
||||
etamin=-1;
|
||||
etamax=2;
|
||||
cout << ":" <<endl;
|
||||
}
|
||||
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;
|
||||
};
|
||||
|
||||
eta2InterpolationBase(eta2InterpolationBase *orig): etaInterpolationBase(orig){ };
|
||||
|
||||
virtual eta2InterpolationBase* Clone()=0;/* {
|
||||
|
||||
return new eta2InterpolationBase(this);
|
||||
|
||||
};
|
||||
|
||||
*/
|
||||
|
||||
|
||||
|
||||
//////////////////////////////////////////////////////////////////////////////
|
||||
//////////// /*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];
|
||||
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 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);
|
||||
}
|
||||
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) {
|
||||
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
|
||||
//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 int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay) {
|
||||
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:
|
||||
;
|
||||
}
|
||||
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 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(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;
|
||||
};
|
||||
|
||||
/* 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
Normal file
294
slsDetectorCalibration/interpolations/eta3InterpolationBase.h
Normal file
@@ -0,0 +1,294 @@
|
||||
#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
|
||||
@@ -0,0 +1,417 @@
|
||||
#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
|
||||
Reference in New Issue
Block a user