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

@@ -1,7 +1,13 @@
-  class MovingStat
+#ifndef MOVINGSTAT_H
+#define MOVINGSTAT_H 
+
+#include <math.h>
+
+
+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

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<uint16_t> {
  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<uint16_t>(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<uint16_t>::getChannel(buff,ix,iy);
 
     if (ix>0 && hc!=0)
-      val-=hc*slsDetectorData::getChannelShort(buff,ix-1,iy);
+      val-=hc*slsDetectorData<uint16_t>::getChannel(buff,ix-1,iy);
     
     if (oldbuff && tc!=0) 
-      val+=tc*slsDetectorData::getChannelShort(oldbuff,ix,iy);
+      val+=tc*slsDetectorData<uint16_t>::getChannel(oldbuff,ix,iy);
 
 
     return val;

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;

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 <iostream>
+
+using namespace std;
+
+
+  enum eventType {
+    PEDESTAL,
+    NEIGHBOUR,
+    PHOTON,
+    PHOTON_MAX,
+    NEGATIVE_PEDESTAL,
+  };
+
+
+
+template <class dataType>
+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<dataType> *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<ny; i++)
+      stat[i]=new MovingStat;
+   
+    if (ny==1)
+      clusterSizeY=1;
+
+    cluster=new single_photon_hit(clusterSize,clusterSizeY);
+    
+  };
+    
+    ~singlePhotonDetector() {delete cluster; delete [] stat;};
+  
+
+
+    void addToPedestal(double val, int ix, int iy){ if (ix>=0 && ix<nx && iy>=0 && iy<ny) stat[iy][ix].Calc(val);};
+    double getPedestal(int ix, int iy){if (ix>=0 && ix<nx && iy>=0 && iy<ny) return stat[iy][ix].Mean(); else return -1;};
+    double getPedestalRMS(int ix, int iy){if (ix>=0 && ix<nx && iy>=0 && iy<ny) return stat[iy][ix].StandardDeviation();else return -1;};
+  
+    double setNSigma(double n=-1){if (n>0) 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<dataType> *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; isc<scmax; isc++) { */
+/*       nCm[isc]=0; */
+/*       cmPed[isc]=0; */
+/*     } */
+    
+    
+
+
+/*     for (ix=xmin-1; ix<xmax+1; ix++) { */
+/*       isc=ix/40; */
+/*       for (iy=ymin-1; iy<ymax+1; iy++) { */
+/* 	if (isGood(ix,iy)) { */
+/* 	  if (getEventType(ix, iy, hc, tc, 1,0)==PEDESTAL) { */
+/* 	    cmPed[isc]+=getChannelShort(ix, iy, hc, tc); */
+/* 	    nCm[isc]++; */
+/* 	  } */
+/* 	} */
+	
+/*       } */
+
+/*     } */
+
+    
+/*     for (isc=scmin; isc<scmax; isc++) { */
+
+/*       if (nCm[isc]>0) */
+/* 	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

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

slsDetectorCalibration/slsDetectorData.h

@@ -6,8 +6,13 @@
 
 using namespace std;
 
+
+template <class dataType>
 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) {
+  virtual dataType getChannel(char *data, int ix, int iy=0) {
-    return (*(data+dataMap[iy][ix]))^dataMask[iy][ix];
+    dataType m=dataMask[iy][ix], d=*((dataType*)(data+dataMap[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]);
     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; 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;
+  };
+
 
 
  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; /**<number of UDP packets 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 **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)   */