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

fixed bug in pedestal cloning and changed cluster finder file format- tested only for 3x3

Browse Source
This commit is contained in:
bergamaschi 2020-03-24 12:03:29 +01:00
parent a108a4df4c
commit ed20e17c3a
9 changed files with 349 additions and 244 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
examples/moench03_T1.config
View File

@ -7,28 +7,28 @@ hostname bchip181+
############################################# #############################################
### edit with hostname or 1Gbs IP address of your server ### edit with hostname or 1Gbs IP address of your server
############################################ ############################################
rx_hostname mpc2608 rx_hostname mpc2011
rx_tcpport 1954 rx_tcpport 1954
############################################# #############################################
### edit with 10 Gbs IP of your server ### edit with 10 Gbs IP of your server
############################################ ############################################
udp_dstip 10.1.2.117 udp_dstip 10.1.1.102
############################################# #############################################
### edit with any number in the subnet of your server (first 3 numbers as above) ### edit with any number in the subnet of your server (first 3 numbers as above)
############################################ ############################################
udp_srcip 10.1.2.19 udp_srcip 10.1.1.19
udp_dstport 32411 udp_dstport 32411
############################################# #############################################
### edit with 10 Gbs IP of your server ### edit with 10 Gbs IP of your server
############################################ ############################################
rx_zmqip 10.1.2.117 rx_zmqip 10.1.1.102
rx_zmqport 50003 rx_zmqport 50003
############################################# #############################################
### edit with 1 Gbs IP of PC where you will run the GUI ### edit with 1 Gbs IP of PC where you will run the GUI
############################################ ############################################
zmqip 129.129.202.86 zmqip 129.129.202.136
zmqport 50001 zmqport 50001
@ -39,8 +39,8 @@ rx_datastream 1
frames 1 frames 100000
period 0.001 period 0.0006
############################################# #############################################
### edit with directory you want to write to ### edit with directory you want to write to
############################################ ############################################
@ -48,8 +48,8 @@ fpath /mnt/moench_data/scratch/
fwrite 0 fwrite 0
rx_datastream 1 rx_datastream 1
rx_jsonpara frameMode newpedestal rx_jsonpara frameMode frame
rx_jsonpara detectorMode analog rx_jsonpara detectorMode counting
rx_discardpolicy discardpartial rx_discardpolicy discardpartial
@ -59,3 +59,5 @@ powerchip 1
vhighvoltage 90 vhighvoltage 90
#adcreg 0x14 0x40 #adcreg 0x14 0x40

45
slsDetectorCalibration/MovingStat.h
View File

@ -14,27 +14,27 @@ class MovingStat
/** constructor /** constructor
\param nn number of samples parameter to be used \param nn number of samples parameter to be used
*/ */
MovingStat(int nn=1000) : n(nn), m_n(0) {} MovingStat(int nn) : n(nn), m_n(0), m_newM(0),m_newM2(0) {}
// void setPointers(double *me, double *va) {mean=me; var=va;}
/** /**
clears the moving average number of samples parameter, mean and standard deviation clears the moving average number of samples parameter, mean and standard deviation
*/ */
void Clear() void Clear()
{ {
m_n = 0; m_n = 0;
m_newM=0; m_newM=0;
m_newM2=0; m_newM2=0;
} }
/** /**
clears the moving average number of samples parameter, mean and standard deviation clears the moving average number of samples parameter, mean and standard deviation
*/ */
void Set(double val, double rms=0, int m=-1) void Set(double val, double rms=0, int m=-1)
{ {
if (m>0) m_n = m; else m_n = n; if (m>=0) m_n = m; else m_n = n;
m_newM=val*m_n; m_newM=val*m_n;
// cout << "set " << val << " " << m << " " << m_n << " " << m_newM << endl;
SetRMS(rms); SetRMS(rms);
// cout << "set " << val << " " << m << " " << m_n << " " << m_newM << endl;
} }
/** /**
clears the moving average number of samples parameter, mean and standard deviation clears the moving average number of samples parameter, mean and standard deviation
@ -42,12 +42,15 @@ class MovingStat
void SetRMS(double rms) void SetRMS(double rms)
{ {
if (rms<=0) { if (rms<=0) {
m_newM2=m_newM*m_newM/n; if (m_n>0)
m_newM2=m_newM*m_newM/m_n;
else
m_newM2=0;
//m_n=0; //m_n=0;
} else { } else {
if (m_n>0) if (m_n>0) {
m_newM2=(m_n*rms*rms+m_newM*m_newM/m_n); m_newM2=(m_n*rms*rms+m_newM*m_newM/m_n);
else { } else {
m_newM2=(m_n*rms*rms+m_newM*m_newM/n); m_newM2=(m_n*rms*rms+m_newM*m_newM/n);
m_n=0; m_n=0;
} }
@ -102,13 +105,14 @@ class MovingStat
m_newM = x; m_newM = x;
m_newM2 = x*x; m_newM2 = x*x;
m_n++; m_n++;
} else { }
m_newM = m_newM + x - m_newM/m_n; else {
m_newM2 = m_newM2 + x*x - m_newM2/m_n; m_newM = m_newM + x - m_newM/m_n;
} m_newM2 = m_newM2 + x*x - m_newM2/m_n;
}
} }
/** returns the current number of elements of the moving average /** returns the current number of elements of the moving average
\returns returns the current number of elements of the moving average \returns returns the current number of elements of the moving average
*/ */
@ -122,7 +126,9 @@ class MovingStat
inline double Mean() const inline double Mean() const
{ {
// cout << "get " << m_n << " " << m_newM << " " << m_newM/m_n << endl; // cout << "get " << m_n << " " << m_newM << " " << m_newM/m_n << endl;
return (m_n > 0) ? m_newM/m_n : 0.0;
return (m_n > 0) ? m_newM/m_n : 0.0;
} }
/** returns the squared mean, 0 if no elements are inside /** returns the squared mean, 0 if no elements are inside
@ -138,7 +144,7 @@ class MovingStat
*/ */
inline double Variance() const inline double Variance() const
{ {
return ( (m_n > 1) ? (M2()-Mean()*Mean()) : 0.0 ); return (m_n > 0) ? m_newM2/m_n-m_newM/m_n*m_newM/m_n : 0.0;
} }
/** returns the standard deviation, 0 if no elements are inside /** returns the standard deviation, 0 if no elements are inside
@ -146,7 +152,8 @@ class MovingStat
*/ */
inline double StandardDeviation() const inline double StandardDeviation() const
{ {
return ( (Variance() > 0) ? sqrt( Variance() ) : -1 );
return sqrt(Variance());//
} }
private: private:

134
slsDetectorCalibration/analogDetector.h
View File

