generalized several funcs for hit finding and looping over meonch data

git-svn-id: file:///afs/psi.ch/project/sls_det_software/svn/slsDetectorCalibration@6 113b152e-814d-439b-b186-022a431db7b5

slsDetectorCalibration/moench02ModuleData.h

@@ -2,13 +2,31 @@
 #define  MOENCH02MODULEDATA_H
 #include "slsDetectorData.h"
 
+#include "RunningStat.h"
+#include "MovingStat.h"
+
 #include <iostream>
 #include <fstream>
 
 using namespace std;
 
+
+
 class moench02ModuleData : public slsDetectorData {
  public:
+
+
+  enum eventType {
+    PEDESTAL,
+    NEIGHBOUR,
+    PHOTON,
+    PHOTON_MAX,
+    NEGATIVE_PEDESTAL,
+  };
+
+
+
+
   /**
 
      Constructor (no error checking if datasize and offsets are compatible!)
@@ -22,30 +40,38 @@ class moench02ModuleData : public slsDetectorData {
   */
   
 
-  moench02ModuleData(): slsDetectorData(160, 160, 1286*40) {
+  moench02ModuleData(int nbg=1000): slsDetectorData(160, 160, 1286*40), nSigma(5), buff(NULL), oldbuff(NULL), pedSub(1) {
     int headerOff=0;
     int footerOff=1284;
     int dataOff=4;
 
-    int ix, iy;
+    int ix, iy, ir, ic, i, isc, ip;
 
     int **dMask;
     int **dMap;
 
+    for (ir=0; ir<160; ir++) {
+      for (ic=0; ic<160; ic++) {
+	stat[ir][ic].Clear();
+	stat[ir][ic].SetN(nbg);
+      }
+    }
+
+
     dMask=new int*[160];
-    for(int i = 0; i < 160; i++) {
+    for(i = 0; i < 160; i++) {
       dMask[i] = new int[160];
     }
     dMap=new int*[160];
-    for(int i = 0; i < 160; i++) {
+    for(i = 0; i < 160; i++) {
       dMap[i] = new int[160];
     }
 
-    for (int isc=0; isc<4; isc++) {
+    for (isc=0; isc<4; isc++) {
-      for (int ip=0; ip<10; ip++) {
+      for (ip=0; ip<10; ip++) {
 
-	for (int ir=0; ir<16; ir++) {
+	for (ir=0; ir<16; ir++) {
-	  for (int ic=0; ic<40; ic++) {
+	  for (ic=0; ic<40; ic++) {
 
 	    ix=isc*40+ic;
 	    iy=ip*16+ir;
@@ -57,12 +83,12 @@ class moench02ModuleData : public slsDetectorData {
       }
     }
 
-    for (int ix=0; ix<120; ix++) {
+    for (ix=0; ix<120; ix++) {
-      for (int iy=0; iy<160; iy++)
+      for (iy=0; iy<160; iy++)
-	dMask[iy][ix]=0x7fff;
+	dMask[iy][ix]=0x3fff;
     }
-    for (int ix=120; ix<160; ix++) {
+    for (ix=120; ix<160; ix++) {
-      for (int iy=0; iy<160; iy++)
+      for (iy=0; iy<160; iy++)
 	dMask[iy][ix]=0x0;
     }
 
@@ -70,9 +96,10 @@ class moench02ModuleData : public slsDetectorData {
     setDataMap(dMap);
     setDataMask(dMask);
   };
-
+    
+    ~moench02ModuleData() {if (buff) delete [] buff; if (oldbuff) delete [] oldbuff; };
   
-  int getFrameNumber(char *buff){
+    int getFrameNumber(char *buff){
     return ((*(int*)buff)&(0xffffff00))>>8;;
   };
 
@@ -82,7 +109,13 @@ class moench02ModuleData : public slsDetectorData {
 
 
   char *readNextFrame(ifstream &filebin) {
-    char *buff=new char[1280*40];
+    
+    if (oldbuff)
+      delete [] oldbuff;
+
+    oldbuff=buff;
+    
+    char *data=new char[1280*40];
     char p[1286];
     
     int fn, fnum=-1,  np=0, pnum=-1;
@@ -110,7 +143,7 @@ class moench02ModuleData : public slsDetectorData {
 	    cout << " **Incomplete frame number " << fnum << " pnum " << pnum << " " << getFrameNumber(p) << endl;
 	  np=0;
 	}
-	filebin.read(buff+1280*(pnum-1),1280); //readdata
+	filebin.read(data+1280*(pnum-1),1280); //readdata
 	np++;
 	filebin.read(p,2); //readfooter
 
@@ -125,13 +158,108 @@ class moench02ModuleData : public slsDetectorData {
     if (np<40) {
       if (np>0)
 	cout << "Too few packets for this frame! "<< fnum << " " << pnum << endl;
-      delete [] buff;
+      delete [] data;
       buff=NULL;
     }
-    return buff;
+    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();};
+  double getPedestalRMS(int ix, int iy){return stat[iy][ix].StandardDeviation();};
+  double setNSigma(double n){if (n>0) nSigma=n; return nSigma;}
+  char *getBuff(){return buff;}
+  char *getOldBuff(){return oldbuff;}
+
+  int setPedestalSubstraction(int i=-1){if (i==0) pedSub=0; if (i==1) pedSub=1; return pedSub;};
+
+
+  double getChannelShort(int ix, int iy, double hc=0, double tc=0) {
+    
+    double val=slsDetectorData::getChannelShort(buff,ix,iy);
+
+    if (ix>0 && hc!=0)
+      val-=hc*slsDetectorData::getChannelShort(buff,ix-1,iy);
+    
+    if (oldbuff && tc!=0)
+      val+=tc*slsDetectorData::getChannelShort(oldbuff,ix,iy);
+
+
+    return val;
+  };
+  
+
+
+  eventType getEventType(int ix, int iy, double hcorr=0, double tcorr=0, int sign=1) {
+
+ /*    PEDESTAL, */
+/*     NEIGHBOUR, */
+/*     PHOTON, */
+/*     PHOTON_MAX, */
+ /*    NEGATIVE_PEDESTAL */
+      eventType ret=PEDESTAL;
+      double tot=0, v, max=0, sig;
+      
+
+      sig=getPedestalRMS(ix,iy);
+
+      for (int ir=-1; ir<2; ir++) {
+	for (int ic=-1; ic<2; ic ++) {
+	  if ((iy+ir)>=0 && (iy+ir)<160 && (ix+ic)>=0 && (ix+ic)<160) {
+	    v=getChannelShort(ix+ic, iy+ir, hcorr, tcorr);
+	    if (pedSub)
+	      v=sign*(v-getPedestal(ix+ic,iy+ir));
+	    cluster[ic+1][ir+1]=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[1][1]>=max)
+	  ret=PHOTON_MAX;
+      }
+      
+   
+
+
+      return ret;
+
+  };
+
+
+  double getClusterElement(int ic, int ir){return cluster[ic+1][ir+1];};
+  double  *getCluster(){return &cluster[0][0];};
+
+
+
+
+ private:
+  
+  MovingStat stat[160][160];
+  double nSigma;
+  double cluster[3][3];
+  char *buff;
+  char *oldbuff;
+  int pedSub;
+
+
 };
 
 

slsDetectorCalibration/moenchReadData.C

@@ -14,26 +14,33 @@
 #include <list>
 #include <queue>
 #include <fstream>
-#include "RunningStat.h"
-#include "MovingStat.h"
 #include "moench02ModuleData.h"
 
 using namespace std;
 
+#define NC 160
+#define NR 160
 
-TH2F *moenchReadData(char *fformat, int runmin, int runmax, int nbins=1500, int hmin=-500, int hmax=1000, int sign=1) {
 
+
+TH2F *moenchReadData(char *fformat, int runmin, int runmax, int nbins=1500, int hmin=-500, int hmax=1000, int sign=1, int nsigma=5, double hc=0, double tc=0, int xmin=0, int xmax=NC, int ymin=0, int ymax=NR) {
+  
   moench02ModuleData *decoder=new moench02ModuleData();
   char *buff;
+  char *oldbuff=NULL;
   char fname[10000];
-  
+  double oldval;
   int nf=0;
 
+
+  moench02ModuleData::eventType thisEvent=moench02ModuleData::PEDESTAL;
+
+
   TH2F *h2=NULL;
-  h2=new TH2F("h2",fformat,nbins,hmin-0.5,hmax-0.5,160*160,-0.5,160*160-0.5);
+  h2=new TH2F("h2",fformat,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
-  
+
   int val, dum;
-  double me, sig, tot;
+  double me=0, sig, tot, vv;
 
   MovingStat stat[160][160];
 
@@ -41,15 +48,11 @@ TH2F *moenchReadData(char *fformat, int runmin, int runmax, int nbins=1500, int
 
   int nbg=1000;
 
-  int ix, iy, ir, ic;
+  int ix=20, iy=20, ir, ic;
 
+  // 6% x-talk from previous pixel
+  // 12% x-talk from previous frame
 
-  for (ir=0; ir<160; ir++) {
-    for (ic=0; ic<160; ic++) {
-      stat[ir][ic].Clear();
-      stat[ir][ic].SetN(nbg);
-    }
-  }
 
   for (int irun=runmin; irun<runmax; irun++) {
     sprintf(fname,fformat,irun);
@@ -58,38 +61,25 @@ TH2F *moenchReadData(char *fformat, int runmin, int runmax, int nbins=1500, int
    
     while ((buff=decoder->readNextFrame(filebin))) {
 	
-	for (ix=0; ix<160; ix++)
+      for (ix=xmin-1; ix<xmax+1; ix++)
-	  for (iy=0; iy<160; iy++) {
+	for (iy=ymin-1; iy<ymax+1; iy++) {
 
-	    dum=0; //no hit
+	  
-	    tot=0;
+	  if (nf>100) {
-
+	    thisEvent= decoder->getEventType(ix, iy, hc, tc, 1);
-	    if (nf>1000) {
+	  } 
-	      me=stat[iy][ix].Mean();
+		
-	      sig=stat[iy][ix].StandardDeviation();
+	  if (thisEvent==moench02ModuleData::PEDESTAL) {	
-	      val=sign*decoder->getChannelShort(buff,ix,iy)-me;
+	    decoder->addToPedestal(decoder->getClusterElement(0,0), ix, iy);
-	      h2->Fill(val,ix*160+iy);
-
-	      dum=0; //no hit
-	      tot=0;
-
-	      for (ir=-1; ir<2; ir++)
-		for (ic=-1; ic<2; ic++){
-		  if ((ix+ic)>=0 && (ix+ic)<160 && (iy+ir)>=0 && (iy+ir)<160) {
-		    if (decoder->getChannelShort(buff,ix+ic,iy+ir)>(stat[iy+ir][ix+ic].Mean()+3.*stat[iy+ir][ix+ic].StandardDeviation())) dum=1; //is a hit or neighbour is a hit!
-		    tot+=decoder->getChannelShort(buff,ix+ic,iy+ir)-stat[iy+ir][ix+ic].Mean();
-		  }
-		}
-
-	      if (tot>3.*sig) dum=3; //discard events where sum of the neighbours is too large.
-
-	      if (val<(me-3.*sig)) dum=2; //discard negative events!
-	    }
-	    if (nf<1000 || dum==0) 
-	      stat[iy][ix].Calc(sign*decoder->getChannelShort(buff,ix,iy));
-	   
 	  }
-      delete [] buff;
+	    
+	      
+	      /***********************************************************
+Add  here the function that you want to call: fill histos, make trees etc.
+	      ************************************************************/
+	      // getClusterElement to access quickly the photon and the neighbours
+	  
+	}
       cout << "=" ;
       nf++;
     }
@@ -99,7 +89,11 @@ TH2F *moenchReadData(char *fformat, int runmin, int runmax, int nbins=1500, int
     else
       cout << "could not open file " << fname << endl;
   }
+  
+  
 
+  
+  
   delete decoder;
   cout << "Read " << nf << " frames" << endl;
   return h2;