some doxy comments added

git-svn-id: file:///afs/psi.ch/project/sls_det_software/svn/slsDetectorCalibration@13 113b152e-814d-439b-b186-022a431db7b5
This commit is contained in:
bergamaschi 2013-12-12 16:05:24 +00:00
parent 1ae2d78bf5
commit 8d09f061da
9 changed files with 599 additions and 258 deletions

View File

@ -6,43 +6,60 @@
class MovingStat class MovingStat
{ {
/** @short approximated moving average structure */
public: public:
/** constructor
\param nn number of samples parameter to be used
*/
MovingStat(int nn=1000) : n(nn), m_n(0) {} MovingStat(int nn=1000) : n(nn), m_n(0) {}
/**
clears the moving average number of samples parameter, mean and standard deviation
*/
void Clear() void Clear()
{ {
m_n = 0; m_n = 0;
m_oldM=0;
m_newM=0; m_newM=0;
m_oldM2=0;
m_newM2=0; m_newM2=0;
} }
void SetN(int i) {n=i;};
/** sets number of samples parameter
\param i number of samples parameter to be set
*/
void SetN(int i) {if (i>=1) n=i;};
/**
gets number of samples parameter
\returns actual number of samples parameter
*/
int GetN() {return n;}; int GetN() {return 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
*/
inline void Calc(double x) { inline void Calc(double x) {
if (m_n<n) Add(x); if (m_n<n) Add(x);
else Push(x); else Push(x);
} }
/** adds the element to the accumulated average and standard deviation
\param x value to add
*/
inline void Add(double x) { inline void Add(double x) {
m_n++; m_n++;
// See Knuth TAOCP vol 2, 3rd edition, page 232
if (m_n == 1) if (m_n == 1)
{ {
m_oldM = m_newM = x; m_newM = x;
m_oldM2 = x*x; m_newM2 = x*x;
} } else {
else m_newM = m_newM + x;
{ m_newM2 = m_newM2 + x*x;
m_newM = m_oldM + x;
m_newM2 = m_oldM2 + x*x;
//m_newM2 = m_oldM2 + (x*x - m_oldM*m_oldM)/m_n;
// set up for next iteration
m_oldM = m_newM;
m_oldM2 = m_newM2;
} }
} }
@ -50,50 +67,64 @@ class MovingStat
inline void Push(double x) inline void Push(double x)
{ {
if (m_n == 1) /** adds the element to the accumulated average and squared mean, while subtracting the current value of the average and squared average
\param x value to push
*/
if (m_n == 0)
{ {
m_oldM = m_newM = x; m_newM = x;
m_oldM2 = x*x; m_newM2 = x*x;
} m_n++;
else } else {
{ m_newM = m_newM + x - m_newM/m_n;
m_newM = m_oldM + (x - m_oldM/m_n); m_newM2 = m_newM2 + x*x - m_newM2/m_n;
m_newM2 = m_oldM2 + (x*x - m_oldM2/m_n);
//m_newM2 = m_oldM2 + (x*x/m_n - m_oldM*m_oldM/(m_n*m_n));
// set up for next iteration
m_oldM = m_newM;
m_oldM2 = m_newM2;
} }
} }
/** returns the current number of elements of the moving average
\returns returns the current number of elements of the moving average
*/
int NumDataValues() const int NumDataValues() const
{ {
return m_n; return m_n;
} }
/** returns the mean, 0 if no elements are inside
\returns returns the mean
*/
inline double Mean() const inline double Mean() const
{ {
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 returns the squared average
*/
double M2() const double M2() const
{ {
return ( (m_n > 1) ? m_newM2/m_n : 0.0 ); return ( (m_n > 1) ? m_newM2/m_n : 0.0 );
} }
/** returns the variance, 0 if no elements are inside
\returns returns the variance
*/
inline double Variance() const inline double Variance() const
{ {
return ( (m_n > 1) ? (M2()-Mean()*Mean()) : 0.0 ); return ( (m_n > 1) ? (M2()-Mean()*Mean()) : 0.0 );
} }
/** returns the standard deviation, 0 if no elements are inside
\returns returns the standard deviation
*/
inline double StandardDeviation() const inline double StandardDeviation() const
{ {
return ( (Variance() > 0) ? sqrt( Variance() ) : -1 ); return ( (Variance() > 0) ? sqrt( Variance() ) : -1 );
} }
private: private:
int n; int n; /**< number of samples parameter */
int m_n; int m_n; /**< current number of elements */
double m_oldM, m_newM, m_oldM2, m_newM2; double m_newM; /**< accumulated average */
double m_newM2; /**< accumulated squared average */
}; };
#endif #endif

View File

@ -5,10 +5,22 @@
class commonModeSubtraction { class commonModeSubtraction {
/** @short class to calculate the common mode of the pedestals based on an approximated moving average*/
public: public:
/** constructor
\param nn number of samples for the moving average to calculate the average common mode
\param iroi number of regions on which one can calculate the common mode separately. Defaults to 1 i.e. whole detector
*/
commonModeSubtraction(int nn=1000, int iroi=1) : cmStat(NULL), cmPed(NULL), nCm(NULL), nROI(iroi) {cmStat=new MovingStat[nROI]; for (int i=0; i<nROI; i++) cmStat[i].SetN(nn); cmPed=new double[nROI]; nCm=new double[nROI];}; commonModeSubtraction(int nn=1000, int iroi=1) : cmStat(NULL), cmPed(NULL), nCm(NULL), nROI(iroi) {cmStat=new MovingStat[nROI]; for (int i=0; i<nROI; i++) cmStat[i].SetN(nn); cmPed=new double[nROI]; nCm=new double[nROI];};
/** destructor - deletes the moving average(s) and the sum of pedestals calculator(s) */
virtual ~commonModeSubtraction() {delete [] cmStat; delete [] cmPed; delete [] nCm;}; virtual ~commonModeSubtraction() {delete [] cmStat; delete [] cmPed; delete [] nCm;};
/** clears the moving average and the sum of pedestals calculation - virtual func*/
virtual void Clear(){ virtual void Clear(){
for (int i=0; i<nROI; i++) { for (int i=0; i<nROI; i++) {
cmStat[i].Clear(); cmStat[i].Clear();
@ -16,6 +28,7 @@ class commonModeSubtraction {
cmPed[i]=0; cmPed[i]=0;
}}; }};
/** adds the average of pedestals to the moving average and reienitialize the calculation of the sum of pedestals for all ROIs. - virtual func*/
virtual void newFrame(){ virtual void newFrame(){
for (int i=0; i<nROI; i++) { for (int i=0; i<nROI; i++) {
if (nCm[i]>0) cmStat[i].Calc(cmPed[i]/nCm[i]); if (nCm[i]>0) cmStat[i].Calc(cmPed[i]/nCm[i]);
@ -23,25 +36,32 @@ class commonModeSubtraction {
cmPed[i]=0; cmPed[i]=0;
}}; }};
virtual void addToCommonMode(double val, int ix=0, int iy=0) { /** adds the pixel to the sum of pedestals -- virtual func
(void)ix; (void)iy; \param isc region of interest index
cmPed[0]+=val; */
nCm[0]++;}; virtual void addToCommonMode(double val, int isc=0) {
if (isc>=0 && isc<nROI) {
cmPed[isc]+=val;
nCm[isc]++;}};
virtual double getCommonMode(int ix=0, int iy=0) { /** gets the common mode i.e. the difference between the current average sum of pedestals mode and the average pedestal -- virtual func
(void)ix; (void)iy; \param isc region of interest index
if (nCm[0]>0) return cmPed[0]/nCm[0]-cmStat[0].Mean(); \return the difference between the current average sum of pedestals and the average pedestal
else return 0;}; */
virtual double getCommonMode(int isc=0) {
if (isc>=0 && isc<nROI)
if (nCm[0]>0) return cmPed[0]/nCm[0]-cmStat[0].Mean();
return 0;};
protected: protected:
MovingStat *cmStat; //stores the common mode average per supercolumn */ MovingStat *cmStat; /**<array of moving average of the pedestal average per region of interest */
double *cmPed; // stores the common mode for this frame per supercolumn */ double *cmPed; /**< array storing the sum of pedestals per region of interest */
double *nCm; // stores the number of pedestals to calculate the cm per supercolumn */ double *nCm; /**< array storing the number of pixels currently contributing to the pedestals */
const int nROI; const int nROI; /**< constant parameter for number of regions on which the common mode should be calculated separately e.g. supercolumns */
}; };

View File

@ -1,39 +1,24 @@
#ifndef MOENCH02MODULEDATA_H #ifndef MOENCH02MODULEDATA_H
#define MOENCH02MODULEDATA_H #define MOENCH02MODULEDATA_H
//#include "slsDetectorData.h" #include "slsReceiverData.h"
#include "singlePhotonDetector.h"
#include "RunningStat.h"
#include "MovingStat.h"
#include <iostream>
#include <fstream>
#include <TTree.h>
using namespace std;
class moench02ModuleData : public slsDetectorData<uint16_t> { class moench02ModuleData : public slsReceiverData<uint16_t> {
public: public:
/** /**
Implements the slsReceiverData structure for the moench02 prototype read out by a module i.e. using the slsReceiver
Constructor (no error checking if datasize and offsets are compatible!) (160x160 pixels, 40 packets 1286 large etc.)
\param npx number of pixels in the x direction \param c crosstalk parameter for the output buffer
\param npy number of pixels in the y direction
\param ds size of the dataset
\param dMap array of size nx*ny storing the pointers to the data in the dataset (as offset)
\param dMask Array of size nx*ny storing the polarity of the data in the dataset (should be 0 if no inversion is required, 0xffffffff is inversion is required)
*/ */
moench02ModuleData(double c=0): slsDetectorData<uint16_t>(160, 160, 40, 1286), xtalk(c) { moench02ModuleData(double c=0): slsReceiverData<uint16_t>(160, 160, 40, 1286), xtalk(c) {
@ -60,7 +45,7 @@ class moench02ModuleData : public slsDetectorData<uint16_t> {
ix=isc*40+ic; ix=isc*40+ic;
iy=ip*16+ir; iy=ip*16+ir;
dMap[iy][ix]=1286*(isc*10+ip)+2*ir*40+2*ic; dMap[iy][ix]=1286*(isc*10+ip)+2*ir*40+2*ic+4;
// cout << ix << " " << iy << " " << dMap[ix][iy] << endl; // cout << ix << " " << iy << " " << dMap[ix][iy] << endl;
} }
} }
@ -85,9 +70,15 @@ class moench02ModuleData : public slsDetectorData<uint16_t> {
}; };
// ~moench02ModuleData() {if (buff) delete [] buff; if (oldbuff) delete [] oldbuff; };
/**
gets the packets number (last packet is labelled with 0 and is replaced with 40)
\param buff pointer to the memory
\returns packet number
*/
int getPacketNumber(char *buff){ int getPacketNumber(char *buff){
int np=(*(int*)buff)&0xff; int np=(*(int*)buff)&0xff;
if (np==0) if (np==0)
@ -95,6 +86,15 @@ class moench02ModuleData : public slsDetectorData<uint16_t> {
return np; return np;
}; };
/**
returns the pixel value as double correcting for the output buffer crosstalk
\param data pointer to the memory
\param ix coordinate in the x direction
\param iy coordinate in the y direction
\returns channel value as double
*/
double getValue(char *data, int ix, int iy=0) { double getValue(char *data, int ix, int iy=0) {
if (xtalk==0 || ix%40==0) if (xtalk==0 || ix%40==0)
return (double)getValue(data, ix, iy); return (double)getValue(data, ix, iy);
@ -102,7 +102,19 @@ class moench02ModuleData : public slsDetectorData<uint16_t> {
return slsDetectorData<uint16_t>::getValue(data, ix, iy)-xtalk*slsDetectorData<uint16_t>::getValue(data, ix-1, iy); return slsDetectorData<uint16_t>::getValue(data, ix, iy)-xtalk*slsDetectorData<uint16_t>::getValue(data, ix-1, iy);
}; };
/** sets the output buffer crosstalk correction parameter
\param c output buffer crosstalk correction parameter to be set
\returns current value for the output buffer crosstalk correction parameter
*/
double setXTalk(double c) {xtalk=c; return xtalk;} double setXTalk(double c) {xtalk=c; return xtalk;}
/** gets the output buffer crosstalk parameter
\returns current value for the output buffer crosstalk correction parameter
*/
double getXTalk() {return xtalk;} double getXTalk() {return xtalk;}
@ -113,7 +125,7 @@ class moench02ModuleData : public slsDetectorData<uint16_t> {
private: private:
double xtalk; double xtalk; /**<output buffer crosstalk correction parameter */
}; };

View File

@ -4,29 +4,32 @@
#include "commonModeSubtraction.h" #include "commonModeSubtraction.h"
class moenchCommonMode : public commonModeSubtraction { class moenchCommonMode : public commonModeSubtraction {
/** @short class to calculate the common mode noise for moench02 i.e. on 4 supercolumns separately */
public: public:
/** constructor - initalizes a commonModeSubtraction with 4 different regions of interest
\param nn number of samples for the moving average
*/
moenchCommonMode(int nn=1000) : commonModeSubtraction(nn,4){} ; moenchCommonMode(int nn=1000) : commonModeSubtraction(nn,4){} ;
/** add value to common mode as a function of the pixel value, subdividing the region of interest in the 4 supercolumns of 40 columns each;
\param val value to add to the common mode
\param ix pixel coordinate in the x direction
\param iy pixel coordinate in the y direction
*/
virtual void addToCommonMode(double val, int ix=0, int iy=0) { virtual void addToCommonMode(double val, int ix=0, int iy=0) {
(void)iy; (void) iy;
int isc=ix/40; commonModeSubtraction::addToCommonMode(val, ix/40);
if (isc>=0 && isc<4) {
cmPed[isc]+=val;
nCm[isc]++;
}
}; };
/**returns common mode value as a function of the pixel value, subdividing the region of interest in the 4 supercolumns of 40 columns each;
\param ix pixel coordinate in the x direction
\param iy pixel coordinate in the y direction
\returns common mode value
*/
virtual double getCommonMode(int ix=0, int iy=0) { virtual double getCommonMode(int ix=0, int iy=0) {
(void)iy; (void) iy;
int isc=ix/40; return commonModeSubtraction::getCommonMode(ix/40);
if (isc>=0 && isc<4) { };
if (nCm[isc]>0) return cmPed[isc]/nCm[isc]-cmStat[isc].Mean();
}
return 0;};
}; };

View File

@ -4,21 +4,45 @@
#include "MovingStat.h" #include "MovingStat.h"
class pedestalSubtraction { class pedestalSubtraction {
/** @short class defining the pedestal subtraction based on an approximated moving average */
public: public:
/** constructor
\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) {};
/** virtual destructorr
*/
virtual ~pedestalSubtraction() {}; virtual ~pedestalSubtraction() {};
/** clears the moving average */
virtual void Clear() {stat.Clear();} virtual void Clear() {stat.Clear();}
/** adds the element to the moving average
\param val value to be added
*/
virtual void addToPedestal(double val){ stat.Calc(val);}; virtual void addToPedestal(double val){ stat.Calc(val);};
/** returns the average value of the pedestal
\returns mean of the moving average
*/
virtual double getPedestal(){return stat.Mean();}; virtual double getPedestal(){return stat.Mean();};
/** returns the standard deviation of the moving average
\returns standard deviation of the moving average
*/
virtual double getPedestalRMS(){return stat.StandardDeviation();}; virtual double getPedestalRMS(){return stat.StandardDeviation();};
/**sets/gets the number of samples for the moving average
\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) {if (i>0) stat.SetN(i); return stat.GetN();};
private: private:
MovingStat stat; MovingStat stat; /**< approximated moving average struct */
}; };
#endif #endif

View File

@ -2,7 +2,6 @@
#define SINGLEPHOTONDETECTOR_H #define SINGLEPHOTONDETECTOR_H
#include <inttypes.h>
#include "slsDetectorData.h" #include "slsDetectorData.h"
#include "single_photon_hit.h" #include "single_photon_hit.h"
@ -37,6 +36,7 @@ using namespace std;
template <class dataType> template <class dataType>
class singlePhotonDetector { class singlePhotonDetector {
/** @short class to perform pedestal subtraction etc. and find single photon clusters for an analog detector */
public: public:
@ -44,11 +44,12 @@ class singlePhotonDetector {
/** /**
Constructor (no error checking if datasize and offsets are compatible!) Constructor (no error checking if datasize and offsets are compatible!)
\param npx number of pixels in the x direction \param s detector data structure to be used
\param npy number of pixels in the y direction \param csize cluster size (should be an odd number). Defaults to 3
\param ds size of the dataset \param nsigma number of rms to discriminate from the noise. Defaults to 5
\param dMap array of size nx*ny storing the pointers to the data in the dataset (as offset) \param cm common mode subtraction algorithm, if any. Defaults to NULL i.e. none
\param dMask Array of size nx*ny storing the polarity of the data in the dataset (should be 0 if no inversion is required, 0xffffffff is inversion is required) \param nped number of samples for pedestal averaging
\param nd number of dark frames to average as pedestals without photon discrimination at the beginning of the measurement
*/ */
@ -63,7 +64,6 @@ class singlePhotonDetector {
int nd=100) : det(d), nx(0), ny(0), stat(NULL), cmSub(cm), nDark(nd), eventMask(NULL),nSigma (nsigma), clusterSize(csize), clusterSizeY(csize), cluster(NULL), iframe(-1), dataSign(sign) { int nd=100) : det(d), nx(0), ny(0), stat(NULL), cmSub(cm), nDark(nd), eventMask(NULL),nSigma (nsigma), clusterSize(csize), clusterSizeY(csize), cluster(NULL), iframe(-1), dataSign(sign) {
det->getDetectorSize(nx,ny); det->getDetectorSize(nx,ny);
@ -80,39 +80,80 @@ class singlePhotonDetector {
clusterSizeY=1; clusterSizeY=1;
cluster=new single_photon_hit(clusterSize,clusterSizeY); cluster=new single_photon_hit(clusterSize,clusterSizeY);
setClusterSize(csize);
}; };
/**
destructor. Deletes the cluster structure and the pdestalSubtraction array
*/
virtual ~singlePhotonDetector() {delete cluster; for (int i=0; i<ny; i++) delete [] stat[i]; delete [] stat;}; virtual ~singlePhotonDetector() {delete cluster; for (int i=0; i<ny; i++) delete [] stat[i]; delete [] stat;};
/** 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(); }; 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(); };
/** resets the eventMask to undefined and the commonModeSubtraction */
void newFrame(){iframe++; for (int iy=0; iy<ny; iy++) for (int ix=0; ix<nx; ix++) eventMask[iy][ix]=UNDEFINED; if (cmSub) cmSub->newFrame();}; void newFrame(){iframe++; for (int iy=0; iy<ny; iy++) for (int ix=0; ix<nx; ix++) eventMask[iy][ix]=UNDEFINED; if (cmSub) cmSub->newFrame();};
/** sets the commonModeSubtraction algorithm to be used
\param cm commonModeSubtraction algorithm to be used (NULL unsets)
\returns pointer to the actual common mode subtraction algorithm
*/
commonModeSubtraction setCommonModeSubtraction(commonModeSubtraction *cm) {cmSub=cm; return cmSub;}; commonModeSubtraction setCommonModeSubtraction(commonModeSubtraction *cm) {cmSub=cm; return cmSub;};
/**
sets the sign of the data
\param sign 1 means positive values for photons, -1 negative, 0 gets
\returns current sign for the data
*/
int setDataSign(int sign=0) {if (sign==1 || sign==-1) dataSign=sign; return dataSign;}; int setDataSign(int sign=0) {if (sign==1 || sign==-1) dataSign=sign; return dataSign;};
/**
adds value to pedestal (and common mode) for the given pixel
\param val value to be added
\param ix pixel x coordinate
\param iy pixel y coordinate
*/
virtual void addToPedestal(double val, int ix, int iy){ if (ix>=0 && ix<nx && iy>=0 && iy<ny) { stat[iy][ix].addToPedestal(val); if (cmSub && det->isGood(ix, iy) ) cmSub->addToCommonMode(val, ix, iy);}}; virtual void addToPedestal(double val, int ix, int iy){ if (ix>=0 && ix<nx && iy>=0 && iy<ny) { stat[iy][ix].addToPedestal(val); if (cmSub && det->isGood(ix, iy) ) cmSub->addToCommonMode(val, ix, iy);}};
/**
gets pedestal (and common mode)
\param ix pixel x coordinate
\param iy pixel y coordinate
\param cm 0 (default) without common mode subtraction, 1 with common mode subtraction (if defined)
*/
virtual double getPedestal(int ix, int iy, int cm=0){if (ix>=0 && ix<nx && iy>=0 && iy<ny) if (cmSub && cm>0) return stat[iy][ix].getPedestal()-cmSub->getCommonMode(); else return stat[iy][ix].getPedestal(); else return -1;}; virtual double getPedestal(int ix, int iy, int cm=0){if (ix>=0 && ix<nx && iy>=0 && iy<ny) if (cmSub && cm>0) return stat[iy][ix].getPedestal()-cmSub->getCommonMode(); else return stat[iy][ix].getPedestal(); else return -1;};
/**
gets pedestal rms (i.e. noise)
\param ix pixel x coordinate
\param iy pixel y coordinate
*/
double getPedestalRMS(int ix, int iy){if (ix>=0 && ix<nx && iy>=0 && iy<ny) return stat[iy][ix].getPedestalRMS();else return -1;}; double getPedestalRMS(int ix, int iy){if (ix>=0 && ix<nx && iy>=0 && iy<ny) return stat[iy][ix].getPedestalRMS();else return -1;};
/** 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
*/
double setNSigma(double n=-1){if (n>0) nSigma=n; return nSigma;} double setNSigma(double n=-1){if (n>0) nSigma=n; return nSigma;}
/** sets/gets cluster size
\param n cluster size to be set, (0 or negative gets). If even is incremented by 1.
\returns actual cluster size
*/
int setClusterSize(int n=-1){ int setClusterSize(int n=-1){
if (n>0 && n!=clusterSize) { if (n>0 && n!=clusterSize) {
if (n%2==0) if (n%2==0)
n+=1; n+=1;
clusterSize=n; clusterSize=n;
delete cluster; if (cluster)
delete cluster;
if (ny>1) if (ny>1)
clusterSizeY=clusterSize; clusterSizeY=clusterSize;
cluster=new single_photon_hit(clusterSize,clusterSizeY); cluster=new single_photon_hit(clusterSize,clusterSizeY);
} }
return clusterSize; return clusterSize;
@ -120,6 +161,17 @@ class singlePhotonDetector {
/** finds event type for pixel and fills cluster structure. The algorithm loops only if the evenMask for this pixel is still undefined.
if pixel or cluster around it are above threshold (nsigma*pedestalRMS) cluster is filled and pixel mask is PHOTON_MAX (if maximum in cluster) or NEIGHBOUR; If PHOTON_MAX, the elements of the cluster are also set as NEIGHBOURs in order to speed up the looping
if below threshold the pixel is either marked as PEDESTAL (and added to the pedestal calculator) or NEGATIVE_PEDESTAL is case it's lower than -threshold, otherwise the pedestal average would drift to negative values while it should be 0.
/param data pointer to the data
/param ix pixel x coordinate
/param iy pixel y coordinate
/param cm enable(1)/disable(0) common mode subtraction (if defined).
/returns event type for the given pixel
*/
eventType getEventType(char *data, int ix, int iy, int cm=0) { eventType getEventType(char *data, int ix, int iy, int cm=0) {
// eventType ret=PEDESTAL; // eventType ret=PEDESTAL;
@ -194,13 +246,32 @@ class singlePhotonDetector {
/** sets/gets number of samples for moving average pedestal calculation
\param i number of samples to be set (0 or negative gets)
\returns actual number of samples
*/
int SetNPedestals(int i=-1) {int ix=0, iy=0; if (i>0) for (ix=0; ix<nx; ix++) for (iy=0; iy<ny; iy++) stat[iy][ix].SetNPedestals(i); return stat[0][0].SetNPedestals();}; int SetNPedestals(int i=-1) {int ix=0, iy=0; if (i>0) for (ix=0; ix<nx; ix++) for (iy=0; iy<ny; iy++) stat[iy][ix].SetNPedestals(i); return stat[0][0].SetNPedestals();};
double getClusterElement(int ic, int ir){return cluster->get_data(ic,ir);}; /** returns value for cluster element in relative coordinates
eventType getEventMask(int ic, int ir){return eventMask[ir][ic];}; \param ic x coordinate (center is (0,0))
\param ir y coordinate (center is (0,0))
\returns cluster element
*/
double getClusterElement(int ic, int ir=0){return cluster->get_data(ic,ir);};
/** returns event mask for the given pixel
\param ic x coordinate (center is (0,0))
\param ir y coordinate (center is (0,0))
\returns event mask enum for the given pixel
*/
eventType getEventMask(int ic, int ir=0){return eventMask[ir][ic];};
#ifdef MYROOT1 #ifdef MYROOT1
/** generates a tree and maps the branches
\param tname name for the tree
\param iframe pointer to the frame number
\returns returns pointer to the TTree
*/
TTree *initEventTree(char *tname, int *iFrame=NULL) { TTree *initEventTree(char *tname, int *iFrame=NULL) {
TTree* tall=new TTree(tname,tname); TTree* tall=new TTree(tname,tname);
@ -223,22 +294,22 @@ class singlePhotonDetector {
private: private:
slsDetectorData<dataType> *det; slsDetectorData<dataType> *det; /**< slsDetectorData to be used */
int nx, ny; int nx; /**< Size of the detector in x direction */
int ny; /**< Size of the detector in y direction */
pedestalSubtraction **stat; pedestalSubtraction **stat; /**< pedestalSubtraction class */
commonModeSubtraction *cmSub; commonModeSubtraction *cmSub;/**< commonModeSubtraction class */
//MovingStat **stat; int nDark; /**< number of frames to be used at the beginning of the dataset to calculate pedestal without applying photon discrimination */
int nDark; eventType **eventMask; /**< matrix of event type or each pixel */
eventType **eventMask; double nSigma; /**< number of sigma parameter for photon discrimination */
double nSigma; int clusterSize; /**< cluster size in the x direction */
int clusterSize, clusterSizeY; int clusterSizeY; /**< cluster size in the y direction i.e. 1 for strips, clusterSize for pixels */
single_photon_hit *cluster; single_photon_hit *cluster; /**< single photon hit data structure */
int iframe; int iframe; /**< frame number (not from file but incremented within the dataset every time newFrame is called */
int dataSign; int dataSign; /**< sign of the data i.e. 1 if photon is positive, -1 if negative */
/* int cx, cy; */
}; };

View File

@ -5,38 +5,56 @@ typedef double double32_t;
typedef float float32_t; typedef float float32_t;
typedef int int32_t; typedef int int32_t;
/* /\** */
/* @short structure for a single photon hit */
/* *\/ */
/* typedef struct{ */
/* double* data; /\**< data size *\/ */
/* int x; /\**< x-coordinate of the center of hit *\/ */
/* int y; /\**< x-coordinate of the center of hit *\/ */
/* double rms; /\**< noise of central pixel *\/ */
/* double ped; /\**< pedestal of the central pixel *\/ */
/* int iframe; /\**< frame number *\/ */
/* }single_photon_hit; */
class single_photon_hit { class single_photon_hit {
/** @short Structure for a single photon hit */
public: public:
/** constructor, instantiates the data array -- all class elements are public!
\param nx cluster size in x direction
\param ny cluster size in y direction (defaults to 1 for 1D detectors)
*/
single_photon_hit(int nx, int ny=1): dx(nx), dy(ny) {data=new double[dx*dy];}; single_photon_hit(int nx, int ny=1): dx(nx), dy(ny) {data=new double[dx*dy];};
~single_photon_hit(){delete [] data;};
~single_photon_hit(){delete [] data;}; /**< destructor, deletes the data array */
/** binary write to file of all elements of the structure, except size of the cluster
\param myFile file descriptor
*/
void write(FILE *myFile) {fwrite((void*)this, 1, 3*sizeof(int)+2*sizeof(double), myFile); fwrite((void*)data, 1, dx*dy*sizeof(double), myFile);}; void write(FILE *myFile) {fwrite((void*)this, 1, 3*sizeof(int)+2*sizeof(double), myFile); fwrite((void*)data, 1, dx*dy*sizeof(double), myFile);};
/**
binary read from file of all elements of the structure, except size of the cluster. The structure is then filled with those args
\param myFile file descriptor
*/
void read(FILE *myFile) {fread((void*)this, 1, 3*sizeof(int)+2*sizeof(double), myFile); fread((void*)data, 1, dx*dy*sizeof(double), myFile);}; void read(FILE *myFile) {fread((void*)this, 1, 3*sizeof(int)+2*sizeof(double), myFile); fread((void*)data, 1, dx*dy*sizeof(double), myFile);};
/**
assign the value to the element of the cluster matrix, with relative coordinates where the center of the cluster is (0,0)
\param v value to be set
\param ix coordinate x within the cluster (center is (0,0))
\param iy coordinate y within the cluster (center is (0,0))
*/
void set_data(double v, int ix, int iy=0){data[(iy+dy/2)*dx+ix+dx/2]=v;}; void set_data(double v, int ix, int iy=0){data[(iy+dy/2)*dx+ix+dx/2]=v;};
/**
gets the value to the element of the cluster matrix, with relative coordinates where the center of the cluster is (0,0)
\param ix coordinate x within the cluster (center is (0,0))
\param iy coordinate y within the cluster (center is (0,0))
\returns value of the cluster element
*/
double get_data(int ix, int iy=0){return data[(iy+dy/2)*dx+ix+dx/2];}; double get_data(int ix, int iy=0){return data[(iy+dy/2)*dx+ix+dx/2];};
int x; /**< x-coordinate of the center of hit */
int x; /**< x-coordinate of the center of hit */ int y; /**< x-coordinate of the center of hit */
int 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 */ double ped; /**< pedestal of the central pixel -- at some point it can be removed*/
double ped; /**< pedestal of the central pixel */ int iframe; /**< frame number */
int iframe; /**< frame number */ double *data; /**< pointer to data */
double *data; /**< data size */ const int dx; /**< size of data cluster in x */
const int dx; const int dy; /**< size of data cluster in y */
const int dy;
}; };

View File

@ -1,6 +1,7 @@
#ifndef SLSDETECTORDATA_H #ifndef SLSDETECTORDATA_H
#define SLSDETECTORDATA_H #define SLSDETECTORDATA_H
#include <inttypes.h>
#include <iostream> #include <iostream>
#include <fstream> #include <fstream>
@ -15,16 +16,19 @@ class slsDetectorData {
/** /**
Constructor (no error checking if datasize and offsets are compatible!)
\param npx number of pixels in the x direction
\param npy number of pixels in the y direction
\param ds size of the dataset
\param dMap array of size nx*ny storing the pointers to the data in the dataset (as offset)
\param dMask Array of size nx*ny storing the polarity of the data in the dataset (should be 0 if no inversion is required, 0xffffffff is inversion is required)
General slsDetectors data structure. Works for data acquired using the slsDetectorReceiver. Can be generalized to other detectors (many virtual funcs).
Constructor (no error checking if datasize and offsets are compatible!)
\param npx number of pixels in the x direction
\param npy number of pixels in the y direction (1 for strips)
\param dsize size of the data
\param dMap array of size nx*ny storing the pointers to the data in the dataset (as offset)
\param dMask Array of size nx*ny storing the polarity of the data in the dataset (should be 0 if no inversion is required, 0xffffffff is inversion is required)
\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 np, int psize, int **dMap=NULL, dataType **dMask=NULL, int **dROI=NULL): nx(npx), ny(npy), nPackets(np), packetSize(psize) { slsDetectorData(int npx, int npy, int dsize, int **dMap=NULL, dataType **dMask=NULL, int **dROI=NULL): nx(npx), ny(npy), dataSize(dsize) {
@ -59,13 +63,20 @@ class slsDetectorData {
delete [] dataROIMask; delete [] dataROIMask;
} }
/**
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);)
*/
void setDataMap(int **dMap=NULL) { void setDataMap(int **dMap=NULL) {
if (dMap==NULL) { if (dMap==NULL) {
for (int iy=0; iy<ny; iy++) for (int iy=0; iy<ny; iy++)
for (int ix=0; ix<nx; ix++) for (int ix=0; ix<nx; ix++)
dataMap[iy][ix]=ix*ny+ix; dataMap[iy][ix]=(iy*nx+ix)*sizeof(dataType);
} else { } else {
for (int iy=0; iy<ny; iy++) for (int iy=0; iy<ny; iy++)
for (int ix=0; ix<nx; ix++) for (int ix=0; ix<nx; ix++)
@ -74,6 +85,14 @@ class slsDetectorData {
}; };
/**
defines the data mask i.e. the polarity of the data
\param dMask Array of size nx*ny storing the polarity of the data in the dataset (should be 0 if no inversion is required, 0xffffffff is inversion is required)
*/
void setDataMask(dataType **dMask=NULL){ void setDataMask(dataType **dMask=NULL){
if (dMask!=NULL) { if (dMask!=NULL) {
@ -81,9 +100,18 @@ class slsDetectorData {
for (int iy=0; iy<ny; iy++) for (int iy=0; iy<ny; iy++)
for (int ix=0; ix<nx; ix++) for (int ix=0; ix<nx; ix++)
dataMask[iy][ix]=dMask[iy][ix]; dataMask[iy][ix]=dMask[iy][ix];
} else {
for (int iy=0; iy<ny; iy++)
for (int ix=0; ix<nx; ix++)
dataMask[iy][ix]=0;
} }
}; };
/**
defines the region of interest and/or the bad channels mask
\param dROI Array of size nx*ny. The lements are 1s if the channel is good or in the ROI, 0 is bad or out of the ROI. NULL (default) means all 1s.
*/
void setDataROIMask(int **dROI=NULL){ void setDataROIMask(int **dROI=NULL){
if (dROI!=NULL) { if (dROI!=NULL) {
@ -102,21 +130,37 @@ class slsDetectorData {
}; };
/**
Define bad channel or roi mask for a single channel
\param ix channel x coordinate
\param iy channel y coordinate (1 for strips)
\param i 1 if pixel is good (or in the roi), 0 if bad
\returns 1 if pixel is good, 0 if it's bad, -1 if pixel is out of range
*/
int setGood(int ix, int iy, int i=1) { if (ix>=0 && ix<nx && iy>=0 && iy<ny) dataROIMask[iy][ix]=i; return isGood(ix,iy);}; int setGood(int ix, int iy, int i=1) { if (ix>=0 && ix<nx && iy>=0 && iy<ny) dataROIMask[iy][ix]=i; return isGood(ix,iy);};
/**
Define bad channel or roi mask for a single channel
\param ix channel x coordinate
\param iy channel y coordinate (1 for strips)
\returns 1 if pixel is good, 0 if it's bad, -1 if pixel is out of range
*/
int isGood(int ix, int iy) { if (ix>=0 && ix<nx && iy>=0 && iy<ny) return dataROIMask[iy][ix]; else return -1;}; int isGood(int ix, int iy) { if (ix>=0 && ix<nx && iy>=0 && iy<ny) return dataROIMask[iy][ix]; else return -1;};
/**
void getDetectorSize(int &npx, int &npy){npx=nx; npy=ny;}; Returns detector size in x,y
\param nx reference to number of channels in x
\param ny reference to number of channels in y (will be 1 for strips)
\returns total number of channels
*/
int getDetectorSize(int &npx, int &npy){npx=nx; npy=ny; return nx*ny;};
/** /**
Returns the value of the selected channel for the given dataset (no error checking if number of pixels are compatible!) Returns the value of the selected channel for the given dataset. Virtual function, can be overloaded.
\param data pointer to the dataset (including headers etc) \param data pointer to the dataset (including headers etc)
\param ix pixel numb er in the x direction \param ix pixel number in the x direction
\param iy pixel number in the y direction \param iy pixel number in the y direction
\returns data for the selected channel, with inversion if required \returns data for the selected channel, with inversion if required
*/ */
@ -124,119 +168,72 @@ class slsDetectorData {
virtual dataType getChannel(char *data, int ix, int iy=0) { virtual dataType getChannel(char *data, int ix, int iy=0) {
dataType m=0, d=0; dataType m=0, d=0;
if (ix>=0 && ix<nx && iy>=0 && iy<ny && dataMap[iy][ix]>=0 && dataMap[iy][ix]<nPackets*packetSize) { if (ix>=0 && ix<nx && iy>=0 && iy<ny && dataMap[iy][ix]>=0 && dataMap[iy][ix]<dataSize) {
m=dataMask[iy][ix]; m=dataMask[iy][ix];
d=*((dataType*)(data+dataMap[iy][ix])); d=*((dataType*)(data+dataMap[iy][ix]));
} }
return d^m; return d^m;
}; };
/**
Returns the value of the selected channel for the given dataset as double.
\param data pointer to the dataset (including headers etc)
\param ix pixel number in the x direction
\param iy pixel number in the y direction
\returns data for the selected channel, with inversion if required as double
*/
virtual double getValue(char *data, int ix, int iy=0) {return (double)getChannel(data, ix, iy);}; virtual double getValue(char *data, int ix, int iy=0) {return (double)getChannel(data, ix, iy);};
/**
virtual int getFrameNumber(char *buff){return ((*(int*)buff)&(0xffffff00))>>8;}; Returns the frame number for the given dataset. Purely virtual func.
virtual int getPacketNumber(char *buff) {return (*(int*)buff)&0xff;}; \param buff pointer to the dataset
\returns frame number
*/
virtual char *findNextFrame(char *data, int &npackets, int dsize) { virtual int getFrameNumber(char *buff)=0;
char *retval=NULL, *p=data;
int dd=0;
int fn, fnum=-1, np=0, pnum=-1;
while (dd<=(dsize-packetSize)) {
pnum=getPacketNumber(p);
fn=getFrameNumber(p);
if (pnum<0 || pnum>=nPackets) {
cout << "Bad packet number " << pnum << " frame "<< fn << endl;
retval=NULL;
continue;
}
if (pnum==1) {
fnum=fn;
retval=p;
if (np>0)
cout << "*Incomplete frame number " << fnum << endl;
np=0;
} else if (fn!=fnum) {
if (fnum!=-1) {
cout << " **Incomplete frame number " << fnum << " pnum " << pnum << " " << getFrameNumber(p) << endl;
retval=NULL;
}
np=0;
}
p+=packetSize;
dd+=packetSize;
np++;
if (np==nPackets)
if (pnum==nPackets)
break;
else {
cout << "Too many packets for this frame! "<< fnum << " " << pnum << endl;
retval=NULL;
}
}
if (np<40) {
if (np>0)
cout << "Too few packets for this frame! "<< fnum << " " << pnum << endl;
}
npackets=np; /**
return retval; Returns the packet number for the given dataset. purely virtual func
}; \param buff pointer to the dataset
\returns packet number number
*/
virtual int getPacketNumber(char *buff)=0;
virtual char *readNextFrame(ifstream &filebin) { /**
Loops over a memory slot until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). purely virtual func
\param data pointer to the memory to be analyzed
\param ndata reference to the amount of data found for the frame, in case the frame is incomplete at the end of the memory slot
\param dsize size of the memory slot to be analyzed
\returns pointer to the beginning of the last good frame (might be incomplete if ndata smaller than dataSize), or NULL if no frame is found
*/
virtual char *findNextFrame(char *data, int &ndata, int dsize)=0;
/* if (oldbuff) */ /**
/* delete [] oldbuff; */
/* oldbuff=buff; */ Loops over a file stream until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). Can be overloaded for different kind of detectors!
\param filebin input file stream (binary)
\returns pointer to the begin of the last good frame, NULL if no frame is found or last frame is incomplete
// char *p=new char[1286]; */
char *data=new char[packetSize*nPackets]; virtual char *readNextFrame(ifstream &filebin)=0;
char *retval=0;
int np=40;
if (filebin.is_open()) {
while (filebin.read(data+np*packetSize,packetSize)) {
if (np==nPackets) {
retval=findNextFrame(data,np,packetSize*nPackets);
if (retval==data && np==nPackets)
return data;
else if (np>nPackets) {
cout << "too many packets!!!!!!!!!!" << endl;
delete [] data;
return NULL;
} else {
for (int ip=0; ip<np; ip++)
memcpy(data+ip*packetSize,retval+ip*packetSize,packetSize);
}
} else if (np>nPackets) {
cout << "*******too many packets!!!!!!!!!!" << endl;
delete [] data;
return NULL;
} else {
// memcpy(data+np*1286,p,1286);
np++;
}
}
}
delete [] data;
return NULL;
};
protected:
private:
const int nx; /**< Number of pixels in the x direction */ const int nx; /**< Number of pixels in the x direction */
const int ny; /**< Number of pixels in the y direction */ const int ny; /**< Number of pixels in the y direction */
const int nPackets; /**<number of UDP packets constituting one frame */ const int dataSize; /**<size of the data constituting one frame */
const int packetSize; /**< size of a udp packet */
int **dataMap; /**< Array of size nx*ny storing the pointers to the data in the dataset (as offset)*/ int **dataMap; /**< Array of size nx*ny storing the pointers to the data in the dataset (as offset)*/
dataType **dataMask; /**< Array of size nx*ny storing the polarity of the data in the dataset (should be 0 if no inversion is required, 0xffffffff is inversion is required) */ dataType **dataMask; /**< Array of size nx*ny storing the polarity of the data in the dataset (should be 0 if no inversion is required, 0xffffffff is inversion is required) */
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) */

View File

@ -0,0 +1,165 @@
#ifndef SLSRECEIVERDATA_H
#define SLSRECEIVERDATA_H
#include "slsDetectorData.h"
template <class dataType>
class slsReceiverData : public slsDetectorData<dataType> {
public:
/**
slsReceiver data structure. Works for data acquired using the slsDetectorReceiver subdivided in different packets with headers and footers.
Inherits and implements slsDetectorData.
Constructor (no error checking if datasize and offsets are compatible!)
\param npx number of pixels in the x direction
\param npy number of pixels in the y direction (1 for strips)
\param np number of packets
\param psize packets size
\param dMap array of size nx*ny storing the pointers to the data in the dataset (as offset)
\param dMask Array of size nx*ny storing the polarity of the data in the dataset (should be 0 if no inversion is required, 0xffffffff is inversion is required)
\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.
*/
slsReceiverData(int npx, int npy, int np, int psize, int **dMap=NULL, dataType **dMask=NULL, int **dROI=NULL): slsDetectorData<dataType>(npx, npy, np*psize, dMap, dMask, dROI), nPackets(np), packetSize(psize) {};
/**
Returns the frame number for the given dataset. Virtual func: works for slsDetectorReceiver data (also for each packet), but can be overloaded.
\param buff pointer to the dataset
\returns frame number
*/
virtual int getFrameNumber(char *buff){return ((*(int*)buff)&(0xffffff00))>>8;};
/**
Returns the packet number for the given dataset. Virtual func: works for slsDetectorReceiver packets, but can be overloaded.
\param buff pointer to the dataset
\returns packet number number
*/
virtual int getPacketNumber(char *buff) {return (*(int*)buff)&0xff;};
/**
Loops over a memory slot until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). Can be overloaded for different kind of detectors!
\param data pointer to the memory to be analyzed
\param ndata
\param dsize size of the memory slot to be analyzed
\returns pointer to the first packet of the last good frame (might be incomplete if npackets lower than the number of packets), or NULL if no frame is found
*/
virtual char *findNextFrame(char *data, int &npackets, int dsize) {
char *retval=NULL, *p=data;
int dd=0;
int fn, fnum=-1, np=0, pnum=-1;
while (dd<=(dsize-packetSize)) {
pnum=getPacketNumber(p);
fn=getFrameNumber(p);
if (pnum<0 || pnum>=nPackets) {
cout << "Bad packet number " << pnum << " frame "<< fn << endl;
retval=NULL;
continue;
}
if (pnum==1) {
fnum=fn;
retval=p;
if (np>0)
cout << "*Incomplete frame number " << fnum << endl;
np=0;
} else if (fn!=fnum) {
if (fnum!=-1) {
cout << " **Incomplete frame number " << fnum << " pnum " << pnum << " " << getFrameNumber(p) << endl;
retval=NULL;
}
np=0;
}
p+=packetSize;
dd+=packetSize;
np++;
if (np==nPackets)
if (pnum==nPackets)
break;
else {
cout << "Too many packets for this frame! "<< fnum << " " << pnum << endl;
retval=NULL;
}
}
if (np<40) {
if (np>0)
cout << "Too few packets for this frame! "<< fnum << " " << pnum << endl;
}
npackets=np*packetSize;
return retval;
};
/**
Loops over a file stream until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). Can be overloaded for different kind of detectors!
\param filebin input file stream (binary)
\returns pointer to the first packet of the last good frame, NULL if no frame is found or last frame is incomplete
*/
virtual char *readNextFrame(ifstream &filebin) {
char *data=new char[packetSize*nPackets];
char *retval=0;
int np=0, nd;
if (filebin.is_open()) {
while (filebin.read(data+np*packetSize,packetSize)) {
if (np==(nPackets-1)) {
retval=findNextFrame(data,nd,packetSize*nPackets);
np=nd/packetSize;
if (retval==data && np==nPackets)
return data;
else if (np>nPackets) {
cout << "too many packets!!!!!!!!!!" << endl;
delete [] data;
return NULL;
} else {
for (int ip=0; ip<np; ip++)
memcpy(data+ip*packetSize,retval+ip*packetSize,packetSize);
}
} else if (np>nPackets) {
cout << "*******too many packets!!!!!!!!!!" << endl;
delete [] data;
return NULL;
} else {
np++;
}
}
}
delete [] data;
return NULL;
};
private:
const int nPackets; /**<number of UDP packets constituting one frame */
const int packetSize; /**< size of a udp packet */
};
#endif