From 01f87cff965efd408d160acda0019b0ffa4ab2c3 Mon Sep 17 00:00:00 2001 From: bergamaschi Date: Mon, 9 Dec 2013 10:26:50 +0000 Subject: [PATCH] added general class to detect photons git-svn-id: file:///afs/psi.ch/project/sls_det_software/svn/slsDetectorCalibration@11 113b152e-814d-439b-b186-022a431db7b5 --- slsDetectorCalibration/MovingStat.h | 11 +- slsDetectorCalibration/moench02ModuleData.h | 90 +------ slsDetectorCalibration/moenchReadData.C | 6 +- slsDetectorCalibration/singlePhotonDetector.h | 225 ++++++++++++++++++ slsDetectorCalibration/single_photon_hit.h | 44 ++++ slsDetectorCalibration/slsDetectorData.h | 144 ++++++++--- 6 files changed, 409 insertions(+), 111 deletions(-) create mode 100644 slsDetectorCalibration/singlePhotonDetector.h create mode 100644 slsDetectorCalibration/single_photon_hit.h diff --git a/slsDetectorCalibration/MovingStat.h b/slsDetectorCalibration/MovingStat.h index 634560000..9c8584d33 100755 --- a/slsDetectorCalibration/MovingStat.h +++ b/slsDetectorCalibration/MovingStat.h @@ -1,7 +1,13 @@ - class MovingStat +#ifndef MOVINGSTAT_H +#define MOVINGSTAT_H + +#include + + +class MovingStat { public: - MovingStat() : m_n(0), n(1000) {} + MovingStat() : n(1000), m_n(0) {} void Clear() { @@ -86,3 +92,4 @@ int m_n; double m_oldM, m_newM, m_oldM2, m_newM2; }; +#endif diff --git a/slsDetectorCalibration/moench02ModuleData.h b/slsDetectorCalibration/moench02ModuleData.h index d7757972b..8a9a2d579 100644 --- a/slsDetectorCalibration/moench02ModuleData.h +++ b/slsDetectorCalibration/moench02ModuleData.h @@ -1,6 +1,7 @@ #ifndef MOENCH02MODULEDATA_H #define MOENCH02MODULEDATA_H -#include "slsDetectorData.h" +//#include "slsDetectorData.h" +#include "singlePhotonDetector.h" #include "RunningStat.h" #include "MovingStat.h" @@ -13,19 +14,10 @@ using namespace std; -class moench02ModuleData : public slsDetectorData { +class moench02ModuleData : public slsDetectorData { public: - enum eventType { - PEDESTAL, - NEIGHBOUR, - PHOTON, - PHOTON_MAX, - NEGATIVE_PEDESTAL, - }; - - /** @@ -41,7 +33,7 @@ class moench02ModuleData : public slsDetectorData { */ - moench02ModuleData(int nbg=1000): slsDetectorData(160, 160, 1286*40), nSigma(5), buff(NULL), oldbuff(NULL), pedSub(1) { + moench02ModuleData(int nbg=1000): slsDetectorData(160, 160, 40, 1286), nSigma(5), buff(NULL), oldbuff(NULL), pedSub(1) { int headerOff=0; int footerOff=1284; int dataOff=4; @@ -77,7 +69,7 @@ class moench02ModuleData : public slsDetectorData { ix=isc*40+ic; iy=ip*16+ir; - dMap[iy][ix]=1280*(isc*10+ip)+2*ir*40+2*ic; + dMap[iy][ix]=1286*(isc*10+ip)+2*ir*40+2*ic; // cout << ix << " " << iy << " " << dMap[ix][iy] << endl; } } @@ -100,12 +92,12 @@ class moench02ModuleData : public slsDetectorData { ~moench02ModuleData() {if (buff) delete [] buff; if (oldbuff) delete [] oldbuff; }; - int getFrameNumber(char *buff){ - return ((*(int*)buff)&(0xffffff00))>>8;; - }; int getPacketNumber(char *buff){ - return (*(int*)buff)&0xff; + int np=(*(int*)buff)&0xff; + if (np==0) + np=40; + return np; }; @@ -166,64 +158,6 @@ class moench02ModuleData : public slsDetectorData { - char *readNextFrame(ifstream &filebin) { - - - if (oldbuff) - delete [] oldbuff; - - oldbuff=buff; - - char *data=new char[1280*40]; - char p[1286]; - - int fn, fnum=-1, np=0, pnum=-1; - - if (filebin.is_open()) { - while (filebin.read(p,4)) { //read header - pnum=getPacketNumber(p); - fn=getFrameNumber(p); - - if (pnum<0 || pnum>39) { - cout << "Bad packet number " << pnum << " frame "<< fn << endl; - continue; - } - - if (pnum==0) - pnum=40; - - if (pnum==1) { - fnum=fn; - 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; - np=0; - } - filebin.read(data+1280*(pnum-1),1280); //readdata - np++; - filebin.read(p,2); //readfooter - - if (np==40) - if (pnum==40) - break; - else - cout << "Too many packets for this frame! "<< fnum << " " << pnum << endl; - - } - } - if (np<40) { - if (np>0) - cout << "Too few packets for this frame! "<< fnum << " " << pnum << endl; - delete [] data; - data=NULL; - } - buff=data; - return data; - }; - void addToPedestal(double val, int ix, int iy){stat[iy][ix].Calc(val);}; double getPedestal(int ix, int iy){return stat[iy][ix].Mean();}; @@ -237,13 +171,13 @@ class moench02ModuleData : public slsDetectorData { double getChannelShort(int ix, int iy, double hc=0, double tc=0) { - double val=slsDetectorData::getChannelShort(buff,ix,iy); + double val=slsDetectorData::getChannel(buff,ix,iy); if (ix>0 && hc!=0) - val-=hc*slsDetectorData::getChannelShort(buff,ix-1,iy); + val-=hc*slsDetectorData::getChannel(buff,ix-1,iy); if (oldbuff && tc!=0) - val+=tc*slsDetectorData::getChannelShort(oldbuff,ix,iy); + val+=tc*slsDetectorData::getChannel(oldbuff,ix,iy); return val; diff --git a/slsDetectorCalibration/moenchReadData.C b/slsDetectorCalibration/moenchReadData.C index 2095e380e..468817e45 100644 --- a/slsDetectorCalibration/moenchReadData.C +++ b/slsDetectorCalibration/moenchReadData.C @@ -60,7 +60,7 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb double oldval; int nf=0; - moench02ModuleData::eventType thisEvent=moench02ModuleData::PEDESTAL; + eventType thisEvent=PEDESTAL; // int iframe; // double *data, ped, sigma; @@ -162,7 +162,7 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb - if (thisEvent==moench02ModuleData::PEDESTAL) { + if (thisEvent==PEDESTAL) { if (cmsub && nf>1000) decoder->addToPedestal(decoder->getChannelShort(ix, iy, hc, tc)-decoder->getCommonMode(ix,iy), ix, iy); else @@ -186,7 +186,7 @@ Add here the function that you want to call: fill histos, make trees etc. // if (nf%1000==0 && ix==20 && iy==20) cout << " val="<< decoder->getClusterElement(0,0)<< endl; - if (thisEvent==moench02ModuleData::PHOTON_MAX ) { + if (thisEvent==PHOTON_MAX ) { #ifdef MY_DEBUG if (iev%100000==0) { cout << "Event " << iev << " Frame "<< nf << endl; diff --git a/slsDetectorCalibration/singlePhotonDetector.h b/slsDetectorCalibration/singlePhotonDetector.h new file mode 100644 index 000000000..2c226a69e --- /dev/null +++ b/slsDetectorCalibration/singlePhotonDetector.h @@ -0,0 +1,225 @@ +#ifndef SINGLEPHOTONDETECTOR_H +#define SINGLEPHOTONDETECTOR_H +#include "slsDetectorData.h" + +#include "MovingStat.h" +#include "single_photon_hit.h" + +#include + +using namespace std; + + + enum eventType { + PEDESTAL, + NEIGHBOUR, + PHOTON, + PHOTON_MAX, + NEGATIVE_PEDESTAL, + }; + + + +template +class singlePhotonDetector { + + + public: + + + /** + + 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) + + + */ + + + singlePhotonDetector(slsDetectorData *d, int csize=3, double nsigma=5) : det(d), clusterSize(csize), clusterSizeY(csize),nSigma(nsigma) { + + + + det->getDetectorSize(nx,ny); + + + + stat=new MovingStat*[ny]; + for (int i=0; i=0 && ix=0 && iy=0 && ix=0 && iy=0 && ix=0 && iy0) nSigma=n; return nSigma;} + + int setClusterSize(int n=-1){ + if (n>0 && n!=clusterSize) { + if (n%2==0) + n+=1; + clusterSize=n; + delete cluster; + if (ny>1) + clusterSizeY=clusterSize; + + cluster=new single_photon_hit(clusterSize,clusterSizeY); + } + return clusterSize; + } + + + + eventType getEventType(char *data, int ix, int iy, double hcorr=0, double tcorr=0, int sign=1, int cm=0) { + + eventType ret=PEDESTAL; + double tot=0, v, max=0, sig; + + + sig=getPedestalRMS(ix,iy); + + cluster->x=ix; + cluster->y=iy; + cluster->rms=sig; + cluster->ped=getPedestal(ix,iy); + + + for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) { + for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) { + if ((iy+ir)>=0 && (iy+ir)<160 && (ix+ic)>=0 && (ix+ic)<160) { + v=sign*(det->getValue(data, ix+ic, iy+ir)-getPedestal(ix+ic,iy+ir)); + // if (cm) + // v-=sign*getCommonMode(ix+ic, iy+ic); + cluster->set_data(ic, ir, v); + tot+=v; + if (v>max) { + max=v; + } + if (ir==0 && ic==0) { + if (v>nSigma*getPedestalRMS(ix,iy)) { + ret=PHOTON; + } else if (v<-nSigma*getPedestalRMS(ix,iy)) + ret=NEGATIVE_PEDESTAL; + } + } + } + } + + if (ret!=PHOTON && tot>3*nSigma*sig) + ret=NEIGHBOUR; + else if (ret==PHOTON) { + if (cluster->get_data(0,0)>=max) + ret=PHOTON_MAX; + } + + return ret; + + }; + + + + + + + + double getClusterElement(int ic, int ir, int cmsub=0){return cluster->get_data(ic,ir);}; + + + private: + + slsDetectorData *det; + int nx, ny; + + + + MovingStat **stat; + double nSigma; + int clusterSize, clusterSizeY; + single_photon_hit *cluster; + + MovingStat cmStat[4]; //stores the common mode average per supercolumn + double cmPed[4]; // stores the common mode for this frame per supercolumn + double nCm[4]; // stores the number of pedestals to calculate the cm per supercolumn + int cx, cy; +}; + + + + + + + + /////////////////////////////////////////////////////// DONE UNTIL HERE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + /* void calculateCommonMode(int xmin=0, int xmax=160, int ymin=0, int ymax=160, double hc=0, double tc=0) { */ + +/* int scmin=xmin/40; */ +/* int scmax=(xmax-1)/40+1; */ +/* int isc=0, ix, iy; */ + +/* for (isc=scmin; isc0) */ +/* cmStat[isc].Calc(cmPed[isc]/nCm[isc]); */ + +/* // cout << "SC " << isc << nCm[isc] << " " << cmPed[isc]/nCm[isc] << " " << cmStat[isc].Mean() << endl; */ + +/* } */ + + +/* }; */ + +/* double getCommonMode(int ix, int iy) { */ +/* int isc=ix/40; */ +/* if (nCm[isc]>0) { */ +/* /\* if (ix==20 && iy==20) *\/ */ +/* /\* cout << cmPed[isc] << " " << nCm[isc] << " " << cmStat[isc].Mean() << " " << cmPed[isc]/nCm[isc]-cmStat[isc].Mean() << endl; *\/ */ +/* return cmPed[isc]/nCm[isc]-cmStat[isc].Mean(); */ +/* } else { */ +/* /\* if (ix==20 && iy==20) *\/ */ +/* /\* cout << "No common mode!" << endl; *\/ */ + +/* return 0; */ +/* } */ +/* }; */ + + + +#endif diff --git a/slsDetectorCalibration/single_photon_hit.h b/slsDetectorCalibration/single_photon_hit.h new file mode 100644 index 000000000..278c9f691 --- /dev/null +++ b/slsDetectorCalibration/single_photon_hit.h @@ -0,0 +1,44 @@ +#ifndef SINGLE_PHOTON_HIT_H +#define SINGLE_PHOTON_HIT_h + +typedef double double32_t; +typedef float float32_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 { + + public: + single_photon_hit(int nx, int ny=1): dx(nx), dy(ny) {data=new double[dx*dy];}; + ~single_photon_hit(){delete [] data;}; + 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 read(FILE *myFile) {fread((void*)this, 1, 3*sizeof(int)+2*sizeof(double), myFile); fread((void*)data, 1, dx*dy*sizeof(double), myFile);}; + void set_data(double v, int ix, int iy=0){data[(iy+dy/2)*dx+ix+dx/2]=v;}; + 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 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 */ + double *data; /**< data size */ + const int dx; + const int dy; +}; + + + +#endif diff --git a/slsDetectorCalibration/slsDetectorData.h b/slsDetectorCalibration/slsDetectorData.h index a8ffb8513..527e4e9bb 100644 --- a/slsDetectorCalibration/slsDetectorData.h +++ b/slsDetectorCalibration/slsDetectorData.h @@ -6,8 +6,13 @@ using namespace std; + +template class slsDetectorData { + + public: + /** Constructor (no error checking if datasize and offsets are compatible!) @@ -19,12 +24,9 @@ class slsDetectorData { */ - slsDetectorData(int npx, int npy=1, int ds=-1, int **dMap=NULL, int **dMask=NULL, int **dROI=NULL): nx(npx), ny(npy) { + slsDetectorData(int npx, int npy, int np, int psize, int **dMap=NULL, int **dMask=NULL, int **dROI=NULL): nx(npx), ny(npy), nPackets(np), packetSize(psize) { - if (ds<=0) - dataSize=nx*ny; - else - dataSize=ds; + dataMask=new int*[ny]; for(int i = 0; i < ny; i++) { @@ -46,7 +48,7 @@ class slsDetectorData { }; - ~slsDetectorData() { + virtual ~slsDetectorData() { for(int i = 0; i < ny; i++) { delete [] dataMap[i]; delete [] dataMask[i]; @@ -100,11 +102,14 @@ class slsDetectorData { }; - int setBad(int ix, int iy, int i=1) {dataROIMask[iy][ix]=i; return isGood(ix,iy);}; + int setGood(int ix, int iy, int i=1) {dataROIMask[iy][ix]=i; return isGood(ix,iy);}; int isGood(int ix, int iy) {return dataROIMask[iy][ix];}; + void getDetectorSize(int &npx, int &npy){npx=nx; npy=ny;}; + + /** Returns the value of the selected channel for the given dataset (no error checking if number of pixels are compatible!) @@ -117,34 +122,117 @@ class slsDetectorData { */ - inline int getChannel(char *data, int ix, int iy=1) { - return (*(data+dataMap[iy][ix]))^dataMask[iy][ix]; - }; - - /** - - Returns the value of the selected channel for the given dataset (no error checking if number of pixels are compatible!) - \param data pointer to the dataset (including headers etc) - \param ix pixel numb er in the x direction - \param iy pixel number in the y direction - - \returns data for the selected channel, with inversion if required - - */ - - inline u_int16_t getChannelShort(char *data, int ix, int iy=1) { - u_int16_t m=dataMask[iy][ix], d=*(u_int16_t*)(data+dataMap[iy][ix]); - // cout << ix << " " << iy << " " << dataMap[ix][iy] << endl; - // return (*(dd+dataMap[ix][iy]))^((u_int16_t)dataMask[ix][iy]); + virtual dataType getChannel(char *data, int ix, int iy=0) { + dataType m=dataMask[iy][ix], d=*((dataType*)(data+dataMap[iy][ix])); return d^m; }; - + + 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;}; + virtual int getPacketNumber(char *buff) {return (*(int*)buff)&0xff;}; + + + 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; + + return retval; + }; + + + virtual char *readNextFrame(ifstream &filebin) { + + +/* if (oldbuff) */ +/* delete [] oldbuff; */ + +/* oldbuff=buff; */ + + // char *p=new char[1286]; + char *data=new char[packetSize*nPackets]; + 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; ipnPackets) { + cout << "*******too many packets!!!!!!!!!!" << endl; + delete [] data; + return NULL; + } else { + // memcpy(data+np*1286,p,1286); + np++; + } + + } + } + delete [] data; + return NULL; + }; + private: const int nx; /**< Number of pixels in the x direction */ const int ny; /**< Number of pixels in the y direction */ - int dataSize; /**< Size of a dataset */ + const int nPackets; /**