mirror of
https://github.com/slsdetectorgroup/slsDetectorPackage.git
synced 2025-06-04 00:50:42 +02:00
moench03OnTheFlyAnalysis multipthreaded with interpolation works
This commit is contained in:
parent
d8803ca5e4
commit
10e57319bf
@ -13,8 +13,8 @@ all: moench03OnTheFlyAnalysis
|
||||
|
||||
|
||||
|
||||
moench03OnTheFlyAnalysis: $(MAIN) $(INCS)
|
||||
g++ -o moench03OnTheFlyAnalysis $(MAIN) -lm -ltiff -lstdc++ $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF)
|
||||
moench03OnTheFlyAnalysis: $(MAIN) $(INCS) clean
|
||||
g++ -o moench03OnTheFlyAnalysis $(MAIN) -lm -ltiff -lstdc++ $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL
|
||||
|
||||
clean:
|
||||
rm -f moench03OnTheFlyAnalysis
|
||||
|
||||
@ -29,28 +29,41 @@ class MovingStat
|
||||
/**
|
||||
clears the moving average number of samples parameter, mean and standard deviation
|
||||
*/
|
||||
void Set(double val, double rms=0)
|
||||
void Set(double val, double rms=0, int m=-1)
|
||||
{
|
||||
m_n = n;
|
||||
m_newM=val*n;
|
||||
if (rms<=0)
|
||||
m_newM2=val*val*n;
|
||||
else
|
||||
m_newM2=(n*rms*rms+m_newM*m_newM/n);
|
||||
if (m>=0) m_n = m; else m_n = n;
|
||||
m_newM=val*m_n;
|
||||
SetRMS(rms);
|
||||
}
|
||||
/**
|
||||
clears the moving average number of samples parameter, mean and standard deviation
|
||||
*/
|
||||
void SetRMS(double rms)
|
||||
{
|
||||
if (rms<=0) {
|
||||
m_newM2=m_newM*m_newM/n;
|
||||
m_n=0;
|
||||
} else {
|
||||
if (m_n>0)
|
||||
m_newM2=(m_n*rms*rms+m_newM*m_newM/m_n);
|
||||
else {
|
||||
m_newM2=(m_n*rms*rms+m_newM*m_newM/n);
|
||||
m_n=0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/** sets number of samples parameter
|
||||
\param i number of samples parameter to be set
|
||||
*/
|
||||
|
||||
void SetN(int i) {if (i>=1) n=i;};
|
||||
int SetN(int i) {if (i>=1) n=i; return n;};
|
||||
|
||||
/**
|
||||
gets number of samples parameter
|
||||
\returns actual number of samples parameter
|
||||
*/
|
||||
int GetN() {return n;};
|
||||
int GetN() {return m_n;};
|
||||
|
||||
/** calculates the moving average i.e. adds if number of elements is lower than number of samples parameter, pushes otherwise
|
||||
\param x value to calculate the moving average
|
||||
|
||||
@ -5,6 +5,11 @@
|
||||
#include "slsDetectorData.h"
|
||||
#include "pedestalSubtraction.h"
|
||||
#include "commonModeSubtraction.h"
|
||||
#include "tiffIO.h"
|
||||
#ifndef FRAMEMODE_DEF
|
||||
#define FRAMEMODE_DEF
|
||||
enum frameMode { eFrame, ePedestal, eFlat };
|
||||
#endif
|
||||
|
||||
|
||||
template <class dataType> class analogDetector {
|
||||
@ -13,6 +18,7 @@ template <class dataType> class analogDetector {
|
||||
|
||||
public:
|
||||
|
||||
|
||||
|
||||
/**
|
||||
|
||||
@ -30,7 +36,7 @@ template <class dataType> class analogDetector {
|
||||
|
||||
|
||||
analogDetector(slsDetectorData<dataType> *d, int sign=1,
|
||||
commonModeSubtraction *cm=NULL, int nnx=-1, int nny=-1) : det(d), nx(nnx), ny(nny), stat(NULL), cmSub(cm), iframe(-1), dataSign(sign){
|
||||
commonModeSubtraction *cm=NULL, int nped=1000, int nnx=-1, int nny=-1, double *gm=NULL) : det(d), nx(nnx), ny(nny), stat(NULL), cmSub(cm), iframe(-1), dataSign(sign), gmap(gm) {
|
||||
|
||||
if (det)
|
||||
det->getDetectorSize(nx,ny);
|
||||
@ -38,6 +44,9 @@ template <class dataType> class analogDetector {
|
||||
stat=new pedestalSubtraction*[ny];
|
||||
for (int i=0; i<ny; i++) {
|
||||
stat[i]=new pedestalSubtraction[nx];
|
||||
for (int ix=0; ix<nx; ix++) {
|
||||
stat[i][ix].SetNPedestals(nped);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -45,8 +54,83 @@ template <class dataType> class analogDetector {
|
||||
/**
|
||||
destructor. Deletes the cluster structure and the pdestalSubtraction array
|
||||
*/
|
||||
virtual ~analogDetector() {for (int i=0; i<ny; i++) delete [] stat[i]; delete [] stat;};
|
||||
virtual ~analogDetector() {for (int i=0; i<ny; i++) delete [] stat[i]; delete [] stat; };
|
||||
|
||||
/**
|
||||
clone
|
||||
*/
|
||||
analogDetector(analogDetector* orig) {
|
||||
/* copy construction from orig*/
|
||||
det=orig->det;
|
||||
nx=orig->nx;
|
||||
ny=orig->ny;
|
||||
dataSign=orig->dataSign;
|
||||
iframe=orig->iframe;
|
||||
gmap=orig->gmap;
|
||||
cmSub=orig->cmSub;
|
||||
|
||||
stat=new pedestalSubtraction*[ny];
|
||||
for (int i=0; i<ny; i++) {
|
||||
stat[i]=new pedestalSubtraction[nx];
|
||||
}
|
||||
|
||||
int nped=orig->SetNPedestals();
|
||||
//cout << nped << " " << orig->getPedestal(ix,iy) << orig->getPedestalRMS(ix,iy) << endl;
|
||||
for (int iy=0; iy<ny; iy++) {
|
||||
for (int ix=0; ix<nx; ix++) {
|
||||
stat[iy][ix].SetNPedestals(nped);
|
||||
setPedestal(ix,iy,orig->getPedestal(ix,iy),orig->getPedestalRMS(ix,iy),orig->GetNPedestals(ix,iy));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
virtual analogDetector *Clone() {
|
||||
return new analogDetector(this);
|
||||
}
|
||||
|
||||
int getDataSize(){return det->getDataSize();};
|
||||
|
||||
/**
|
||||
set gain map
|
||||
*/
|
||||
double *setGainMap(double *gm) {gmap=gm; return gmap;};
|
||||
|
||||
/**
|
||||
return gain map
|
||||
*/
|
||||
double *getGainMap() {return gmap;};
|
||||
|
||||
double *readGainMap(const char * imgname) {
|
||||
uint32 nnx, nny;
|
||||
float *gm=ReadFromTiff( imgname, nny, nnx);
|
||||
if (gm) {
|
||||
if (gmap) delete [] gmap;
|
||||
gmap=new double[nnx*nny];
|
||||
for (int ix=0; ix<nnx; ix++) {
|
||||
for (int iy=0; iy<nny; iy++) {
|
||||
gmap[iy*nnx+ix]=gm[iy*nnx+ix];
|
||||
}
|
||||
}
|
||||
return gmap;
|
||||
}
|
||||
return NULL;
|
||||
}
|
||||
|
||||
void *writeGainMap(const char * imgname) {
|
||||
float *gm=NULL;
|
||||
if (gmap) {
|
||||
gm=new float[nx*ny];
|
||||
for (int ix=0; ix<nx; ix++) {
|
||||
for (int iy=0; iy<ny; iy++) {
|
||||
gm[iy*nx+ix]=gmap[iy*nx+ix];
|
||||
}
|
||||
}
|
||||
return WriteToTiff(gm, imgname, ny, nx);
|
||||
}
|
||||
return NULL;
|
||||
}
|
||||
|
||||
|
||||
|
||||
/** resets the pedestalSubtraction array and the commonModeSubtraction */
|
||||
void newDataSet(){iframe=-1; for (int iy=0; iy<ny; iy++) for (int ix=0; ix<nx; ix++) stat[iy][ix].Clear(); if (cmSub) cmSub->Clear(); };
|
||||
@ -108,11 +192,91 @@ template <class dataType> class analogDetector {
|
||||
\param iy pixel y coordinate
|
||||
\param val value to set
|
||||
*/
|
||||
virtual void setPedestal(int ix, int iy, double val, double rms=0){if (ix>=0 && ix<nx && iy>=0 && iy<ny) stat[iy][ix].setPedestal(val,rms);};
|
||||
virtual void setPedestal(int ix, int iy, double val, double rms=0, int m=-1){if (ix>=0 && ix<nx && iy>=0 && iy<ny) stat[iy][ix].setPedestal(val,rms, m);};
|
||||
|
||||
/**
|
||||
sets pedestal
|
||||
\param ix pixel x coordinate
|
||||
\param iy pixel y coordinate
|
||||
\param val value to set
|
||||
*/
|
||||
virtual void setPedestalRMS(int ix, int iy, double rms=0){if (ix>=0 && ix<nx && iy>=0 && iy<ny) stat[iy][ix].setPedestalRMS(rms);};
|
||||
|
||||
|
||||
|
||||
void *writePedestals(const char * imgname) {
|
||||
float *gm=NULL;
|
||||
gm=new float[nx*ny];
|
||||
for (int ix=0; ix<nx; ix++) {
|
||||
for (int iy=0; iy<ny; iy++) {
|
||||
if (cmSub)
|
||||
gm[iy*nx+ix]=stat[iy][ix].getPedestal()-cmSub->getCommonMode();
|
||||
else
|
||||
gm[iy*nx+ix]=stat[iy][ix].getPedestal();
|
||||
}
|
||||
}
|
||||
WriteToTiff(gm, imgname, ny, nx);
|
||||
delete [] gm;
|
||||
return NULL;
|
||||
}
|
||||
|
||||
int readPedestals(const char * imgname) {
|
||||
uint32 nnx, nny;
|
||||
float *gm=ReadFromTiff( imgname, nny, nnx);
|
||||
if (nnx>nx) nnx=nx;
|
||||
if (nny>ny) nny=ny;
|
||||
|
||||
|
||||
|
||||
if (gm) {
|
||||
for (int ix=0; ix<nnx; ix++) {
|
||||
for (int iy=0; iy<nny; iy++) {
|
||||
stat[iy][ix].setPedestal(gm[iy*nx+ix],-1,-1);
|
||||
}
|
||||
}
|
||||
delete [] gm;
|
||||
return 1;
|
||||
}
|
||||
return NULL;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
void *writePedestalRMS(const char * imgname) {
|
||||
float *gm=NULL;
|
||||
gm=new float[nx*ny];
|
||||
for (int ix=0; ix<nx; ix++) {
|
||||
for (int iy=0; iy<ny; iy++) {
|
||||
gm[iy*nx+ix]=stat[iy][ix].getPedestalRMS();
|
||||
}
|
||||
}
|
||||
WriteToTiff(gm, imgname, ny, nx);
|
||||
delete [] gm;
|
||||
return NULL;
|
||||
}
|
||||
|
||||
int readPedestalRMS(const char * imgname) {
|
||||
uint32 nnx, nny;
|
||||
float *gm=ReadFromTiff( imgname, nny, nnx);
|
||||
if (nnx>nx) nnx=nx;
|
||||
if (nny>ny) nny=ny;
|
||||
if (gm) {
|
||||
for (int ix=0; ix<nnx; ix++) {
|
||||
for (int iy=0; iy<nny; iy++) {
|
||||
stat[iy][ix].setPedestalRMS(gm[iy*nx+ix]);
|
||||
}
|
||||
}
|
||||
delete [] gm;
|
||||
return 1;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
virtual void addToPedestal(char *data) {
|
||||
|
||||
@ -137,6 +301,7 @@ template <class dataType> class analogDetector {
|
||||
|
||||
double val;
|
||||
if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
|
||||
|
||||
if (det)
|
||||
val=dataSign*det->getValue(data, ix, iy);
|
||||
else
|
||||
@ -165,30 +330,42 @@ template <class dataType> class analogDetector {
|
||||
return val;
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
virtual double subtractPedestal(char *data, int ix, int iy=0) {
|
||||
double g=1.;
|
||||
if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
|
||||
if (gmap) {
|
||||
g=gmap[iy*nx+ix];
|
||||
if (g==0) g=-1.;
|
||||
}
|
||||
|
||||
if (det)
|
||||
return dataSign*det->getValue(data, ix, iy)-getPedestal(ix,iy);
|
||||
return (dataSign*det->getValue(data, ix, iy)-getPedestal(ix,iy))/g;
|
||||
else
|
||||
return ((double*)data)[iy*nx+ix]-getPedestal(ix,iy);
|
||||
return (((double*)data)[iy*nx+ix]-getPedestal(ix,iy))/g;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
|
||||
virtual int getNPhotons(char *data, int ix, int iy=0, int thr=-1) {
|
||||
|
||||
if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
|
||||
if (thr<=0) thr=-1*thr*getPedestalRMS(ix,iy);
|
||||
double g=1.;
|
||||
if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
|
||||
|
||||
if (gmap) {
|
||||
g=gmap[iy*nx+ix];
|
||||
if (g==0) g=-1.;
|
||||
}
|
||||
|
||||
if (thr<=0) thr=-1*thr*getPedestalRMS(ix,iy)/g;
|
||||
|
||||
if (det)
|
||||
return (dataSign*det->getValue(data, ix, iy)-getPedestal(ix,iy))/thr;
|
||||
return (dataSign*det->getValue(data, ix, iy)-getPedestal(ix,iy))/g/thr;
|
||||
else
|
||||
return (((double*)data)[(iy)*nx+ix]-getPedestal(ix,iy))/thr;
|
||||
return (((double*)data)[(iy)*nx+ix]-getPedestal(ix,iy))/g/thr;
|
||||
|
||||
}
|
||||
return 0;
|
||||
@ -226,6 +403,13 @@ template <class dataType> class analogDetector {
|
||||
return stat[0][0].SetNPedestals();
|
||||
};
|
||||
|
||||
int GetNPedestals(int ix, int iy) {
|
||||
if (ix>=0 && ix<nx && iy>=0 && iy<ny)
|
||||
return stat[iy][ix].GetNPedestals();
|
||||
else
|
||||
return -1;
|
||||
};
|
||||
|
||||
|
||||
double getROI(char *data, int xmin, int xmax, int ymin=0, int ymax=1) {
|
||||
double val=0;
|
||||
@ -235,9 +419,19 @@ template <class dataType> class analogDetector {
|
||||
val+=subtractPedestal(data, ix, iy);
|
||||
return val;
|
||||
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
virtual void processData(char *data, frameMode i=eFrame, int *val=NULL) {
|
||||
switch(i) {
|
||||
case ePedestal:
|
||||
addToPedestal(data);
|
||||
break;
|
||||
default:
|
||||
getNPhotons(data,-1,val);
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
protected:
|
||||
|
||||
@ -248,7 +442,7 @@ template <class dataType> class analogDetector {
|
||||
commonModeSubtraction *cmSub;/**< commonModeSubtraction class */
|
||||
int dataSign; /**< sign of the data i.e. 1 if photon is positive, -1 if negative */
|
||||
int iframe; /**< frame number (not from file but incremented within the dataset every time newFrame is called */
|
||||
|
||||
double *gmap;
|
||||
};
|
||||
|
||||
#endif
|
||||
|
||||
@ -9,6 +9,7 @@
|
||||
#endif
|
||||
|
||||
#include "slsInterpolation.h"
|
||||
#include "tiffIO.h"
|
||||
|
||||
class etaInterpolationBase : public slsInterpolation {
|
||||
|
||||
@ -35,6 +36,38 @@ class etaInterpolationBase : public slsInterpolation {
|
||||
|
||||
};
|
||||
|
||||
etaInterpolationBase(etaInterpolationBase *orig): slsInterpolation(orig){
|
||||
nbeta=orig->nbeta;
|
||||
etamin=orig->etamin;
|
||||
etamax=orig->etamax;
|
||||
|
||||
etastep=(etamax-etamin)/nbeta;
|
||||
#ifdef MYROOT1
|
||||
heta=(TH2D*)(orig->heta)->Clone("heta");
|
||||
hhx=(TH2D*)(orig->hhx)->Clone("hhx");
|
||||
hhy=(TH2D*)(orig->hhy)->Clone("hhy");
|
||||
#endif
|
||||
|
||||
#ifndef MYROOT1
|
||||
heta=new int[nbeta*nbeta];
|
||||
memcpy(heta,orig->heta,nbeta*nbeta*sizeof(int));
|
||||
hhx=new int[nbeta*nbeta];
|
||||
memcpy(hhx,orig->hhx,nbeta*nbeta*sizeof(int));
|
||||
hhy=new int[nbeta*nbeta];
|
||||
memcpy(hhy,orig->hhy,nbeta*nbeta*sizeof(int));
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
};
|
||||
virtual etaInterpolationBase* Clone() {
|
||||
|
||||
return new etaInterpolationBase(this);
|
||||
|
||||
};
|
||||
|
||||
|
||||
|
||||
#ifdef MYROOT1
|
||||
TH2D *setEta(TH2D *h, int nb=-1, double emin=1, double emax=0)
|
||||
{
|
||||
@ -65,12 +98,74 @@ class etaInterpolationBase : public slsInterpolation {
|
||||
}
|
||||
return heta;
|
||||
};
|
||||
|
||||
int *getFlatField(){return setEta(NULL);};
|
||||
|
||||
|
||||
|
||||
|
||||
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=-0.1;
|
||||
etamax=1.1;
|
||||
}
|
||||
|
||||
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 NULL;
|
||||
}
|
||||
nbeta=nnx;
|
||||
if (gm) {
|
||||
if (heta) {
|
||||
delete [] heta;
|
||||
delete [] hhx;
|
||||
delete [] hhy;
|
||||
}
|
||||
|
||||
heta=new int[nbeta*nbeta];
|
||||
hhx=new int[nbeta*nbeta];
|
||||
hhy=new int[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 NULL;
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
virtual void prepareInterpolation(int &ok)=0;
|
||||
virtual void prepareInterpolation(int &ok){};
|
||||
|
||||
/* ////////////////////////////////////////////////////////////////////////////// */
|
||||
|
||||
@ -125,7 +220,7 @@ int *gethhx()
|
||||
|
||||
double xpos_eta,ypos_eta;
|
||||
double dX,dY;
|
||||
double ex,ey;
|
||||
int ex,ey;
|
||||
switch (corner)
|
||||
{
|
||||
case TOP_LEFT:
|
||||
@ -159,13 +254,13 @@ int *gethhx()
|
||||
ex=(etax-etamin)/etastep;
|
||||
ey=(etay-etamin)/etastep;
|
||||
if (ex<0) ex=0;
|
||||
if (ex>=nSubPixels) ex=nSubPixels-1;
|
||||
if (ex>=nbeta) ex=nbeta-1;
|
||||
if (ey<0) ey=0;
|
||||
if (ey>=nSubPixels) ey=nSubPixels-1;
|
||||
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);
|
||||
xpos_eta=(((double)hhx[(ey*nbeta+ex)]))/((double)nSubPixels);
|
||||
ypos_eta=(((double)hhy[(ey*nbeta+ex)]))/((double)nSubPixels);
|
||||
//else
|
||||
//return 0;
|
||||
|
||||
@ -173,6 +268,7 @@ int *gethhx()
|
||||
|
||||
int_x=((double)x) + 0.5*dX + xpos_eta;
|
||||
int_y=((double)y) + 0.5*dY + ypos_eta;
|
||||
// cout << etax << " " << ex << " " << etay << " " << ey << " " << xpos_eta << " " << int_x << " " << ypos_eta << " " << int_y << endl;
|
||||
//return 1;
|
||||
|
||||
}
|
||||
@ -205,11 +301,12 @@ int *gethhx()
|
||||
ey=(eta3y-etamin)/etastep;
|
||||
|
||||
if (ex<0) ex=0;
|
||||
if (ex>=nSubPixels) ex=nSubPixels-1;
|
||||
if (ex>=nbeta) ex=nbeta-1;
|
||||
if (ey<0) ey=0;
|
||||
if (ey>=nSubPixels) ey=nSubPixels-1;
|
||||
xpos_eta=(((double)hhx[(int)(ey*nbeta+ex)]))/((double)nSubPixels);
|
||||
ypos_eta=(((double)hhy[(int)(ey*nbeta+ex)]))/((double)nSubPixels);
|
||||
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;
|
||||
@ -246,6 +343,7 @@ int *gethhx()
|
||||
#ifndef MYROOT1
|
||||
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]++;
|
||||
#endif
|
||||
|
||||
@ -2,16 +2,24 @@
|
||||
#define ETA_INTERPOLATION_POSXY_H
|
||||
|
||||
|
||||
#include "tiffIO.h"
|
||||
#include "etaInterpolationBase.h"
|
||||
|
||||
class etaInterpolationPosXY : public 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(etaInterpolationPosXY *orig): etaInterpolationBase(orig){};
|
||||
|
||||
virtual etaInterpolationPosXY* Clone() {
|
||||
|
||||
return new etaInterpolationPosXY(this);
|
||||
|
||||
};
|
||||
|
||||
virtual void prepareInterpolation(int &ok)
|
||||
{
|
||||
cout <<"?"<< endl;
|
||||
ok=1;
|
||||
#ifdef MYROOT1
|
||||
if (hhx) delete hhx;
|
||||
@ -40,7 +48,6 @@ class etaInterpolationPosXY : public etaInterpolationBase{
|
||||
double hix[nbeta]; //integral of projection x
|
||||
double hiy[nbeta]; //integral of projection y
|
||||
int ii=0;
|
||||
|
||||
for (int ib=0; ib<nbeta; ib++) {
|
||||
|
||||
tot_eta_x=0;
|
||||
@ -58,7 +65,7 @@ class etaInterpolationPosXY : public etaInterpolationBase{
|
||||
|
||||
for (int iby=1; iby<nbeta; iby++) {
|
||||
hix[iby]=hix[iby-1]+hx[iby];
|
||||
hiy[iby]=hiy[iby-1]+hx[iby];
|
||||
hiy[iby]=hiy[iby-1]+hy[iby];
|
||||
}
|
||||
|
||||
ii=0;
|
||||
@ -90,8 +97,32 @@ class etaInterpolationPosXY : public etaInterpolationBase{
|
||||
|
||||
|
||||
}
|
||||
|
||||
return ;
|
||||
|
||||
#ifdef SAVE_ALL
|
||||
char tit[10000];
|
||||
|
||||
float *etah=new float[nbeta*nbeta];
|
||||
int etabins=nbeta;
|
||||
|
||||
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);
|
||||
#endif
|
||||
return ;
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
@ -11,12 +11,20 @@ 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_t x, Int_t y, double *data, double &int_x, double &int_y)
|
||||
virtual void getInterpolatedPosition(int x, int y, double *data, double &int_x, double &int_y)
|
||||
{
|
||||
double sDum[2][2];
|
||||
double tot, totquad;
|
||||
@ -30,9 +38,9 @@ class linearInterpolation : public slsInterpolation{
|
||||
|
||||
return;
|
||||
};
|
||||
|
||||
|
||||
|
||||
virtual void getInterpolatedPosition(Int_t x, Int_t y, double etax, double etay, int corner, double &int_x, double &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;
|
||||
@ -77,7 +85,7 @@ class linearInterpolation : public slsInterpolation{
|
||||
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
virtual void getPositionETA3(Int_t x, Int_t y, double *data, double &int_x, double &int_y)
|
||||
virtual void getPositionETA3(int x, int y, double *data, double &int_x, double &int_y)
|
||||
{
|
||||
double sDum[2][2];
|
||||
double tot, totquad;
|
||||
|
||||
@ -17,12 +17,20 @@
|
||||
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 void getInterpolatedPosition(Int_t x, Int_t y, double *data, double &int_x, double &int_y)
|
||||
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);
|
||||
@ -30,13 +38,13 @@ class noInterpolation : public slsInterpolation{
|
||||
|
||||
return ;
|
||||
};
|
||||
virtual void getInterpolatedPosition(Int_t x, Int_t y, double etax, double etay, int corner, double &int_x, double &int_y)
|
||||
virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int corner, double &int_x, double &int_y)
|
||||
{
|
||||
return getInterpolatedPosition(x, y, NULL, int_x, int_y);
|
||||
};
|
||||
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
virtual void getPositionETA3(Int_t x, Int_t y, double *data, double &int_x, double &int_y)
|
||||
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);
|
||||
|
||||
@ -7,6 +7,7 @@
|
||||
#include <TH2F.h>
|
||||
#endif
|
||||
|
||||
#include "tiffIO.h"
|
||||
#ifndef DEF_QUAD
|
||||
#define DEF_QUAD
|
||||
enum quadrant {
|
||||
@ -25,7 +26,7 @@ class slsInterpolation
|
||||
{
|
||||
|
||||
public:
|
||||
slsInterpolation(int nx=400, int ny=400, int ns=25) :nPixelsX(nx), nPixelsY(ny), nSubPixels(ns) {
|
||||
slsInterpolation(int nx=400, int ny=400, int ns=25) :nPixelsX(nx), nPixelsY(ny), nSubPixels(ns), id(0) {
|
||||
|
||||
#ifdef MYROOT1
|
||||
hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
|
||||
@ -37,7 +38,33 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
|
||||
|
||||
};
|
||||
|
||||
slsInterpolation(slsInterpolation *orig){
|
||||
nPixelsX=orig->nPixelsX;
|
||||
nPixelsY=orig->nPixelsY;
|
||||
nSubPixels=orig->nSubPixels;
|
||||
#ifdef MYROOT1
|
||||
hint=(TH2F*)(orig->hint)->Clone("hint");
|
||||
#endif
|
||||
|
||||
#ifndef MYROOT1
|
||||
hint=new int[nSubPixels*nPixelsX*nSubPixels*nPixelsY];
|
||||
memcpy(hint, orig->hint,nSubPixels*nPixelsX*nSubPixels*nPixelsY*sizeof(int));
|
||||
#endif
|
||||
|
||||
};
|
||||
|
||||
virtual int setId(int i) {id=i; return id;};
|
||||
|
||||
virtual slsInterpolation* Clone() = 0;
|
||||
|
||||
int getNSubPixels() {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.
|
||||
@ -53,18 +80,52 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
|
||||
#ifndef MYROOT1
|
||||
virtual int *getInterpolatedImage(){return hint;};
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
|
||||
void *writeInterpolatedImage(const char * imgname) {
|
||||
cout << "!" <<endl;
|
||||
float *gm=NULL;
|
||||
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]=hint[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;
|
||||
//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;
|
||||
|
||||
|
||||
|
||||
|
||||
#ifdef MYROOT1
|
||||
TH2F *addToImage(double int_x, double int_y){hint->Fill(int_x, int_y); return hint;};
|
||||
#endif
|
||||
|
||||
#ifndef MYROOT1
|
||||
virtual int *addToImage(double int_x, double int_y){ int iy=nSubPixels*int_y; int ix=nSubPixels*int_x;
|
||||
if (ix>=0 && ix<(nPixelsX*nSubPixels) && iy<(nSubPixels*nPixelsY) && iy>=0 )(*(hint+ix+iy*nPixelsX))+=1;
|
||||
if (ix>=0 && ix<(nPixelsX*nSubPixels) && iy<(nSubPixels*nPixelsY) && iy>=0 )(*(hint+ix+iy*nPixelsX*nSubPixels))+=1;
|
||||
return hint;
|
||||
};
|
||||
#endif
|
||||
@ -75,10 +136,14 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
|
||||
|
||||
#ifdef MYROOT1
|
||||
virtual TH2D *getFlatField(){return NULL;};
|
||||
virtual TH2D *setFlatField(TH2D *h){return NULL;};
|
||||
#endif
|
||||
|
||||
#ifndef MYROOT1
|
||||
virtual int *getFlatField(){return NULL;};
|
||||
virtual int *setFlatField(int *h){return NULL;};
|
||||
void *writeFlatField(const char * imgname){return NULL;};
|
||||
void *readFlatField(const char * imgname, int nb=-1, double emin=1, double emax=0){return NULL;};
|
||||
#endif
|
||||
|
||||
//virtual void Streamer(TBuffer &b);
|
||||
@ -165,40 +230,44 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
|
||||
if (totquad>0) {
|
||||
switch(corner) {
|
||||
case TOP_LEFT:
|
||||
t = sDum[1][1] ;
|
||||
t = sDum[1][1];
|
||||
r = sDum[0][1] ;
|
||||
toth=sDum[1][1]+sDum[1][0];
|
||||
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][0]+t;
|
||||
totv=sDum[0][0]+r;
|
||||
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]+t;
|
||||
totv=sDum[0][1]+r;
|
||||
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][1]+t;
|
||||
totv=sDum[0][1]+r;
|
||||
toth=sDum[1][0]+sDum[1][1];
|
||||
totv=sDum[1][0]+sDum[0][0];
|
||||
break;
|
||||
default:
|
||||
etax=-1;
|
||||
etay=-1;
|
||||
etax=-1000;
|
||||
etay=-1000;
|
||||
return 0;
|
||||
}
|
||||
etax=r/totv;
|
||||
etay=t/toth;
|
||||
//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);
|
||||
@ -263,6 +332,7 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
|
||||
#ifndef MYROOT1
|
||||
int *hint;
|
||||
#endif
|
||||
int id;
|
||||
|
||||
};
|
||||
|
||||
|
||||
@ -46,10 +46,50 @@ class interpolatingDetector : public singlePhotonDetector {
|
||||
commonModeSubtraction *cm=NULL,
|
||||
int nped=1000,
|
||||
int nd=100, int nnx=-1, int nny=-1) :
|
||||
singlePhotonDetector(d, 3,nsigma,sign, cm, nped, nd, nnx, nny) , interp(inte) {};
|
||||
singlePhotonDetector(d, 3,nsigma,sign, cm, nped, nd, nnx, nny) , interp(inte), id(0) {};
|
||||
|
||||
|
||||
|
||||
interpolatingDetector(interpolatingDetector *orig) : singlePhotonDetector(orig) {
|
||||
interp=(orig->interp)->Clone();
|
||||
id=orig->id;
|
||||
}
|
||||
|
||||
|
||||
virtual interpolatingDetector *Clone() {
|
||||
return new interpolatingDetector(this);
|
||||
}
|
||||
|
||||
virtual int setId(int i) {id=i; interp->setId(id); return id;};
|
||||
|
||||
virtual void prepareInterpolation(int &ok) {
|
||||
cout << "*"<< endl;
|
||||
#ifdef SAVE_ALL
|
||||
char tit[1000];
|
||||
sprintf(tit,"/scratch/ped_%d.tiff",id);
|
||||
writePedestals(tit);
|
||||
sprintf(tit,"/scratch/ped_rms_%d.tiff",id);
|
||||
writePedestalRMS(tit);
|
||||
if (gmap) {
|
||||
sprintf(tit,"/scratch/gmap_%d.tiff",id);
|
||||
writeGainMap(tit);
|
||||
}
|
||||
#endif
|
||||
interp->prepareInterpolation(ok);
|
||||
}
|
||||
|
||||
int getImageSize(int &nnx, int &nny, int &ns) {return interp->getImageSize(nnx, nny, ns);};
|
||||
#ifdef MYROOT1
|
||||
virtual TH2F *getInterpolatedImage()
|
||||
#endif
|
||||
#ifndef MYROOT1
|
||||
virtual int *getInterpolatedImage()
|
||||
#endif
|
||||
{
|
||||
if (interp)
|
||||
return interp->getInterpolatedImage();
|
||||
}
|
||||
|
||||
#ifdef MYROOT1
|
||||
virtual TH2F *addToInterpolatedImage(char *data, single_photon_hit *clusters, int &nph)
|
||||
#endif
|
||||
@ -79,13 +119,23 @@ class interpolatingDetector : public singlePhotonDetector {
|
||||
return NULL;
|
||||
};
|
||||
|
||||
void *writeInterpolatedImage(const char * imgname) {
|
||||
cout << id << "=" << imgname<< endl;
|
||||
interp->writeInterpolatedImage(imgname);
|
||||
return NULL;
|
||||
}
|
||||
|
||||
|
||||
|
||||
int addFrame(char *data, single_photon_hit *clusters, int ff=0) {
|
||||
|
||||
int addFrame(char *data, single_photon_hit *clusters=NULL, int ff=0) {
|
||||
|
||||
|
||||
single_photon_hit *cl;
|
||||
single_photon_hit clust;
|
||||
if (clusters)
|
||||
cl=clusters;
|
||||
else
|
||||
cl=&clust;
|
||||
|
||||
double int_x,int_y, eta_x, eta_y;
|
||||
int nph=0;
|
||||
double val[ny][nx];
|
||||
@ -118,8 +168,8 @@ class interpolatingDetector : public singlePhotonDetector {
|
||||
|
||||
|
||||
|
||||
|
||||
(clusters+nph)->rms=getPedestalRMS(ix,iy);
|
||||
cl->rms=getPedestalRMS(ix,iy);
|
||||
//(clusters+nph)->rms=getPedestalRMS(ix,iy);
|
||||
|
||||
// cout << iframe << " " << nph << " " << ix << " " << iy << endl;
|
||||
|
||||
@ -147,7 +197,7 @@ class interpolatingDetector : public singlePhotonDetector {
|
||||
v=&(val[iy+ir][ix+ic]);
|
||||
|
||||
if (ir==0 && ic==0) {
|
||||
if (*v<-nSigma*(clusters+nph)->rms) {
|
||||
if (*v<-nSigma*cl->rms) {
|
||||
eventMask[iy][ix]=NEGATIVE_PEDESTAL;
|
||||
// cout << "neg ped" << endl;
|
||||
}
|
||||
@ -173,74 +223,79 @@ class interpolatingDetector : public singlePhotonDetector {
|
||||
}
|
||||
|
||||
if (bl>=br && bl>=tl && bl>=tr) {
|
||||
(clusters+nph)->quad=BOTTOM_LEFT;
|
||||
(clusters+nph)->quadTot=bl;
|
||||
cl->quad=BOTTOM_LEFT;
|
||||
cl->quadTot=bl;
|
||||
xoff=0;
|
||||
yoff=0;
|
||||
} else if (br>=bl && br>=tl && br>=tr) {
|
||||
(clusters+nph)->quad=BOTTOM_RIGHT;
|
||||
(clusters+nph)->quadTot=br;
|
||||
cl->quad=BOTTOM_RIGHT;
|
||||
cl->quadTot=br;
|
||||
xoff=1;
|
||||
yoff=0;
|
||||
} else if (tl>=br && tl>=bl && tl>=tr) {
|
||||
(clusters+nph)->quad=TOP_LEFT;
|
||||
(clusters+nph)->quadTot=tl;
|
||||
cl->quad=TOP_LEFT;
|
||||
cl->quadTot=tl;
|
||||
xoff=0;
|
||||
yoff=1;
|
||||
} else if (tr>=bl && tr>=tl && tr>=br) {
|
||||
(clusters+nph)->quad=TOP_RIGHT;
|
||||
(clusters+nph)->quadTot=tr;
|
||||
cl->quad=TOP_RIGHT;
|
||||
cl->quadTot=tr;
|
||||
xoff=1;
|
||||
yoff=1;
|
||||
}
|
||||
|
||||
if (max>nSigma*(clusters+nph)->rms || tot>sqrt(clusterSizeY*clusterSize)*nSigma*(clusters+nph)->rms || ((clusters+nph)->quadTot)>sqrt(cy*cs)*nSigma*(clusters+nph)->rms) {
|
||||
if (max>nSigma*cl->rms || tot>sqrt(clusterSizeY*clusterSize)*nSigma*cl->rms || (cl->quadTot)>sqrt(cy*cs)*nSigma*cl->rms) {
|
||||
if (val[iy][ix]>=max) {
|
||||
// cout << "max" << endl;
|
||||
eventMask[iy][ix]=PHOTON_MAX;
|
||||
(clusters+nph)->tot=tot;
|
||||
(clusters+nph)->x=ix;
|
||||
(clusters+nph)->y=iy;
|
||||
(clusters+nph)->ped=getPedestal(ix,iy, 0);
|
||||
cl->tot=tot;
|
||||
cl->x=ix;
|
||||
cl->y=iy;
|
||||
cl->ped=getPedestal(ix,iy, 0);
|
||||
// cout << iframe << " " << ix << " " << iy << " "<< (clusters+nph)->tot << " " << (clusters+nph)->quadTot << " " << (clusters+nph)->ped<< " " << (clusters+nph)->rms << endl;
|
||||
for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) {
|
||||
for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) {
|
||||
if ((iy+ir)>=0 && (iy+ir)<ny && (ix+ic)>=0 && (ix+ic)<nx) {
|
||||
(clusters+nph)->set_data(val[iy+ir][ix+ic],ic,ir);
|
||||
cout << val[iy+ir][ix+ic] << " " ;
|
||||
cl->set_data(val[iy+ir][ix+ic],ic,ir);
|
||||
|
||||
|
||||
// cout << val[iy+ir][ix+ic] << " " ;
|
||||
}
|
||||
}
|
||||
cout << endl << " " ;
|
||||
// cout << endl << " " ;
|
||||
}
|
||||
|
||||
cout << endl << " " ;
|
||||
// cout << endl << " " ;
|
||||
|
||||
cc[0][0]=(clusters+nph)->get_data(-1+xoff,-1+yoff);
|
||||
cc[1][0]=(clusters+nph)->get_data(-1+xoff,0+yoff);
|
||||
cc[0][1]=(clusters+nph)->get_data(0+xoff,-1+yoff);
|
||||
cc[1][1]=(clusters+nph)->get_data(0+xoff,0+yoff);
|
||||
cc[0][0]=cl->get_data(-1+xoff,-1+yoff);
|
||||
cc[1][0]=cl->get_data(-1+xoff,0+yoff);
|
||||
cc[0][1]=cl->get_data(0+xoff,-1+yoff);
|
||||
cc[1][1]=cl->get_data(0+xoff,0+yoff);
|
||||
|
||||
cout << cc[0][0] << " " << cc[0][1] << endl;
|
||||
cout << cc[1][0] << " " << cc[1][1] << endl;
|
||||
// cout << cc[0][0] << " " << cc[0][1] << endl;
|
||||
// cout << cc[1][0] << " " << cc[1][1] << endl;
|
||||
|
||||
|
||||
cout << endl << " " ;
|
||||
// cout << endl << " " ;
|
||||
|
||||
if (interp) {
|
||||
interp->calcEta((clusters+nph)->quadTot,cc,eta_x, eta_y);
|
||||
cout << eta_x << " " << eta_y << endl;
|
||||
cout << "**************************************************************************"<< endl;
|
||||
//interp->calcEta((clusters+nph)->quadTot,cc,eta_x, eta_y);
|
||||
interp->calcEtaL(cl->quadTot,cl->quad,cc,eta_x, eta_y);
|
||||
// cout << eta_x << " " << eta_y << endl;
|
||||
// cout << "eta" << endl;
|
||||
if (ff) {
|
||||
interp->addToFlatField(eta_x,eta_y);
|
||||
// cout << "**************************************************************************"<< endl;
|
||||
} else {
|
||||
// cout << "interp" << endl;
|
||||
interp->getInterpolatedPosition(ix,iy,eta_x,eta_y,(clusters+nph)->quad,int_x,int_y);
|
||||
interp->getInterpolatedPosition(ix,iy,eta_x,eta_y,cl->quad,int_x,int_y);
|
||||
// cout << "add" << endl;
|
||||
interp->addToImage(int_x, int_y);
|
||||
}
|
||||
// cout << "done" << endl;
|
||||
}
|
||||
if (clusters) cl=(clusters+nph);
|
||||
|
||||
nph++;
|
||||
} else {
|
||||
// cout << "ph" << endl;
|
||||
@ -259,12 +314,25 @@ class interpolatingDetector : public singlePhotonDetector {
|
||||
};
|
||||
|
||||
|
||||
virtual void processData(char *data, frameMode i=eFrame, int *val=NULL) {
|
||||
switch(i) {
|
||||
case ePedestal:
|
||||
addToPedestal(data);
|
||||
break;
|
||||
case eFlat:
|
||||
addFrame(data,NULL,1);
|
||||
break;
|
||||
default:
|
||||
addFrame(data);
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
protected:
|
||||
|
||||
|
||||
slsInterpolation *interp;
|
||||
|
||||
int id;
|
||||
};
|
||||
|
||||
|
||||
|
||||
@ -188,6 +188,13 @@ class moench03Ctb10GbT1Data : public slsReceiverData<uint16_t> {
|
||||
|
||||
virtual char *readNextFrame(ifstream &filebin, int& ff, int &np) {
|
||||
char *data=new char[packetSize*nPackets];
|
||||
char *d=readNextFrame(filebin, ff, np, data);
|
||||
if (d==NULL) {delete [] data; data=NULL;}
|
||||
return data;
|
||||
|
||||
}
|
||||
|
||||
virtual char *readNextFrame(ifstream &filebin, int& ff, int &np, char *data) {
|
||||
char *retval=0;
|
||||
int nd;
|
||||
int fnum = -1;
|
||||
@ -215,7 +222,7 @@ class moench03Ctb10GbT1Data : public slsReceiverData<uint16_t> {
|
||||
if (getFrameNumber(packet) !=fnum) {
|
||||
|
||||
if (np==0){
|
||||
delete [] data;
|
||||
// delete [] data;
|
||||
return NULL;
|
||||
} else
|
||||
filebin.seekg(-8208,ios_base::cur);
|
||||
@ -240,7 +247,7 @@ class moench03Ctb10GbT1Data : public slsReceiverData<uint16_t> {
|
||||
}
|
||||
|
||||
if (np==0){
|
||||
delete [] data;
|
||||
// delete [] data;
|
||||
return NULL;
|
||||
}
|
||||
|
||||
@ -249,6 +256,23 @@ class moench03Ctb10GbT1Data : public slsReceiverData<uint16_t> {
|
||||
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
int getPacketNumber(int x, int y) {return dataMap[y][x]/8208;};
|
||||
|
||||
|
||||
|
||||
@ -13,7 +13,9 @@
|
||||
|
||||
#include "interpolatingDetector.h"
|
||||
#include "etaInterpolationPosXY.h"
|
||||
|
||||
#include "linearInterpolation.h"
|
||||
#include "noInterpolation.h"
|
||||
#include "multiThreadedDetector.h"
|
||||
|
||||
#include <ctime>
|
||||
|
||||
@ -27,13 +29,14 @@ void *moenchProcessFrame() {
|
||||
char fname[10000];
|
||||
strcpy(fname,"/mnt/moench_data/m03-15_mufocustube/plant_40kV_10uA/m03-15_100V_g4hg_300us_dtf_0.raw");
|
||||
|
||||
|
||||
int nthreads=3;
|
||||
|
||||
int nph, nph1;
|
||||
single_photon_hit clusters[NR*NC];
|
||||
// cout << "hits "<< endl;
|
||||
int etabins=250;
|
||||
double etamin=-0.1, etamax=1.1;
|
||||
int nsubpix=5;
|
||||
int etabins=550;
|
||||
double etamin=-1, etamax=2;
|
||||
int nsubpix=4;
|
||||
float *etah=new float[etabins*etabins];
|
||||
// cout << "etah "<< endl;
|
||||
cout << "image size "<< nsubpix*nsubpix*NC*NR << endl;
|
||||
@ -44,11 +47,15 @@ void *moenchProcessFrame() {
|
||||
// cout << "decoder "<< endl;
|
||||
etaInterpolationPosXY *interp=new etaInterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax);
|
||||
// cout << "interp "<< endl;
|
||||
|
||||
interpolatingDetector *filter=new interpolatingDetector(decoder,interp, 5, 1, 0, 10000, 10000);
|
||||
//linearInterpolation *interp=new linearInterpolation(NC, NR, nsubpix);
|
||||
//noInterpolation *interp=new noInterpolation(NC, NR, nsubpix);
|
||||
interp->readFlatField("/scratch/eta_100.tiff",etamin,etamax);
|
||||
interpolatingDetector *filter=new interpolatingDetector(decoder,interp, 5, 1, 0, 1000, 10);
|
||||
filter->readPedestals("/scratch/ped_100.tiff");
|
||||
cout << "filter "<< endl;
|
||||
|
||||
|
||||
|
||||
char *buff;
|
||||
int nf=0;
|
||||
int ok=0;
|
||||
@ -62,8 +69,10 @@ void *moenchProcessFrame() {
|
||||
filter->newDataSet();
|
||||
|
||||
|
||||
multiThreadedDetector *mt=new multiThreadedDetector(filter,nthreads,100);
|
||||
nph=0;
|
||||
nph1=0;
|
||||
//int np;
|
||||
int iph;
|
||||
|
||||
cout << "file name " << fname << endl;
|
||||
@ -72,55 +81,76 @@ void *moenchProcessFrame() {
|
||||
cout << "Opened file " << fname<< endl;
|
||||
else
|
||||
cout << "Could not open file " << fname<< endl;
|
||||
while ((buff=decoder->readNextFrame(filebin, iFrame)) && nf<5E4) {
|
||||
if (nf<1.1E4)
|
||||
filter->addToFlatField(buff,clusters,nph1);
|
||||
else {
|
||||
if (ok==0) {
|
||||
cout << "**************************************************************************"<< endl;
|
||||
heta=interp->getFlatField();
|
||||
for (int ii=0; ii<etabins*etabins; ii++) {
|
||||
etah[ii]=(float)heta[ii];
|
||||
}
|
||||
|
||||
|
||||
std::time(&end_time);
|
||||
cout << std::ctime(&end_time) << " " << nf << " " << nph1 << " " << nph << endl;
|
||||
WriteToTiff(etah, "/scratch/eta.tiff", etabins, etabins);
|
||||
mt->setFrameMode(eFrame);
|
||||
mt->prepareInterpolation(ok);
|
||||
mt->StartThreads();
|
||||
mt->popFree(buff);
|
||||
|
||||
while ((decoder->readNextFrame(filebin, iFrame, np, buff)) && nf<1.5E4) {
|
||||
if (nf<9E3)
|
||||
;
|
||||
else {
|
||||
|
||||
// if (nf>1.1E4 && ok==0) {
|
||||
// mt->prepareInterpolation(ok);
|
||||
// mt->setFrameMode(eFrame);
|
||||
// //ok=1;
|
||||
// }
|
||||
|
||||
interp->prepareInterpolation(ok);
|
||||
cout << "**************************************************************************"<< endl;
|
||||
std::time(&end_time);
|
||||
cout << std::ctime(&end_time) << " " << nf << " " << nph1 << " " << nph << endl;
|
||||
}
|
||||
filter->addToInterpolatedImage(buff,clusters,nph1);
|
||||
}
|
||||
mt->pushData(buff);
|
||||
mt->nextThread();
|
||||
// cout << " " << (void*)buff;
|
||||
mt->popFree(buff);
|
||||
|
||||
// if (ok==0) {
|
||||
// cout << "**************************************************************************"<< endl;
|
||||
// heta=interp->getFlatField();
|
||||
// // for (int ii=0; ii<etabins*etabins; ii++) {
|
||||
// // etah[ii]=(float)heta[ii];
|
||||
// // }
|
||||
|
||||
|
||||
// std::time(&end_time);
|
||||
// cout << std::ctime(&end_time) << " " << nf << endl;
|
||||
// // WriteToTiff(etah, "/scratch/eta.tiff", etabins, etabins);
|
||||
|
||||
// interp->prepareInterpolation(ok);
|
||||
// cout << "**************************************************************************"<< endl;
|
||||
// std::time(&end_time);
|
||||
// cout << std::ctime(&end_time) << " " << nf << endl;
|
||||
// }
|
||||
// filter->processData(buff,eFrame);
|
||||
// }
|
||||
|
||||
nph+=nph1;
|
||||
// nph+=nph1;
|
||||
}
|
||||
if (nf%1000==0) {
|
||||
std::time(&end_time);
|
||||
cout << std::ctime(&end_time) << " " << nf << " " << nph1 << " " << nph << endl;
|
||||
cout << std::ctime(&end_time) << " " << nf << endl;
|
||||
}
|
||||
|
||||
nf++;
|
||||
delete [] buff;
|
||||
//delete [] buff;
|
||||
iFrame=-1;
|
||||
}
|
||||
}
|
||||
|
||||
if (filebin.is_open())
|
||||
filebin.close();
|
||||
else
|
||||
cout << "could not open file " << fname << endl;
|
||||
if (filebin.is_open())
|
||||
filebin.close();
|
||||
else
|
||||
cout << "could not open file " << fname << endl;
|
||||
|
||||
|
||||
himage=interp->getInterpolatedImage();
|
||||
for (int ii=0; ii<nsubpix*nsubpix*NC*NR; ii++) {
|
||||
image[ii]=(float)himage[ii];
|
||||
}
|
||||
WriteToTiff(image, "/scratch/int_image.tiff", nsubpix*NR, nsubpix*NC);
|
||||
delete interp;
|
||||
delete decoder;
|
||||
cout << "Read " << nf << " frames" << endl;
|
||||
return NULL;
|
||||
mt->StopThreads();
|
||||
|
||||
char tit[10000];
|
||||
sprintf(tit,"/scratch/int_image_mt%d.tiff",nthreads);
|
||||
|
||||
mt->writeInterpolatedImage(tit);
|
||||
// delete [] etah;
|
||||
|
||||
// delete interp;
|
||||
//delete decoder;
|
||||
//cout << "Read " << nf << " frames" << endl;
|
||||
return NULL;
|
||||
}
|
||||
|
||||
int main(int argc, char *argv[]){
|
||||
|
||||
427
slsDetectorCalibration/multiThreadedDetector.h
Normal file
427
slsDetectorCalibration/multiThreadedDetector.h
Normal file
@ -0,0 +1,427 @@
|
||||
#include <vector>
|
||||
#include <string>
|
||||
#include <sstream>
|
||||
#include <iomanip>
|
||||
#include <fstream>
|
||||
#include <stdio.h>
|
||||
//#include <deque>
|
||||
//#include <list>
|
||||
//#include <queue>
|
||||
#include <fstream>
|
||||
|
||||
#include <pthread.h>
|
||||
|
||||
#include "interpolatingDetector.h"
|
||||
#include "circularFifo.h"
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
class threadedDetector
|
||||
{
|
||||
public:
|
||||
threadedDetector(interpolatingDetector *d, int fs=10000) {
|
||||
char *mem, *mm;
|
||||
det=d;
|
||||
fifoFree=new CircularFifo<char>(fs);
|
||||
fifoData=new CircularFifo<char>(fs);
|
||||
mem=(char*)malloc(fs*det->getDataSize());
|
||||
for (int i=0; i<fs; i++) {
|
||||
mm=mem+i*det->getDataSize();
|
||||
fifoFree->push(mm);
|
||||
}
|
||||
stop=1;
|
||||
fMode=eFrame;
|
||||
}
|
||||
|
||||
|
||||
int setFrameMode(int fm) {if (fMode>=0) fMode=fm; return fMode;}
|
||||
|
||||
void prepareInterpolation(int &ok) {
|
||||
cout << "-" << endl;
|
||||
det->prepareInterpolation(ok);
|
||||
};
|
||||
|
||||
int *getInterpolatedImage() {
|
||||
return det->getInterpolatedImage();
|
||||
}
|
||||
int getImageSize(int &nnx, int &nny, int &ns) {return det->getImageSize(nnx, nny, ns);};
|
||||
virtual ~threadedDetector() {StopThread(); free(mem); delete fifoFree; delete fifoData;}
|
||||
|
||||
/** Returns true if the thread was successfully started, false if there was an error starting the thread */
|
||||
bool StartThread()
|
||||
{ stop=0;
|
||||
return (pthread_create(&_thread, NULL, processData, this) == 0);
|
||||
}
|
||||
|
||||
void StopThread()
|
||||
{ stop=1;
|
||||
(void) pthread_join(_thread, NULL);
|
||||
}
|
||||
|
||||
|
||||
bool pushData(char* &ptr) {
|
||||
fifoData->push(ptr);
|
||||
}
|
||||
bool popFree(char* &ptr) {
|
||||
fifoFree->pop(ptr);
|
||||
}
|
||||
|
||||
//protected:
|
||||
/** Implement this method in your subclass with the code you want your thread to run. */
|
||||
//virtual void InternalThreadEntry() = 0;
|
||||
void *writeInterpolatedImage(const char * imgname) {cout << "a" <<endl; return det->writeInterpolatedImage(imgname);};
|
||||
private:
|
||||
interpolatingDetector *det;
|
||||
int fMode;
|
||||
int *dataSize;
|
||||
pthread_t _thread;
|
||||
char *mem;
|
||||
CircularFifo<char> *fifoFree;
|
||||
CircularFifo<char> *fifoData;
|
||||
int stop;
|
||||
char *data;
|
||||
|
||||
static void * processData(void * ptr) {
|
||||
threadedDetector *This=((threadedDetector *)ptr);
|
||||
return This->processData();
|
||||
}
|
||||
|
||||
void * processData() {
|
||||
while (!stop) {
|
||||
if (fifoData->isEmpty()) {
|
||||
usleep(100);
|
||||
} else {
|
||||
fifoData->pop(data); //blocking!
|
||||
det->processData(data,(frameMode)fMode);
|
||||
fifoFree->push(data);
|
||||
}
|
||||
}
|
||||
return NULL;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
};
|
||||
|
||||
|
||||
|
||||
class multiThreadedDetector
|
||||
{
|
||||
public:
|
||||
multiThreadedDetector(interpolatingDetector *d, int n, int fs=1000) : nThreads(n), ithread(0) {
|
||||
dd[0]=d;
|
||||
if (nThreads==1)
|
||||
dd[0]->setId(100);
|
||||
else
|
||||
dd[0]->setId(0);
|
||||
for (int i=1; i<nThreads; i++) {
|
||||
dd[i]=d->Clone();
|
||||
dd[i]->setId(i);
|
||||
}
|
||||
|
||||
for (int i=0; i<nThreads; i++) {
|
||||
dets[i]=new threadedDetector(dd[i], fs);
|
||||
}
|
||||
|
||||
int nnx, nny, ns;
|
||||
int nn=dets[0]->getImageSize(nnx, nny,ns);
|
||||
image=new int[nn];
|
||||
|
||||
}
|
||||
|
||||
~multiThreadedDetector() {
|
||||
StopThreads();
|
||||
for (int i=0; i<nThreads; i++)
|
||||
delete dets[i];
|
||||
for (int i=1; i<nThreads; i++)
|
||||
delete dd[i];
|
||||
delete [] image;
|
||||
}
|
||||
|
||||
|
||||
int setFrameMode(int fm) { int ret; for (int i=0; i<nThreads; i++) ret=dets[i]->setFrameMode(fm); return ret;};
|
||||
|
||||
void prepareInterpolation(int &ok) {
|
||||
int oo;
|
||||
ok=1;
|
||||
for (int i=0; i<nThreads; i++) {
|
||||
cout << i << endl;
|
||||
dets[i]->prepareInterpolation(oo);
|
||||
//if (oo<1) ok=0;
|
||||
}
|
||||
};
|
||||
|
||||
int *getInterpolatedImage() {
|
||||
int *img;
|
||||
int nnx, nny, ns;
|
||||
int nn=dets[0]->getImageSize(nnx, nny, ns);
|
||||
//for (i=0; i<nn; i++) image[i]=0;
|
||||
|
||||
for (int ii=0; ii<nThreads; ii++) {
|
||||
img=dets[ii]->getInterpolatedImage();
|
||||
for (int i=0; i<nn; i++) {
|
||||
if (ii==0)
|
||||
image[i]=img[i];
|
||||
else
|
||||
image[i]+=img[i];
|
||||
}
|
||||
|
||||
}
|
||||
return image;
|
||||
|
||||
}
|
||||
|
||||
void *writeInterpolatedImage(const char * imgname) {
|
||||
#ifdef SAVE_ALL
|
||||
for (int ii=0; ii<nThreads; ii++) {
|
||||
char tit[10000];cout << "m" <<endl;
|
||||
sprintf(tit,"/scratch/int_%d.tiff",ii);
|
||||
dets[ii]->writeInterpolatedImage(tit);
|
||||
}
|
||||
#endif
|
||||
getInterpolatedImage();
|
||||
int nnx, nny, ns;
|
||||
int nn=dets[0]->getImageSize(nnx, nny, ns);
|
||||
float *gm=new float[ nn];
|
||||
if (gm) {
|
||||
for (int ix=0; ix<nn; ix++) {
|
||||
gm[ix]=image[ix];
|
||||
}
|
||||
|
||||
WriteToTiff(gm,imgname ,nnx, nny);
|
||||
delete [] gm;
|
||||
} else cout << "Could not allocate float image " << endl;
|
||||
return NULL;
|
||||
}
|
||||
|
||||
|
||||
void StartThreads() {
|
||||
for (int i=0; i<nThreads; i++)
|
||||
dets[i]->StartThread();
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
void StopThreads() {
|
||||
for (int i=0; i<nThreads; i++)
|
||||
dets[i]->StopThread();
|
||||
|
||||
}
|
||||
|
||||
|
||||
bool pushData(char* &ptr) {
|
||||
dets[ithread]->pushData(ptr);
|
||||
}
|
||||
|
||||
bool popFree(char* &ptr) {
|
||||
dets[ithread]->popFree(ptr);
|
||||
}
|
||||
|
||||
int nextThread() {
|
||||
ithread++;
|
||||
if (ithread==nThreads) ithread=0;
|
||||
return ithread;
|
||||
}
|
||||
|
||||
private:
|
||||
bool stop;
|
||||
const int nThreads;
|
||||
threadedDetector *dets[20];
|
||||
interpolatingDetector *dd[20];
|
||||
int ithread;
|
||||
int *image;
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
/* void *moenchProcessFrame() { */
|
||||
/* char fname[10000]; */
|
||||
/* strcpy(fname,"/mnt/moench_data/m03-15_mufocustube/plant_40kV_10uA/m03-15_100V_g4hg_300us_dtf_0.raw"); */
|
||||
|
||||
|
||||
/* int nph, nph1; */
|
||||
/* single_photon_hit clusters[NR*NC]; */
|
||||
/* // cout << "hits "<< endl; */
|
||||
/* int etabins=550; */
|
||||
/* double etamin=-1, etamax=2; */
|
||||
/* int nsubpix=2; */
|
||||
/* float *etah=new float[etabins*etabins]; */
|
||||
/* // cout << "etah "<< endl; */
|
||||
/* cout << "image size "<< nsubpix*nsubpix*NC*NR << endl; */
|
||||
/* float *image=new float[nsubpix*nsubpix*NC*NR]; */
|
||||
/* int *heta, *himage; */
|
||||
|
||||
/* moench03Ctb10GbT1Data *decoder=new moench03Ctb10GbT1Data(); */
|
||||
/* // cout << "decoder "<< endl; */
|
||||
/* etaInterpolationPosXY *interp=new etaInterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax); */
|
||||
/* // cout << "interp "<< endl; */
|
||||
/* // linearInterpolation *interp=new linearInterpolation(NC, NR, nsubpix); */
|
||||
|
||||
/* interpolatingDetector *filter=new interpolatingDetector(decoder,interp, 5, 1, 0, 1000, 100); */
|
||||
/* cout << "filter "<< endl; */
|
||||
|
||||
|
||||
/* char *buff; */
|
||||
/* int nf=0; */
|
||||
/* int ok=0; */
|
||||
/* ifstream filebin; */
|
||||
/* std::time_t end_time; */
|
||||
|
||||
/* int iFrame=-1; */
|
||||
/* int np=-1; */
|
||||
|
||||
|
||||
/* filter->newDataSet(); */
|
||||
|
||||
|
||||
/* nph=0; */
|
||||
/* nph1=0; */
|
||||
/* int iph; */
|
||||
|
||||
/* cout << "file name " << fname << endl; */
|
||||
/* filebin.open((const char *)(fname), ios::in | ios::binary); */
|
||||
/* if (filebin.is_open()) */
|
||||
/* cout << "Opened file " << fname<< endl; */
|
||||
/* else */
|
||||
/* cout << "Could not open file " << fname<< endl; */
|
||||
/* while ((buff=decoder->readNextFrame(filebin, iFrame)) && nf<1.5E4) { */
|
||||
/* if (nf<9E3) */
|
||||
/* ; */
|
||||
/* else if (nf<1.1E4) { */
|
||||
/* filter->processData(buff,eFlat); */
|
||||
/* } else { */
|
||||
/* if (ok==0) { */
|
||||
/* cout << "**************************************************************************"<< endl; */
|
||||
/* heta=interp->getFlatField(); */
|
||||
/* // for (int ii=0; ii<etabins*etabins; ii++) { */
|
||||
/* // etah[ii]=(float)heta[ii]; */
|
||||
/* // } */
|
||||
|
||||
|
||||
/* std::time(&end_time); */
|
||||
/* cout << std::ctime(&end_time) << " " << nf << endl; */
|
||||
/* // WriteToTiff(etah, "/scratch/eta.tiff", etabins, etabins); */
|
||||
|
||||
/* interp->prepareInterpolation(ok); */
|
||||
/* cout << "**************************************************************************"<< endl; */
|
||||
/* std::time(&end_time); */
|
||||
/* cout << std::ctime(&end_time) << " " << nf << endl; */
|
||||
/* } */
|
||||
/* filter->processData(buff,eFrame); */
|
||||
/* } */
|
||||
|
||||
/* nph+=nph1; */
|
||||
|
||||
/* if (nf%1000==0) { */
|
||||
/* std::time(&end_time); */
|
||||
/* cout << std::ctime(&end_time) << " " << nf << endl; */
|
||||
/* } */
|
||||
|
||||
/* nf++; */
|
||||
/* delete [] buff; */
|
||||
/* iFrame=-1; */
|
||||
/* } */
|
||||
|
||||
/* if (filebin.is_open()) */
|
||||
/* filebin.close(); */
|
||||
/* else */
|
||||
/* cout << "could not open file " << fname << endl; */
|
||||
|
||||
|
||||
/* himage=interp->getInterpolatedImage(); */
|
||||
/* for (int ii=0; ii<nsubpix*nsubpix*NC*NR; ii++) { */
|
||||
/* image[ii]=(float)himage[ii]; */
|
||||
/* } */
|
||||
/* WriteToTiff(image, "/scratch/int_image.tiff", nsubpix*NR, nsubpix*NC); */
|
||||
/* delete [] etah; */
|
||||
/* delete [] image; */
|
||||
/* delete interp; */
|
||||
/* delete decoder; */
|
||||
/* cout << "Read " << nf << " frames" << endl; */
|
||||
/* return NULL; */
|
||||
/* } */
|
||||
|
||||
/* int main(int argc, char *argv[]){ */
|
||||
|
||||
/* moenchProcessFrame(); */
|
||||
|
||||
/* } */
|
||||
@ -37,10 +37,19 @@ class pedestalSubtraction {
|
||||
\param i number of elements for the moving average. If -1 (default) or negative, gets.
|
||||
\returns actual number of samples for the moving average
|
||||
*/
|
||||
virtual int SetNPedestals(int i=-1) {if (i>0) stat.SetN(i); return stat.GetN();};
|
||||
virtual int SetNPedestals(int i=-1) {return stat.SetN(i);};
|
||||
|
||||
/**sets/gets the number of samples for the moving average
|
||||
\returns actual number of samples for the moving average
|
||||
*/
|
||||
virtual int GetNPedestals() {return stat.GetN();};
|
||||
|
||||
/** sets the moving average */
|
||||
virtual void setPedestal(double val, double rms=0) {stat.Set(val, rms);}
|
||||
virtual void setPedestal(double val, double rms=0, int m=-1) {stat.Set(val, rms, m);}
|
||||
|
||||
|
||||
/** sets the moving average */
|
||||
virtual void setPedestalRMS(double rms) {stat.SetRMS(rms);}
|
||||
|
||||
|
||||
|
||||
|
||||
@ -57,7 +57,7 @@ public analogDetector<uint16_t> {
|
||||
int sign=1,
|
||||
commonModeSubtraction *cm=NULL,
|
||||
int nped=1000,
|
||||
int nd=100, int nnx=-1, int nny=-1) : analogDetector<uint16_t>(d, sign, cm, nnx, nny), nDark(nd), eventMask(NULL),nSigma (nsigma), clusterSize(csize), clusterSizeY(csize), cluster(NULL), quad(UNDEFINED_QUADRANT), tot(0), quadTot(0) {
|
||||
int nd=100, int nnx=-1, int nny=-1, double *gm=NULL) : analogDetector<uint16_t>(d, sign, cm, nped, nnx, nny, gm), nDark(nd), eventMask(NULL),nSigma (nsigma), clusterSize(csize), clusterSizeY(csize), cluster(NULL), quad(UNDEFINED_QUADRANT), tot(0), quadTot(0) {
|
||||
|
||||
|
||||
|
||||
@ -65,9 +65,6 @@ public analogDetector<uint16_t> {
|
||||
eventMask=new eventType*[ny];
|
||||
for (int i=0; i<ny; i++) {
|
||||
eventMask[i]=new eventType[nx];
|
||||
for (int ix=0; ix<nx; ix++) {
|
||||
stat[i][ix].SetNPedestals(nped);
|
||||
}
|
||||
}
|
||||
|
||||
if (ny==1)
|
||||
@ -80,7 +77,7 @@ public analogDetector<uint16_t> {
|
||||
/**
|
||||
destructor. Deletes the cluster structure and the pdestalSubtraction array
|
||||
*/
|
||||
virtual ~singlePhotonDetector() {delete cluster;};
|
||||
virtual ~singlePhotonDetector() {delete cluster; for (int i=0; i<ny; i++) delete [] eventMask[i]; delete [] eventMask; };
|
||||
|
||||
|
||||
|
||||
@ -89,9 +86,32 @@ public analogDetector<uint16_t> {
|
||||
|
||||
|
||||
|
||||
singlePhotonDetector(singlePhotonDetector *orig) : analogDetector<uint16_t>(orig) {
|
||||
|
||||
nDark=orig->nDark;
|
||||
|
||||
|
||||
eventMask=new eventType*[ny];
|
||||
for (int i=0; i<ny; i++) {
|
||||
eventMask[i]=new eventType[nx];
|
||||
|
||||
}
|
||||
|
||||
|
||||
nSigma=orig->nSigma;
|
||||
clusterSize=orig->clusterSize;
|
||||
clusterSizeY=orig->clusterSizeY;
|
||||
cluster=new single_photon_hit(clusterSize,clusterSizeY);
|
||||
quad=UNDEFINED_QUADRANT;
|
||||
tot=0;
|
||||
quadTot=0;
|
||||
|
||||
}
|
||||
|
||||
|
||||
virtual singlePhotonDetector *Clone() {
|
||||
return new singlePhotonDetector(this);
|
||||
}
|
||||
/** sets/gets number of rms threshold to detect photons
|
||||
\param n number of sigma to be set (0 or negative gets)
|
||||
\returns actual number of sigma parameter
|
||||
@ -118,22 +138,23 @@ public analogDetector<uint16_t> {
|
||||
|
||||
|
||||
|
||||
virtual int *getNPhotons(char *data, double thr=-1) {
|
||||
virtual int *getNPhotons(char *data, double thr=-1, int *nph=NULL) {
|
||||
|
||||
double val;
|
||||
int *nph=new int[nx*ny];
|
||||
if (nph==NULL)
|
||||
nph=new int[nx*ny];
|
||||
double rest[ny][nx];
|
||||
int cy=(clusterSizeY+1)/2;
|
||||
int cs=(clusterSize+1)/2;
|
||||
|
||||
|
||||
double g=-1.;
|
||||
|
||||
int ccs=clusterSize;
|
||||
int ccy=clusterSizeY;
|
||||
|
||||
|
||||
double tthr=thr;
|
||||
int nn;
|
||||
int nn=0;
|
||||
double max=0, tl=0, tr=0, bl=0,br=0, v;
|
||||
|
||||
|
||||
@ -148,13 +169,23 @@ public analogDetector<uint16_t> {
|
||||
for (int ix=0; ix<nx; ix++) {
|
||||
for (int iy=0; iy<ny; iy++) {
|
||||
|
||||
if (thr<=0) tthr=nSigma*getPedestalRMS(ix,iy);
|
||||
|
||||
val=subtractPedestal(data,ix,iy);
|
||||
if (thr<=0) {
|
||||
|
||||
if (gmap) {
|
||||
g=gmap[iy*nx+ix];
|
||||
if (g==0) g=-1.;
|
||||
}
|
||||
|
||||
nph[ix+nx*iy]=analogDetector<uint16_t>::getNPhotons(data,ix,iy,tthr);
|
||||
tthr=nSigma*getPedestalRMS(ix,iy)/g;
|
||||
|
||||
}
|
||||
|
||||
|
||||
val=subtractPedestal(data,ix,iy);
|
||||
nn=analogDetector<uint16_t>::getNPhotons(data,ix,iy,tthr);
|
||||
nph[ix+nx*iy]+=nn;
|
||||
|
||||
rest[iy][ix]=subtractPedestal(data,ix,iy)-nph[ix+nx*iy]*tthr;
|
||||
rest[iy][ix]=(val-nn*tthr);
|
||||
|
||||
|
||||
}
|
||||
@ -185,7 +216,7 @@ public analogDetector<uint16_t> {
|
||||
if (ir>=0 && ic<=0)
|
||||
tl+=v;
|
||||
if (ir>=0 && ic>=0)
|
||||
tr+=v;
|
||||
tr+=v;
|
||||
|
||||
if (v>max) {
|
||||
max=v;
|
||||
@ -205,7 +236,7 @@ public analogDetector<uint16_t> {
|
||||
quadTot=bl;
|
||||
} else if (br>=bl && br>=tl && br>=tr) {
|
||||
quad=BOTTOM_RIGHT;
|
||||
quadTot=br;
|
||||
quadTot=br;
|
||||
} else if (tl>=br && tl>=bl && tl>=tr) {
|
||||
quad=TOP_LEFT;
|
||||
quadTot=tl;
|
||||
@ -519,6 +550,16 @@ int getClusters(char *data, single_photon_hit *clusters) {
|
||||
|
||||
#endif
|
||||
|
||||
virtual void processData(char *data, frameMode i=eFrame, int *val=NULL) {
|
||||
switch(i) {
|
||||
case ePedestal:
|
||||
addToPedestal(data);
|
||||
break;
|
||||
default:
|
||||
getNPhotons(data,-1,val);
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
protected:
|
||||
|
||||
|
||||
@ -32,13 +32,13 @@ void *WriteToTiff(float * imgData, const char * imgname, int nrow, int ncol){
|
||||
int sampleperpixel=1;
|
||||
// unsigned char * buff=NULL;
|
||||
tsize_t linebytes;
|
||||
|
||||
cout << "--" <<endl;
|
||||
TIFF * tif = TIFFOpen(imgname,"w");
|
||||
if (tif) {
|
||||
TIFFSetField(tif,TIFFTAG_IMAGEWIDTH,ncol);
|
||||
TIFFSetField(tif, TIFFTAG_IMAGELENGTH, nrow);
|
||||
TIFFSetField(tif, TIFFTAG_SAMPLESPERPIXEL,sampleperpixel);
|
||||
TIFFSetField(tif, TIFFTAG_BITSPERSAMPLE, sizeof(float)*8);
|
||||
TIFFSetField(tif, TIFFTAG_BITSPERSAMPLE, 32);
|
||||
TIFFSetField(tif, TIFFTAG_ORIENTATION, ORIENTATION_BOTLEFT);
|
||||
TIFFSetField(tif, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);
|
||||
TIFFSetField(tif, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_MINISBLACK);
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user