@ -68,10 +68,15 @@ template <class dataType> class analogDetector {
det->getDetectorSize(nx,ny); det->getDetectorSize(nx,ny);
stat=new pedestalSubtraction*[ny]; stat=new pedestalSubtraction*[ny];
/* pedMean=new double*[ny]; */
/* pedVariance=new double*[ny]; */
for (int i=0; i<ny; i++) { for (int i=0; i<ny; i++) {
stat[i]=new pedestalSubtraction[nx]; stat[i]=new pedestalSubtraction[nx];
for (int ix=0; ix<nx; ix++) { /* pedMean[i]=new double[nx]; */
/* pedVariance[i]=new double[nx]; */
for (ix=0; ix<nx; ++ix) {
stat[i][ix].SetNPedestals(nped); stat[i][ix].SetNPedestals(nped);
/* stat[i][ix].setPointers(&(pedMean[iy][ix]),&(pedVariance[iy][ix])); */
} }
} }
image=new int[nx*ny]; image=new int[nx*ny];
@ -95,7 +100,16 @@ template <class dataType> class analogDetector {
/** /**
destructor. Deletes the pdestalSubtraction array and the image destructor. Deletes the pdestalSubtraction array and the image
*/ */
virtual ~analogDetector() {for (int i=0; i<ny; i++) delete [] stat[i]; delete [] stat; delete [] image; virtual ~analogDetector() {
for (int i=0; i<ny; i++) {
delete [] stat[i];
/* delete [] pedMean[i]; */
/* delete [] pedVariance[i]; */
};
/* delete [] pedMean; */
/* delete [] pedVariance; */
delete [] stat;
delete [] image;
#ifdef ROOTSPECTRUM #ifdef ROOTSPECTRUM
delete hs; delete hs;
#ifdef ROOTCLUST #ifdef ROOTCLUST
@ -132,16 +146,23 @@ template <class dataType> class analogDetector {
stat=new pedestalSubtraction*[ny]; stat=new pedestalSubtraction*[ny];
/* pedMean=new double*[ny]; */
/* pedVariance=new double*[ny]; */
for (int i=0; i<ny; i++) { for (int i=0; i<ny; i++) {
stat[i]=new pedestalSubtraction[nx]; stat[i]=new pedestalSubtraction[nx];
/* pedMean[i]=new double[nx]; */
/* pedVariance[i]=new double[nx]; */
} }
int nped=orig->SetNPedestals(); int nped=orig->SetNPedestals();
//cout << nped << " " << orig->getPedestal(ix,iy) << orig->getPedestalRMS(ix,iy) << endl; //cout << nped << " " << orig->getPedestal(ix,iy) << orig->getPedestalRMS(ix,iy) << endl;
for (int iy=0; iy<ny; iy++) { for (iy=0; iy<ny; ++iy) {
for (int ix=0; ix<nx; ix++) { for (ix=0; ix<nx; ++ix) {
//stat[iy][ix].setPointers(&(pedMean[iy][ix]),&(pedVariance[iy][ix]));
stat[iy][ix].SetNPedestals(nped); stat[iy][ix].SetNPedestals(nped);
setPedestal(ix,iy,orig->getPedestal(ix,iy),orig->getPedestalRMS(ix,iy),orig->GetNPedestals(ix,iy)); setPedestal(ix,iy,orig->getPedestal(ix,iy),orig->getPedestalRMS(ix,iy),orig->GetNPedestals(ix,iy));
/* if (ix==50 && iy==50) */
/* cout << "clone ped " << " " << ix << " " << iy << " " << getPedestal(ix,iy) << " " << getPedestalRMS(ix,iy)<< " " << GetNPedestals(ix,iy) << endl; */
} }
} }
image=new int[nx*ny]; image=new int[nx*ny];
@ -231,8 +252,8 @@ template <class dataType> class analogDetector {
if (gm) { if (gm) {
if (gmap) delete [] gmap; if (gmap) delete [] gmap;
gmap=new double[nnx*nny]; gmap=new double[nnx*nny];
for (int iy=0; iy<nny; iy++) { for (iy=0; iy<nny; ++iy) {
for (int ix=0; ix<nnx; ix++) { for (ix=0; ix<nnx; ++ix) {
gmap[iy*nnx+ix]=gm[iy*nnx+ix]; gmap[iy*nnx+ix]=gm[iy*nnx+ix];
// cout << gmap[iy*nnx+ix] << " " ; // cout << gmap[iy*nnx+ix] << " " ;
} }
@ -251,8 +272,8 @@ template <class dataType> class analogDetector {
void *ret; void *ret;
if (gmap) { if (gmap) {
gm=new float[nx*ny]; gm=new float[nx*ny];
for (int iy=0; iy<ny; iy++) { for (iy=0; iy<ny; ++iy) {
for (int ix=0; ix<nx; ix++) { for (ix=0; ix<nx; ++ix) {
gm[iy*nx+ix]=gmap[iy*nx+ix]; gm[iy*nx+ix]=gmap[iy*nx+ix];
} }
} }
@ -267,8 +288,8 @@ template <class dataType> class analogDetector {
virtual void newDataSet(){ virtual void newDataSet(){
iframe=-1; iframe=-1;
for (int iy=0; iy<ny; iy++) for (iy=0; iy<ny; ++iy)
for (int ix=0; ix<nx; ix++) { for (ix=0; ix<nx; ++ix) {
stat[iy][ix].Clear(); stat[iy][ix].Clear();
image[iy*nx+ix]=0; image[iy*nx+ix]=0;
} }
@ -285,15 +306,17 @@ template <class dataType> class analogDetector {
}; };
/** resets the commonModeSubtraction and increases the frame index */ /** resets the commonModeSubtraction and increases the frame index */
virtual void newFrame(){iframe++; if (cmSub) cmSub->newFrame();}; virtual void newFrame(){iframe++; if (cmSub) cmSub->newFrame(); det->newFrame();};
/** resets the commonModeSubtraction and increases the frame index */ /** resets the commonModeSubtraction and increases the frame index */
virtual void newFrame(char *data){ virtual void newFrame(char *data){
iframe++; iframe++;
if (cmSub) cmSub->newFrame(); if (cmSub) cmSub->newFrame();
det->newFrame();
// det->getData(data);
calcGhost(data); calcGhost(data);
// cout << getId() << " Calc ghost " << getGhost(15,15) << endl; // cout << getId() << " Calc ghost " << getGhost(15,15) << endl;
}; };
/** sets the commonModeSubtraction algorithm to be used /** sets the commonModeSubtraction algorithm to be used
@ -359,8 +382,8 @@ template <class dataType> class analogDetector {
// cout << "+"<< getId() << endl; // cout << "+"<< getId() << endl;
if (cmSub) { if (cmSub) {
//cout << "*" << endl; //cout << "*" << endl;
for (int iy=ymin; iy<ymax; iy++) { for (iy=ymin; iy<ymax; ++iy) {
for (int ix=xmin; ix<xmax; ix++) { for (ix=xmin; ix<xmax; ++ix) {
// if (getNumpedestals(ix,iy)>0) // if (getNumpedestals(ix,iy)>0)
// if (det->isGood(ix,iy)) { // if (det->isGood(ix,iy)) {
addToCommonMode(data, ix, iy); addToCommonMode(data, ix, iy);
@ -403,8 +426,10 @@ template <class dataType> class analogDetector {
if (ix>=0 && ix<nx && iy>=0 && iy<ny) { if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
if (cmSub && cm>0) { if (cmSub && cm>0) {
return stat[iy][ix].getPedestal()+getCommonMode(ix,iy); return stat[iy][ix].getPedestal()+getCommonMode(ix,iy);
//return pedMean[iy][ix]+getCommonMode(ix,iy);
} else return stat[iy][ix].getPedestal(); }
//return pedMean[iy][ix];
return stat[iy][ix].getPedestal();
} else return -1; } else return -1;
}; };
@ -422,7 +447,7 @@ template <class dataType> class analogDetector {
g=gmap[iy*nx+ix]; g=gmap[iy*nx+ix];
if (g==0) g=-1.; if (g==0) g=-1.;
} }
// return sqrt(pedVariance[iy][ix])/g;
return stat[iy][ix].getPedestalRMS()/g;//divide by gain? return stat[iy][ix].getPedestalRMS()/g;//divide by gain?
} }
return -1; return -1;
@ -444,8 +469,8 @@ template <class dataType> class analogDetector {
if (ped==NULL) { if (ped==NULL) {
ped=new double[nx*ny]; ped=new double[nx*ny];
} }
for (int iy=0; iy<ny; iy++) { for (iy=0; iy<ny; ++iy) {
for (int ix=0; ix<nx; ix++) { for (ix=0; ix<nx; ++ix) {
ped[iy*nx+ix]=stat[iy][ix].getPedestal(); ped[iy*nx+ix]=stat[iy][ix].getPedestal();
//cout << ped[iy*nx+ix] << " " ; //cout << ped[iy*nx+ix] << " " ;
} }
@ -463,8 +488,8 @@ template <class dataType> class analogDetector {
if (ped==NULL) { if (ped==NULL) {
ped=new double[nx*ny]; ped=new double[nx*ny];
} }
for (int iy=0; iy<ny; iy++) { for (iy=0; iy<ny; ++iy) {
for (int ix=0; ix<nx; ix++) { for (ix=0; ix<nx; ++ix) {
ped[iy*nx+ix]=stat[iy][ix].getPedestalRMS(); ped[iy*nx+ix]=stat[iy][ix].getPedestalRMS();
} }
} }
@ -503,8 +528,8 @@ template <class dataType> class analogDetector {
*/ */
virtual void setPedestal(double *ped, double *rms=NULL, int m=-1){ virtual void setPedestal(double *ped, double *rms=NULL, int m=-1){
double rr=0; double rr=0;
for (int iy=ymin; iy<ymax; iy++) { for (iy=ymin; iy<ymax; ++iy) {
for (int ix=xmin; ix<xmax; ix++) { for (ix=xmin; ix<xmax; ++ix) {
if (rms) rr=rms[iy*nx+ix]; if (rms) rr=rms[iy*nx+ix];
stat[iy][ix].setPedestal(ped[iy*nx+ix],rr, m); stat[iy][ix].setPedestal(ped[iy*nx+ix],rr, m);
// cout << ix << " " << iy << " " << ped[iy*nx+ix] << " " << stat[iy][ix].getPedestal() << endl; // cout << ix << " " << iy << " " << ped[iy*nx+ix] << " " << stat[iy][ix].getPedestal() << endl;
@ -533,8 +558,8 @@ template <class dataType> class analogDetector {
\param rms pointer to array of pedestal rms \param rms pointer to array of pedestal rms
*/ */
virtual void setPedestalRMS(double *rms){ virtual void setPedestalRMS(double *rms){
for (int iy=ymin; iy<ymax; iy++) { for (iy=ymin; iy<ymax; ++iy) {
for (int ix=xmin; ix<xmax; ix++) { for (ix=xmin; ix<xmax; ++ix) {
stat[iy][ix].setPedestalRMS(rms[iy*nx+ix]); stat[iy][ix].setPedestalRMS(rms[iy*nx+ix]);
}; };
}; };
@ -574,8 +599,8 @@ template <class dataType> class analogDetector {
#endif #endif
gm=new float[nx*ny]; gm=new float[nx*ny];
for (int iy=0; iy<ny; iy++) { for (iy=0; iy<ny; ++iy) {
for (int ix=0; ix<nx; ix++) { for (ix=0; ix<nx; ++ix) {
gm[iy*nx+ix]=image[iy*nx+ix]; gm[iy*nx+ix]=image[iy*nx+ix];
#ifdef ROOTSPECTRUM #ifdef ROOTSPECTRUM
hmap->SetBinContent(ix+1, iy+1,image[iy*nx+ix]); hmap->SetBinContent(ix+1, iy+1,image[iy*nx+ix]);
@ -622,8 +647,8 @@ template <class dataType> class analogDetector {
TH2F *hmap=new TH2F("hmap","hmap",nx, -0.5,nx-0.5, ny, -0.5, ny-0.5); TH2F *hmap=new TH2F("hmap","hmap",nx, -0.5,nx-0.5, ny, -0.5, ny-0.5);
#endif #endif
for (int iy=0; iy<ny; iy++) { for (iy=0; iy<ny; ++iy) {
for (int ix=0; ix<nx; ix++) { for (ix=0; ix<nx; ++ix) {
/* if (cmSub) */ /* if (cmSub) */
/* gm[iy*nx+ix]=stat[iy][ix].getPedestal()-cmSub->getCommonMode(); */ /* gm[iy*nx+ix]=stat[iy][ix].getPedestal()-cmSub->getCommonMode(); */
/* else */ /* else */
@ -669,8 +694,8 @@ template <class dataType> class analogDetector {
if (gm) { if (gm) {
for (int iy=0; iy<nny; iy++) { for (iy=0; iy<nny; ++iy) {
for (int ix=0; ix<nnx; ix++) { for (ix=0; ix<nnx; ++ix) {
stat[iy][ix].setPedestal(gm[iy*nx+ix],-1,-1); stat[iy][ix].setPedestal(gm[iy*nx+ix],-1,-1);
} }
} }
@ -694,8 +719,8 @@ template <class dataType> class analogDetector {
if (gm) { if (gm) {
for (int iy=0; iy<nny; iy++) { for (iy=0; iy<nny; ++iy) {
for (int ix=0; ix<nnx; ix++) { for (ix=0; ix<nnx; ++ix) {
image[iy*nx+ix]=gm[iy*nx+ix]; image[iy*nx+ix]=gm[iy*nx+ix];
} }
} }
@ -720,9 +745,9 @@ template <class dataType> class analogDetector {
float *gm=NULL; float *gm=NULL;
void *ret; void *ret;
gm=new float[nx*ny]; gm=new float[nx*ny];
for (int iy=0; iy<ny; iy++) { for (iy=0; iy<ny; ++iy) {
for (int ix=0; ix<nx; ix++) { for (ix=0; ix<nx; ++ix) {
gm[iy*nx+ix]=stat[iy][ix].getPedestalRMS(); gm[iy*nx+ix]=stat[iy][ix].getPedestalRMS();
} }
} }
ret=WriteToTiff(gm, imgname, ny, nx); ret=WriteToTiff(gm, imgname, ny, nx);
@ -742,8 +767,8 @@ template <class dataType> class analogDetector {
if (nnx>nx) nnx=nx; if (nnx>nx) nnx=nx;
if (nny>ny) nny=ny; if (nny>ny) nny=ny;
if (gm) { if (gm) {
for (int iy=0; iy<nny; iy++) { for (iy=0; iy<nny; ++iy) {
for (int ix=0; ix<nnx; ix++) { for (ix=0; ix<nnx; ++ix) {
stat[iy][ix].setPedestalRMS(gm[iy*nx+ix]); stat[iy][ix].setPedestalRMS(gm[iy*nx+ix]);
} }
} }
@ -779,11 +804,13 @@ template <class dataType> class analogDetector {
//cout << xmin << " " << xmax << endl; //cout << xmin << " " << xmax << endl;
// cout << ymin << " " << ymax << endl; // cout << ymin << " " << ymax << endl;
for (int iy=ymin; iy<ymax; iy++) { for (iy=ymin; iy<ymax; ++iy) {
for (int ix=xmin; ix<xmax; ix++) { for (ix=xmin; ix<xmax; ++ix) {
if (det->isGood(ix,iy)) { if (det->isGood(ix,iy)) {
// addToPedestal(data,ix,iy,1); // addToPedestal(data,ix,iy,1);
addToPedestal(data,ix,iy,cm); addToPedestal(data,ix,iy,cm);
/* if (ix==50 && iy==50) */
/* cout<< "*ped* " << id << " " << ix << " " << iy << " " << det->getChannel(data,ix,iy) << " " << stat[iy][ix].getPedestal() << " " << stat[iy][ix].getPedestalRMS() <<" " << stat[iy][ix]. GetNPedestals() << endl; */
//if (ix==10 && iy==10) //if (ix==10 && iy==10)
// cout <<ix << " " << iy << " " << getPedestal(ix,iy)<< endl; // cout <<ix << " " << iy << " " << getPedestal(ix,iy)<< endl;
#ifdef ROOTSPECTRUM #ifdef ROOTSPECTRUM
@ -907,8 +934,8 @@ template <class dataType> class analogDetector {
//calcGhost(data); //calcGhost(data);
for (int iy=ymin; iy<ymax; iy++) { for (iy=ymin; iy<ymax; ++iy) {
for (int ix=xmin; ix<xmax; ix++) { for (ix=xmin; ix<xmax; ++ix) {
if (det->isGood(ix,iy)) if (det->isGood(ix,iy))
val[iy*nx+ix]+=subtractPedestal(data, ix, iy,cm); val[iy*nx+ix]+=subtractPedestal(data, ix, iy,cm);
} }
@ -954,8 +981,8 @@ template <class dataType> class analogDetector {
hs->Fill(val,(iy-ymin)*(xmax-xmin)+(ix-xmin)); hs->Fill(val,(iy-ymin)*(xmax-xmin)+(ix-xmin));
#ifdef ROOTCLUST #ifdef ROOTCLUST
double v3=0,v5=0,v7=0,v9=0; double v3=0,v5=0,v7=0,v9=0;
for (int iix=-4; iix<5; iix++) for (int iix=-4; iix<5; i++ix)
for (int iiy=-4; iiy<5; iiy++) { for (int iiy=-4; iiy<5; i++iy) {
if (det) if (det)
val= (dataSign*det->getValue(data, ix+iix, iy+iiy)-getPedestal(ix+iix,iy+iiy,cm))/g; val= (dataSign*det->getValue(data, ix+iix, iy+iiy)-getPedestal(ix+iix,iy+iiy,cm))/g;
else else
@ -1049,8 +1076,8 @@ template <class dataType> class analogDetector {
//calcGhost(data); //calcGhost(data);
addToCommonMode(data); addToCommonMode(data);
for (int iy=ymin; iy<ymax; iy++) { for (iy=ymin; iy<ymax; ++iy) {
for (int ix=xmin; ix<xmax; ix++) { for (ix=xmin; ix<xmax; ++ix) {
if (det->isGood(ix,iy)) if (det->isGood(ix,iy))
nph[iy*nx+ix]+=getNPhotons(data, ix, iy); nph[iy*nx+ix]+=getNPhotons(data, ix, iy);
} }
@ -1064,8 +1091,8 @@ template <class dataType> class analogDetector {
*/ */
virtual void clearImage(){ virtual void clearImage(){
for (int iy=0; iy<ny; iy++) { for (iy=0; iy<ny; ++iy) {
for (int ix=0; ix<nx; ix++) { for (ix=0; ix<nx; ++ix) {
image[iy*nx+ix]=0; image[iy*nx+ix]=0;
} }
} }
@ -1095,8 +1122,8 @@ template <class dataType> class analogDetector {
int SetNPedestals(int i=-1) { int SetNPedestals(int i=-1) {
int ix=0, iy=0; int ix=0, iy=0;
if (i>0) if (i>0)
for (iy=0; iy<ny; iy++) for (iy=0; iy<ny; ++iy)
for (ix=0; ix<nx; ix++) for (ix=0; ix<nx; ++ix)
stat[iy][ix].SetNPedestals(i); stat[iy][ix].SetNPedestals(i);
return stat[0][0].SetNPedestals(); return stat[0][0].SetNPedestals();
}; };
@ -1131,8 +1158,8 @@ template <class dataType> class analogDetector {
newFrame(data); newFrame(data);
//calcGhost(data); //calcGhost(data);
addToCommonMode(data); addToCommonMode(data);
for (int iy=ymi; iy<yma; iy++) for (iy=ymi; iy<yma; ++iy)
for (int ix=xmi; ix<xma; ix++) for (ix=xmi; ix<xma; ++ix)
if (det->isGood(ix,iy)) { if (det->isGood(ix,iy)) {
if (ix>=0 && ix<nx && iy>=0 && iy<ny) if (ix>=0 && ix<nx && iy>=0 && iy<ny)
val+=getNPhotons(data, ix, iy); val+=getNPhotons(data, ix, iy);
@ -1235,6 +1262,8 @@ FILE *getFilePointer(){return myFile;};
int nx; /**< Size of the detector in x direction */ int nx; /**< Size of the detector in x direction */
int ny; /**< Size of the detector in y direction */ int ny; /**< Size of the detector in y direction */
pedestalSubtraction **stat; /**< pedestalSubtraction class */ pedestalSubtraction **stat; /**< pedestalSubtraction class */
/* double **pedMean; /\**< pedestalSubtraction class *\/ */
/* double **pedVariance; /\**< pedestalSubtraction class *\/ */
commonModeSubtraction *cmSub;/**< commonModeSubtraction class */ commonModeSubtraction *cmSub;/**< commonModeSubtraction class */
int dataSign; /**< sign of the data i.e. 1 if photon is positive, -1 if negative */ 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 */ int iframe; /**< frame number (not from file but incremented within the dataset every time newFrame is called */
@ -1252,6 +1281,7 @@ FILE *getFilePointer(){return myFile;};
frameMode fMode; /**< current detector frame mode */ frameMode fMode; /**< current detector frame mode */
detectorMode dMode; /**< current detector frame mode */ detectorMode dMode; /**< current detector frame mode */
FILE *myFile; /**< file pointer to write to */ FILE *myFile; /**< file pointer to write to */
int ix, iy;
#ifdef ROOTSPECTRUM #ifdef ROOTSPECTRUM
TH2F *hs; TH2F *hs;
#ifdef ROOTCLUST #ifdef ROOTCLUST

68
slsDetectorCalibration/dataStructures/slsDetectorData.h
View File

@ -20,6 +20,8 @@ class slsDetectorData {
int **dataROIMask; /**< Array of size nx*ny 1 if channel is good (or in the ROI), 0 if bad channel (or out of ROI) */ int **dataROIMask; /**< Array of size nx*ny 1 if channel is good (or in the ROI), 0 if bad channel (or out of ROI) */
int *xmap; int *xmap;
int *ymap; int *ymap;
dataType **orderedData;
int isOrdered;
public: public:
@ -37,34 +39,27 @@ class slsDetectorData {
\param dROI Array of size nx*ny. The elements are 1s if the channel is good or in the ROI, 0 is bad or out of the ROI. NULL (default) means all 1s. \param dROI Array of size nx*ny. The elements are 1s if the channel is good or in the ROI, 0 is bad or out of the ROI. NULL (default) means all 1s.
*/ */
slsDetectorData(int npx, int npy, int dsize, int **dMap=NULL, dataType **dMask=NULL, int **dROI=NULL): nx(npx), ny(npy), dataSize(dsize) { slsDetectorData(int npx, int npy, int dsize, int **dMap=NULL, dataType **dMask=NULL, int **dROI=NULL): nx(npx), ny(npy), dataSize(dsize), orderedData(NULL), isOrdered(0) {
int el=dsize/sizeof(dataType); int el=dsize/sizeof(dataType);
xmap=new int[el]; xmap=new int[el];
ymap=new int[el]; ymap=new int[el];
// if (dataMask==NULL) {
dataMask=new dataType*[ny];
for(int i = 0; i < ny; i++) {
dataMask[i] = new dataType[nx];
}
// }
// if (dataMap==NULL) { orderedData=new dataType*[ny];
dataMap=new int*[ny]; dataMap=new int*[ny];
for(int i = 0; i < ny; i++) { dataMask=new dataType*[ny];
dataMap[i] = new int[nx];
}
// }
// if (dataROIMask==NULL) {
dataROIMask=new int*[ny]; dataROIMask=new int*[ny];
for(int i = 0; i < ny; i++) { for(int i = 0; i < ny; i++) {
dataMap[i] = new int[nx];
orderedData[i]=new dataType[nx];
dataMask[i] = new dataType[nx];
dataROIMask[i] = new int[nx]; dataROIMask[i] = new int[nx];
for (int j=0; j<nx; j++) for (int j=0; j<nx; j++)
dataROIMask[i][j]=1; dataROIMask[i][j]=1;
}
// } }
for (int ip=0; ip<el; ip++){ for (int ip=0; ip<el; ip++){
xmap[ip]=-1; xmap[ip]=-1;
@ -75,22 +70,24 @@ class slsDetectorData {
setDataMask(dMask); setDataMask(dMask);
setDataROIMask(dROI); setDataROIMask(dROI);
}; }
virtual ~slsDetectorData() { virtual ~slsDetectorData() {
for(int i = 0; i < ny; i++) { for(int i = 0; i < ny; i++) {
delete [] dataMap[i]; delete [] dataMap[i];
delete [] dataMask[i]; delete [] dataMask[i];
delete [] dataROIMask[i]; delete [] dataROIMask[i];
delete [] orderedData[i];
} }
delete [] dataMap; delete [] dataMap;
delete [] dataMask; delete [] dataMask;
delete [] dataROIMask; delete [] dataROIMask;
delete [] orderedData;
delete [] xmap; delete [] xmap;
delete [] ymap; delete [] ymap;
}; };
virtual int getPointer(int ix,int iy) {return dataMap[iy][ix];};
/** /**
defines the data map (as offset) - no error checking if datasize and offsets are compatible! defines the data map (as offset) - no error checking if datasize and offsets are compatible!
\param dMap array of size nx*ny storing the pointers to the data in the dataset (as offset). If NULL (default),the data are arranged as if read out row by row (dataMap[iy][ix]=(iy*nx+ix)*sizeof(dataType);) \param dMap array of size nx*ny storing the pointers to the data in the dataset (as offset). If NULL (default),the data are arranged as if read out row by row (dataMap[iy][ix]=(iy*nx+ix)*sizeof(dataType);)
@ -203,27 +200,32 @@ class slsDetectorData {
int setDataSize(int d) {dataSize=d; return dataSize;}; int setDataSize(int d) {dataSize=d; return dataSize;};
virtual void getPixel(int ip, int &x, int &y) {x=xmap[ip]; y=ymap[ip];}; virtual void getPixel(int ip, int &x, int &y) {x=xmap[ip]; y=ymap[ip];};
virtual dataType **getData(char *ptr, int dsize=-1) { virtual dataType **getData(char *ptr, int dsize=-1) {
int el=dsize/sizeof(dataType); int el=dsize/sizeof(dataType);
dataType **data; //dataType **data;
int ix,iy; int ix,iy;
data=new dataType*[ny]; //data=new dataType*[ny];
for(int i = 0; i < ny; i++) { //for(int i = 0; i < ny; i++) {
data[i]=new dataType[nx]; // data[i]=new dataType[nx];
} //}
isOrdered=0;
if (dsize<=0 || dsize>dataSize) dsize=dataSize; if (dsize<=0 || dsize>dataSize) dsize=dataSize;
for (int ip=0; ip<(el); ip++) { for (int ip=0; ip<(el); ip++) {
getPixel(ip,ix,iy); getPixel(ip,ix,iy);
if (ix>=0 && ix<nx && iy>=0 && iy<ny) { if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
data[iy][ix]=getChannel(ptr,ix,iy); //data[iy][ix]=getChannel(ptr,ix,iy);
orderedData[iy][ix]=*(ptr+ip);//getChannel(ptr,ix,iy);
} }
} }
return data; isOrdered=1;
return orderedData;
}; }
void newFrame(){isOrdered=0;};
virtual double **getImage(char *ptr, int dsize=-1) { virtual double **getImage(char *ptr, int dsize=-1) {
double **data; double **data;
@ -240,8 +242,7 @@ class slsDetectorData {
data[iy][ix]=getValue(ptr,ix,iy); data[iy][ix]=getValue(ptr,ix,iy);
} }
} }
return data; return data;
}; };
/** /**
@ -258,7 +259,10 @@ class slsDetectorData {
// cout << ix << " " << iy << " " ; // cout << ix << " " << iy << " " ;
//cout << dataMap[ix][iy] << " " << (void*)data << " " << dataSize<< endl; //cout << dataMap[ix][iy] << " " << (void*)data << " " << dataSize<< endl;
m=dataMask[iy][ix]; m=dataMask[iy][ix];
d=*((dataType*)(data+dataMap[iy][ix])); if (isOrdered==0)
d=*((dataType*)(data+getPointer(ix,iy)));
else
d=orderedData[iy][ix];
} }
return d^m; return d^m;
}; };

12
slsDetectorCalibration/moenchExecutables/Makefile.zmq
View File

@ -1,18 +1,18 @@
INCDIR= -I. -I../dataStructures ../tiffIO.cpp -I../ -I../interpolations/ -I../../slsSupportLib/include/ -I../../slsReceiverSoftware/include/ -I../../libs/rapidjson/ INCDIR= -I. -I../dataStructures ../tiffIO.cpp -I../ -I../interpolations/ -I../../slsSupportLib/include/ -I../../slsReceiverSoftware/include/ -I../../libs/rapidjson/
LDFLAG= -L/usr/lib64/ -lpthread -lm -lstdc++ -lzmq -pthread -lrt -ltiff -O3 -g -std=c++11 -Wall -L../../build/bin/ -lSlsSupport LDFLAG= -L/usr/lib64/ -lpthread -lm -lstdc++ -lzmq -pthread -lrt -ltiff -O3 -std=c++11 -Wall -L../../build/bin/ -lSlsSupport
#-L../../bin -lhdf5 -L. #-L../../bin -lhdf5 -L.
#DESTDIR?=../bin #DESTDIR?=../bin
all: moenchZmqProcess moenchZmqProcessCtbGui all: moenchZmqProcess
#moenchZmqProcessCtbGui
moenchZmqProcess: moenchZmqProcess.cpp clean moenchZmqProcess: moenchZmqProcess.cpp clean
g++ -o moenchZmqProcess moenchZmqProcess.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DNEWZMQ -DINTERP g++ -o moenchZmqProcess moenchZmqProcess.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DNEWZMQ -DINTERP
moenchZmqProcessCtbGui: moenchZmqProcess.cpp clean #moenchZmqProcessCtbGui: moenchZmqProcess.cpp clean
g++ -o moenchZmqProcessCtbGui moenchZmqProcess.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DNEWZMQ -DINTERP -DCTBGUI # g++ -o moenchZmqProcessCtbGui moenchZmqProcess.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DNEWZMQ -DINTERP -DCTBGUI
clean: clean:
rm -f moenchZmqProcess rm -f moenchZmqProcess

50
slsDetectorCalibration/moenchExecutables/moenchZmqProcess.cpp
View File

@ -1,6 +1,6 @@
#define WRITE_QUAD //#define WRITE_QUAD
#define DEVELOPER #define DEVELOPER
#undef CORR
#define C_GHOST 0.0004 #define C_GHOST 0.0004
#define CM_ROWS 20 #define CM_ROWS 20
@ -36,12 +36,12 @@
#include "ansi.h" #include "ansi.h"
#include <iostream> #include <iostream>
//#include <chrono> #include <chrono>
#include <ctime> // time_t #include <ctime> // time_t
#include <cstdio> #include <cstdio>
using namespace std; using namespace std;
//using namespace std::chrono; using namespace std::chrono;
//#define SLS_DETECTOR_JSON_HEADER_VERSION 0x2 //#define SLS_DETECTOR_JSON_HEADER_VERSION 0x2
@ -77,9 +77,10 @@ int main(int argc, char *argv[]) {
int ok; int ok;
// high_resolution_clock::time_point t1; high_resolution_clock::time_point t1;
// high_resolution_clock::time_point t2 ; high_resolution_clock::time_point t2 ;
time_t begin,end,finished; std::chrono::steady_clock::time_point begin,end,finished;
//time_t begin,end,finished;
int rms=0; int rms=0;
@ -141,8 +142,12 @@ int main(int argc, char *argv[]) {
int ncol_cm=CM_ROWS; int ncol_cm=CM_ROWS;
double xt_ghost=C_GHOST; double xt_ghost=C_GHOST;
moench03CommonMode *cm=new moench03CommonMode(ncol_cm); moench03CommonMode *cm=NULL;
moench03GhostSummation *gs=new moench03GhostSummation(det, xt_ghost); moench03GhostSummation *gs=NULL;
#ifdef CORR
cm=new moench03CommonMode(ncol_cm);
gs=new moench03GhostSummation(det, xt_ghost);
#endif
double *gainmap=NULL; double *gainmap=NULL;
float *gm; float *gm;
double *gmap=NULL; double *gmap=NULL;
@ -170,7 +175,7 @@ int main(int argc, char *argv[]) {
//analogDetector<uint16_t> *filter=new analogDetector<uint16_t>(det,1,NULL,1000); //analogDetector<uint16_t> *filter=new analogDetector<uint16_t>(det,1,NULL,1000);
#ifndef INTERP #ifndef INTERP
singlePhotonDetector *filter=new singlePhotonDetector(det,3, nSigma, 1, cm, 1000, 10, -1, -1, gainmap, gs); singlePhotonDetector *filter=new singlePhotonDetector(det,3, nSigma, 1, cm, 1000, 100, -1, -1, gainmap, gs);
multiThreadedCountingDetector *mt=new multiThreadedCountingDetector(filter,nthreads,fifosize); multiThreadedCountingDetector *mt=new multiThreadedCountingDetector(filter,nthreads,fifosize);
@ -356,10 +361,13 @@ int main(int argc, char *argv[]) {
#endif #endif
// if (!zmqsocket->ReceiveHeader(0, acqIndex, frameIndex, subframeIndex, filename, fileindex)) { // if (!zmqsocket->ReceiveHeader(0, acqIndex, frameIndex, subframeIndex, filename, fileindex)) {
// cprintf(RED, "Got Dummy\n"); cprintf(RED, "Got Dummy\n");
// t1=high_resolution_clock::now(); // t1=high_resolution_clock::now();
time(&end); //time(&end);
//cout << "Measurement lasted " << difftime(end,begin) << endl;
end = std::chrono::steady_clock::now();
cout << "Measurement lasted " << (end-begin).count()*0.000001 << " ms" << endl;
while (mt->isBusy()) {;}//wait until all data are processed from the queues while (mt->isBusy()) {;}//wait until all data are processed from the queues
@ -521,14 +529,16 @@ int main(int argc, char *argv[]) {
mt->clearImage(); mt->clearImage();
newFrame=1; newFrame=1;
//t2 = high_resolution_clock::now();
time(&finished);
// auto meas_duration = duration_cast<microseconds>( t2 - t0 ).count();
// auto real_duration = duration_cast<microseconds>( t2 - t1 ).count();
cout << "Measurement lasted " << difftime(end,begin) << endl; //time(&finished);
cout << "Processing lasted " << difftime(finished,begin) << endl; //cout << "Processing lasted " << difftime(finished,begin) << endl;
finished = std::chrono::steady_clock::now();
cout << "Processing lasted " << (finished-begin).count()*0.000001 << " ms" << endl;
#ifdef OPTIMIZE
return 0;
#endif
continue; //continue to not get out continue; //continue to not get out
@ -536,7 +546,8 @@ int main(int argc, char *argv[]) {
#ifdef NEWZMQ #ifdef NEWZMQ
if (newFrame) { if (newFrame) {
time(&begin); begin = std::chrono::steady_clock::now();
//time(&begin);
// t0 = high_resolution_clock::now(); // t0 = high_resolution_clock::now();
//cout <<"new frame" << endl; //cout <<"new frame" << endl;
@ -808,6 +819,7 @@ int main(int argc, char *argv[]) {
// timestamp=doc["timestamp"].GetUint(); // timestamp=doc["timestamp"].GetUint();
packetNumber=doc["packetNumber"].GetUint(); packetNumber=doc["packetNumber"].GetUint();
// cout << acqIndex << " " << frameIndex << " " << subFrameIndex << " "<< bunchId << " " << timestamp << " " << packetNumber << endl; // cout << acqIndex << " " << frameIndex << " " << subFrameIndex << " "<< bunchId << " " << timestamp << " " << packetNumber << endl;
//cprintf(GREEN, "frame\n");
if (packetNumber>=40) { if (packetNumber>=40) {
//*((int*)buff)=frameIndex; //*((int*)buff)=frameIndex;
if (insubframe==0) f0=frameIndex; if (insubframe==0) f0=frameIndex;

2
slsDetectorCalibration/pedestalSubtraction.h
View File

@ -10,7 +10,7 @@ class pedestalSubtraction {
\param nn number of samples to calculate the moving average (defaults to 1000) \param nn number of samples to calculate the moving average (defaults to 1000)
*/ */
pedestalSubtraction (int nn=1000) : stat(nn) {}; pedestalSubtraction (int nn=1000) : stat(nn) {};
/* void setPointers(double *me, double *va) {mean=me; var=va; stat.setPointers(mean, var);}; */
/** virtual destructorr /** virtual destructorr
*/ */
virtual ~pedestalSubtraction() {}; virtual ~pedestalSubtraction() {};

242
slsDetectorCalibration/singlePhotonDetector.h
View File

@ -50,7 +50,6 @@ public analogDetector<uint16_t> {
*/ */
singlePhotonDetector(slsDetectorData<uint16_t> *d, singlePhotonDetector(slsDetectorData<uint16_t> *d,
int csize=3, int csize=3,
@ -58,7 +57,7 @@ public analogDetector<uint16_t> {
int sign=1, int sign=1,
commonModeSubtraction *cm=NULL, commonModeSubtraction *cm=NULL,
int nped=1000, int nped=1000,
int nd=100, int nnx=-1, int nny=-1, double *gm=NULL, ghostSummation<uint16_t> *gs=NULL) : analogDetector<uint16_t>(d, sign, cm, nped, nnx, nny, gm, gs), nDark(nd), eventMask(NULL),nSigma (nsigma), eMin(-1), eMax(-1), clusterSize(csize), clusterSizeY(csize), clusters(NULL), quad(UNDEFINED_QUADRANT), tot(0), quadTot(0) { int nd=100, int nnx=-1, int nny=-1, double *gm=NULL, ghostSummation<uint16_t> *gs=NULL) : analogDetector<uint16_t>(d, sign, cm, nped, nnx, nny, gm, gs), nDark(nd), eventMask(NULL),nSigma (nsigma), eMin(-1), eMax(-1), clusterSize(csize), clusterSizeY(csize), c2(1),c3(1), clusters(NULL), quad(UNDEFINED_QUADRANT), tot(0), quadTot(0) {
@ -66,13 +65,18 @@ public analogDetector<uint16_t> {
fm=new pthread_mutex_t ; fm=new pthread_mutex_t ;
eventMask=new eventType*[ny]; eventMask=new eventType*[ny];
// val=new double*[ny];
for (int i=0; i<ny; i++) { for (int i=0; i<ny; i++) {
eventMask[i]=new eventType[nx]; eventMask[i]=new eventType[nx];
// val[i]=new double[nx];
} }
if (ny==1) if (ny==1)
clusterSizeY=1; clusterSizeY=1;
c2=sqrt((clusterSizeY+1)/2* (clusterSize+1)/2);
c3=sqrt(clusterSizeY*clusterSize);
// cluster=new single_photon_hit(clusterSize,clusterSizeY); // cluster=new single_photon_hit(clusterSize,clusterSizeY);
clusters=new single_photon_hit[nx*ny]; clusters=new single_photon_hit[nx*ny];
@ -100,8 +104,10 @@ public analogDetector<uint16_t> {
myFile=orig->myFile; myFile=orig->myFile;
eventMask=new eventType*[ny]; eventMask=new eventType*[ny];
// val=new double*[ny];
for (int i=0; i<ny; i++) { for (int i=0; i<ny; i++) {
eventMask[i]=new eventType[nx]; eventMask[i]=new eventType[nx];
// val[i]=new double[nx];
} }
eMin=orig->eMin; eMin=orig->eMin;
eMax=orig->eMax; eMax=orig->eMax;
@ -111,6 +117,10 @@ public analogDetector<uint16_t> {
clusterSize=orig->clusterSize; clusterSize=orig->clusterSize;
clusterSizeY=orig->clusterSizeY; clusterSizeY=orig->clusterSizeY;
// cluster=new single_photon_hit(clusterSize,clusterSizeY); // cluster=new single_photon_hit(clusterSize,clusterSizeY);
c2=sqrt((clusterSizeY+1)/2* (clusterSize+1)/2);
c3=sqrt(clusterSizeY*clusterSize);
clusters=new single_photon_hit[nx*ny]; clusters=new single_photon_hit[nx*ny];
// cluster=clusters; // cluster=clusters;
@ -124,7 +134,9 @@ public analogDetector<uint16_t> {
gmap=orig->gmap; gmap=orig->gmap;
nphTot=0; nphTot=0;
nphFrame=0; nphFrame=0;
nphTot=0;
nphFrame=0;
} }
@ -196,7 +208,6 @@ public analogDetector<uint16_t> {
int nn=0; int nn=0;
double max=0, tl=0, tr=0, bl=0,br=0, v; double max=0, tl=0, tr=0, bl=0,br=0, v;
double rms=0; double rms=0;
int cm=0; int cm=0;
if (cmSub) cm=1; if (cmSub) cm=1;
@ -218,8 +229,8 @@ public analogDetector<uint16_t> {
cout << "add to common mode?"<< endl; cout << "add to common mode?"<< endl;
addToCommonMode(data); addToCommonMode(data);
} }
for (int iy=ymin; iy<ymax; iy++) { for (iy=ymin; iy<ymax; ++iy) {
for (int ix=xmin; ix<xmax; ix++) { for (ix=xmin; ix<xmax; ++ix) {
if (det->isGood(ix,iy)) { if (det->isGood(ix,iy)) {
val=subtractPedestal(data,ix,iy, cm); val=subtractPedestal(data,ix,iy, cm);
@ -236,8 +247,8 @@ public analogDetector<uint16_t> {
} }
} }
for (int iy=ymin; iy<ymax; iy++) { for (iy=ymin; iy<ymax; ++iy) {
for (int ix=xmin; ix<xmax; ix++) { for (ix=xmin; ix<xmax; ++ix) {
if (det->isGood(ix,iy)) { if (det->isGood(ix,iy)) {
eventMask[iy][ix]=PEDESTAL; eventMask[iy][ix]=PEDESTAL;
@ -346,18 +357,23 @@ int *getClusters(char *data, int *ph=NULL) {
int nph=0; int nph=0;
double val[ny][nx]; // const int cy=(clusterSizeY+1)/2;
int cy=(clusterSizeY+1)/2; //const int cs=(clusterSize+1)/2;
int cs=(clusterSize+1)/2;
//int ir, ic; //int ir, ic;
eventType ee;
double max=0, tl=0, tr=0, bl=0,br=0, *v; double max=0, tl=0, tr=0, bl=0,br=0, *v;
int cm=0; int cm=0;
int good=1; int good=1;
if (cmSub) cm=1; int ir, ic;
// double quadTot;
//quadrant quad;
double rms;
//if (cmSub) cm=1;
double val[ny][nx];
if (ph==NULL) if (ph==NULL)
ph=image; ph=image;
if (iframe<nDark) { if (iframe<nDark) {
addToPedestal(data); addToPedestal(data);
return 0; return 0;
@ -366,13 +382,15 @@ int *getClusters(char *data, int *ph=NULL) {
if (cm) { if (cmSub) {
addToCommonMode(data); addToCommonMode(data);
cm=1;
} }
for (int iy=ymin; iy<ymax; iy++) { for (iy=ymin; iy<ymax; ++iy) {
for (int ix=xmin; ix<xmax; ix++) { for (ix=xmin; ix<xmax; ++ix) {
if (det->isGood(ix,iy)) { if (det->isGood(ix,iy)==0) continue;
max=0; max=0;
tl=0; tl=0;
tr=0; tr=0;
@ -384,97 +402,120 @@ int *getClusters(char *data, int *ph=NULL) {
eventMask[iy][ix]=PEDESTAL; //eventMask[iy][ix]
ee=PEDESTAL;
(clusters+nph)->rms=getPedestalRMS(ix,iy); rms=getPedestalRMS(ix,iy);
// cluster=clusters+nph; // cluster=clusters+nph;
// cout << ix << " " << iy << endl;
for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) { for (ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) {
for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) { for (ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) {
if ((iy+ir)>=iy && (iy+ir)<ny && (ix+ic)>=ix && (ix+ic)<nx) { if ((iy+ir)>=iy && (iy+ir)<ny && (ix+ic)>=ix && (ix+ic)<nx) {
val[iy+ir][ix+ic]=subtractPedestal(data,ix+ic,iy+ir, cm); val[iy+ir][ix+ic]=subtractPedestal(data,ix+ic,iy+ir, cm);
} v=&(val[iy+ir][ix+ic]);
tot+=*v;
if (ir<=0 && ic<=0)
bl+=*v;
if (ir<=0 && ic>=0)
br+=*v;
if (ir>=0 && ic<=0)
tl+=*v;
if (ir>=0 && ic>=0)
tr+=*v;
if (*v>max) {
max=*v;
}
}
v=&(val[iy+ir][ix+ic]);
tot+=*v;
if (ir<=0 && ic<=0)
bl+=*v;
if (ir<=0 && ic>=0)
br+=*v;
if (ir>=0 && ic<=0)
tl+=*v;
if (ir>=0 && ic>=0)
tr+=*v;
if (*v>max) {
max=*v;
}
if (ir==0 && ic==0) {
if (*v<-nSigma*(clusters+nph)->rms)
eventMask[iy][ix]=NEGATIVE_PEDESTAL;
else if (*v>nSigma*(clusters+nph)->rms)
eventMask[iy][ix]=PHOTON;
}
} }
} }
if (eventMask[iy][ix]==PHOTON && val[iy][ix]<max) /* if (ix==50 && iy==50) */
continue; /* cout << id << " " << ix << " " << iy << " " << det->getValue(data,ix,iy)<< " " << val[iy][ix] << " " << getPedestal(ix,iy) << " " << rms << endl; */
if (val[iy][ix]<-nSigma*rms){
if (bl>=br && bl>=tl && bl>=tr) { ee=NEGATIVE_PEDESTAL;
(clusters+nph)->quad=BOTTOM_LEFT; continue;
(clusters+nph)->quadTot=bl; }
} else if (br>=bl && br>=tl && br>=tr) { if (max>nSigma*rms){
(clusters+nph)->quad=BOTTOM_RIGHT; // cout << "ph1 " << max << " " << nSigma*rms << endl;
(clusters+nph)->quadTot=br; ee=PHOTON;
} else if (tl>=br && tl>=bl && tl>=tr) { if (val[iy][ix]<max)
(clusters+nph)->quad=TOP_LEFT; continue;
(clusters+nph)->quadTot=tl;
} else if (tr>=bl && tr>=tl && tr>=br) {
(clusters+nph)->quad=TOP_RIGHT;
(clusters+nph)->quadTot=tr;
} }
else if (tot>c3*nSigma*rms) {
if (max>nSigma*(clusters+nph)->rms || tot>sqrt(clusterSizeY*clusterSize)*nSigma*(clusters+nph)->rms || ((clusters+nph)->quadTot)>sqrt(cy*cs)*nSigma*(clusters+nph)->rms) { // cout << "ph3 " << tot << " " << c3*nSigma*rms << endl;
if (val[iy][ix]>=max) { ee=PHOTON;
eventMask[iy][ix]=PHOTON_MAX; }
(clusters+nph)->tot=tot; #ifndef WRITE_QUAD
(clusters+nph)->x=ix; else {
(clusters+nph)->y=iy; #endif
// (clusters+nph)->iframe=det->getFrameNumber(data); quad=BOTTOM_RIGHT;
// cout << det->getFrameNumber(data) << " " << (clusters+nph)->iframe << endl; quadTot=br;
(clusters+nph)->ped=getPedestal(ix,iy,0); if (bl>=quadTot) {
for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) { quad=BOTTOM_LEFT;
for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) { quadTot=bl;
(clusters+nph)->set_data(val[iy+ir][ix+ic],ic,ir); }
if (tl>=quadTot) {
quad=TOP_LEFT;
quadTot=tl;
}
if (tr>=quadTot) {
quad=TOP_RIGHT;
quadTot=tr;
}
if (quadTot>c2*nSigma*rms) {
// cout << "ph2 " << quadTot << " " << c2*nSigma*rms << endl;
ee=PHOTON;
} }
}
good=1;
if (eMin>0 && tot<eMin) good=0;
if (eMax>0 && tot>eMax) good=0;
if (good) {
nph++;
image[iy*nx+ix]++;
}
#ifndef WRITE_QUAD
} else {
eventMask[iy][ix]=PHOTON;
}
} else if (eventMask[iy][ix]==PEDESTAL) {
addToPedestal(data,ix,iy,cm);
} }
#endif
if (ee==PHOTON && val[iy][ix]==max) {
ee=PHOTON_MAX;
// cout << "**" <<id<< " " << iframe << " " << nDark << " " << ix << " " << iy << " " << rms << " " << max << " " << quadTot << " " << tot << endl;
(clusters+nph)->tot=tot;
(clusters+nph)->x=ix;
(clusters+nph)->y=iy;
(clusters+nph)->quad=quad;
(clusters+nph)->quadTot=quadTot;
//(clusters+nph)->rms=rms;
// (clusters+nph)->iframe=det->getFrameNumber(data);
// cout << det->getFrameNumber(data) << " " << (clusters+nph)->iframe << endl;
// (clusters+nph)->ped=getPedestal(ix,iy,0);
for (ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) {
for (ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) {
if ((iy+ir)>=iy && (iy+ir)<ny && (ix+ic)>=ix && (ix+ic)<nx)
(clusters+nph)->set_data(val[iy+ir][ix+ic],ic,ir);
}
} }
} good=1;
if (eMin>0 && tot<eMin) good=0;
if (eMax>0 && tot>eMax) good=0;
if (good) {
nph++;
image[iy*nx+ix]++;
}
} else if (ee==PEDESTAL) {
addToPedestal(data,ix,iy,cm);
} /*else {
eventMask[iy][ix]=PHOTON;
}*/
//eventMask[iy][ix]=ee;
} }
}
nphFrame=nph; nphFrame=nph;
nphTot+=nph; nphTot+=nph;
//cout << nphFrame << endl; //cout << nphFrame << endl;
// cout <<"**********************************"<< det->getFrameNumber(data) << " " << nphFrame << endl; //cout <<id << " **********************************"<< iframe << " " << det->getFrameNumber(data) << " " << nphFrame << endl;
writeClusters(det->getFrameNumber(data)); writeClusters(det->getFrameNumber(data));
return image; return image;
@ -561,13 +602,20 @@ int *getClusters(char *data, int *ph=NULL) {
*/ */
static void writeClusters(FILE *f, single_photon_hit *cl, int nph, int fn=0){ static void writeClusters(FILE *f, single_photon_hit *cl, int nph, int fn=0){
if (nph>0) {
/* #ifndef OLDFORMAT */ #ifndef OLDFORMAT
/* if (fwrite((void*)&fn, 1, sizeof(int), f)) */ if (fwrite((void*)&fn, 1, sizeof(int), f))
/* if (fwrite((void*)&nph, 1, sizeof(int), f)) */ if (fwrite((void*)&nph, 1, sizeof(int), f))
/* #endif */ #endif
for (int i=0; i<nph; i++) (cl+i)->write(f); for (int i=0; i<nph; i++) (cl+i)->write(f);
}; }
};
void writeClusters(FILE *f, int fn=0){ void writeClusters(FILE *f, int fn=0){
writeClusters(f,clusters,nphFrame, fn); writeClusters(f,clusters,nphFrame, fn);
//for (int i=0; i<nphFrame; i++) //for (int i=0; i<nphFrame; i++)
@ -623,6 +671,7 @@ int *getClusters(char *data, int *ph=NULL) {
double eMin, eMax; double eMin, eMax;
int clusterSize; /**< cluster size in the x direction */ int clusterSize; /**< cluster size in the x direction */
int clusterSizeY; /**< cluster size in the y direction i.e. 1 for strips, clusterSize for pixels */ int clusterSizeY; /**< cluster size in the y direction i.e. 1 for strips, clusterSize for pixels */
double c2, c3;
// single_photon_hit *cluster; /**< single photon hit data structure */ // single_photon_hit *cluster; /**< single photon hit data structure */
single_photon_hit *clusters; /**< single photon hit data structure */ single_photon_hit *clusters; /**< single photon hit data structure */
quadrant quad; /**< quadrant where the photon is located */ quadrant quad; /**< quadrant where the photon is located */
@ -631,6 +680,7 @@ int *getClusters(char *data, int *ph=NULL) {
int nphTot; int nphTot;
int nphFrame; int nphFrame;
// double **val;
pthread_mutex_t *fm; pthread_mutex_t *fm;
}; };

20
slsDetectorCalibration/single_photon_hit.h
View File

@ -39,13 +39,13 @@ class single_photon_hit {
//fwrite((void*)this, 1, 3*sizeof(int)+4*sizeof(double)+sizeof(quad), myFile); //fwrite((void*)this, 1, 3*sizeof(int)+4*sizeof(double)+sizeof(quad), myFile);
// if (fwrite((void*)this, 1, sizeof(int)+2*sizeof(int16_t), myFile)) // if (fwrite((void*)this, 1, sizeof(int)+2*sizeof(int16_t), myFile))
//#ifdef OLDFORMAT #ifdef OLDFORMAT
if (fwrite((void*)&iframe, 1, sizeof(int), myFile)) {}; if (fwrite((void*)&iframe, 1, sizeof(int), myFile)) {};
//#endif #endif
#ifndef WRITE_QUAD #ifndef WRITE_QUAD
//printf("no quad "); //printf("no quad ");
if (fwrite((void*)&x, 2, sizeof(int16_t), myFile)) //if (fwrite((void*)&x, 2, sizeof(int16_t), myFile))
return fwrite((void*)data, 1, dx*dy*sizeof(int), myFile); return fwrite((void*)&x, 1, dx*dy*sizeof(int)+2*sizeof(int16_t), myFile);
#endif #endif
#ifdef WRITE_QUAD #ifdef WRITE_QUAD
// printf("quad "); // printf("quad ");
@ -91,7 +91,7 @@ class single_photon_hit {
default: default:
; ;
} }
if (fwrite((void*)&x, 2, sizeof(int16_t), myFile)) if (fwrite((void*)&x, 2, sizeof(int16_t), myFile))
return fwrite((void*)qq, 1, 4*sizeof(int), myFile); return fwrite((void*)qq, 1, 4*sizeof(int), myFile);
#endif #endif
return 0; return 0;
@ -104,9 +104,9 @@ class single_photon_hit {
size_t read(FILE *myFile) { size_t read(FILE *myFile) {
//fread((void*)this, 1, 3*sizeof(int)+4*sizeof(double)+sizeof(quad), myFile); //fread((void*)this, 1, 3*sizeof(int)+4*sizeof(double)+sizeof(quad), myFile);
//#ifdef OLDFORMAT #ifdef OLDFORMAT
if (fread((void*)&iframe, 1, sizeof(int), myFile)) {} if (fread((void*)&iframe, 1, sizeof(int), myFile)) {}
//#endif #endif
#ifndef WRITE_QUAD #ifndef WRITE_QUAD
// printf( "no quad \n"); // printf( "no quad \n");
if (fread((void*)&x, 2, sizeof(int16_t), myFile)) if (fread((void*)&x, 2, sizeof(int16_t), myFile))
@ -246,15 +246,15 @@ class single_photon_hit {
int *get_cluster() {return data;}; int *get_cluster() {return data;};
int iframe; /**< frame number */ int iframe; /**< frame number */
int16_t x; /**< x-coordinate of the center of hit */
int16_t y; /**< x-coordinate of the center of hit */
double rms; /**< noise of central pixel l -- at some point it can be removed*/ double rms; /**< noise of central pixel l -- at some point it can be removed*/
double ped; /**< pedestal of the central pixel -- at some point it can be removed*/ double ped; /**< pedestal of the central pixel -- at some point it can be removed*/
double tot; /**< sum of the 3x3 cluster */ double tot; /**< sum of the 3x3 cluster */
quadrant quad; /**< quadrant where the photon is located */ quadrant quad; /**< quadrant where the photon is located */
double quadTot; /**< sum of the maximum 2x2cluster */ double quadTot; /**< sum of the maximum 2x2cluster */
int dx; /**< size of data cluster in x */ int dx; /**< size of data cluster in x */
int dy; /**< size of data cluster in y */ int dy; /**< size of data cluster in y */
int16_t x; /**< x-coordinate of the center of hit */
int16_t y; /**< x-coordinate of the center of hit */
int *data; /**< pointer to data */ int *data; /**< pointer to data */
}; };