moench data classes

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

slsDetectorCalibration/energyCalibration.cpp

@@ -22,14 +22,65 @@
 using namespace std;
 
 #ifdef MYROOT
+
+Double_t energyCalibrationFunctions::pedestal(Double_t *x, Double_t *par) { 
+    return par[0]-par[1]*sign*x[0];                                      
+}
+
+
 Double_t energyCalibrationFunctions::gaussChargeSharing(Double_t *x, Double_t *par) { 
   Double_t f, arg=0;
-  if (par[3]!=0) arg=(x[0]-par[2])/par[3];
+  if (par[3]!=0) arg=sign*(x[0]-par[2])/par[3];
-    f=par[4]*TMath::Exp(-1*arg*arg/2.); 
+    f=TMath::Exp(-1*arg*arg/2.); 
-    f=f+par[5]*(par[4]/2*(TMath::Erfc(sign*arg/(TMath::Sqrt(2.)))))+par[0]-par[1]*sign*x[0]; 
+    f=f+par[5]/2.*(TMath::Erfc(arg/(TMath::Sqrt(2.)))); 
-    return f;                                       
+    return par[4]*f+pedestal(x,par);                                       
 }
 
+Double_t energyCalibrationFunctions::gaussChargeSharingPixel(Double_t *x, Double_t *par) { 
+  Double_t f;
+  if (par[3]<=0 || par[2]*(*x)<=0 || par[5]<0 || par[4]<=0) return 0;
+
+  Double_t pp[3];
+
+  pp[0]=0;
+  pp[1]=par[2];
+  pp[2]=par[3];
+  
+
+  f=(par[5]-par[6]*(TMath::Log(*x/par[2])))*erfBox(x,pp);
+  f+=par[4]*TMath::Gaus(*x, par[2], par[3], kTRUE);
+  return f+pedestal(x,par);                                       
+}
+
+Double_t energyCalibrationFunctions::erfBox(Double_t *z, Double_t *par) {
+
+  
+  
+  Double_t m=par[0];
+  Double_t M=par[1];
+  
+  if (par[0]>par[1]) {
+    m=par[1];
+    M=par[0];
+  }
+  
+  if (m==M)
+    return 0;
+
+
+  if (par[2]<=0) {
+    if (*z>=m && *z<=M)
+      return 1./(M-m);
+    else
+      return 0;
+
+  }
+
+  return (TMath::Erfc((z[0]-M)/par[2])-TMath::Erfc((z[0]-m)/par[2]))*0.5/(M-m);
+
+}
+
+
 // basic erf function
 Double_t energyCalibrationFunctions::erfFunction(Double_t *x, Double_t *par) {	
   double arg=0;
@@ -155,6 +206,10 @@ Double_t energyCalibrationFunctions::spectrum(Double_t *x, Double_t *par) {
   return gaussChargeSharing(x,par);	
 }
 
+Double_t energyCalibrationFunctions::spectrumPixel(Double_t *x, Double_t *par) {
+  return gaussChargeSharingPixel(x,par);	
+}
+
 
 Double_t energyCalibrationFunctions::scurve(Double_t *x, Double_t *par) {
   return erfFunctionChargeSharing(x,par);
@@ -193,6 +248,9 @@ energyCalibration::energyCalibration() :
   fspectrum=new TF1("fspectrum",funcs,&energyCalibrationFunctions::spectrum,0,1000,6,"energyCalibrationFunctions","spectrum");	
   fspectrum->SetParNames("Background Pedestal","Background slope", "Peak position","Noise RMS", "Number of Photons","Charge Sharing Pedestal");
 
+  fspixel=new TF1("fspixel",funcs,&energyCalibrationFunctions::spectrumPixel,0,1000,7,"energyCalibrationFunctions","spectrumPixel");	
+  fspixel->SetParNames("Background Pedestal","Background slope", "Peak position","Noise RMS", "Number of Photons","Charge Sharing Pedestal","Corner");
+
 #endif
 
 

slsDetectorCalibration/energyCalibration.h

@@ -18,6 +18,7 @@
 #define MYROOT
 #endif
 
+#define MYROOT
 
 #ifdef MYROOT
 #include <TROOT.h>
@@ -80,6 +81,13 @@ class energyCalibrationFunctions {
   
 
 #ifdef MYROOT
+  /** 
+      Gaussian Function with charge sharing pedestal
+      par[0] is the absolute height of the background pedestal
+      par[1] is the slope of the background pedestal
+  */
+  Double_t pedestal(Double_t *x, Double_t *par);
+
   /** 
       Gaussian Function with charge sharing pedestal
       par[0] is the absolute height of the background pedestal
@@ -90,6 +98,16 @@ class energyCalibrationFunctions {
       par[5] is the fractional height of the charge sharing pedestal (scales with par[3])
   */
   Double_t gaussChargeSharing(Double_t *x, Double_t *par);
+  /** 
+      Gaussian Function with charge sharing pedestal
+      par[0] is the absolute height of the background pedestal
+      par[1] is the slope of the background pedestal
+      par[2] is the gaussian peak position
+      par[3] is the RMS of the gaussian (and of the pedestal)
+      par[4] is the height of the function
+      par[5] is the fractional height of the charge sharing pedestal (scales with par[3])
+  */
+  Double_t gaussChargeSharingPixel(Double_t *x, Double_t *par);
 
   /** 
       Basic erf function
@@ -98,7 +116,7 @@ class energyCalibrationFunctions {
       par[2] is the amplitude
   */
 Double_t erfFunction(Double_t *x, Double_t *par) ;
-
+ Double_t erfBox(Double_t *z, Double_t *par);
   /** Erf function with charge sharing slope
       par[0] is the  pedestal 
       par[1] is the  slope of the pedestal
@@ -128,13 +146,25 @@ Double_t erfFuncFluo(Double_t *x, Double_t *par);
   /** 
       static function Gaussian with charge sharing pedestal with the correct scan sign
       par[0] is the absolute height of the background pedestal
-      par[1] is the fractional height of the charge sharing pedestal (scales with par[3]
+      par[1]  is the  slope of the pedestal 
       par[2] is the gaussian peak position
       par[3] is the RMS of the gaussian (and of the pedestal)
       par[4] is the height of the function
+      par[5] is the fractional height of the charge sharing pedestal (scales with par[4]
   */
   Double_t spectrum(Double_t *x, Double_t *par);
 
+  /** 
+      static function Gaussian with charge sharing pedestal with the correct scan sign
+      par[0] is the absolute height of the background pedestal
+      par[1]  is the  slope of the pedestal 
+      par[2] is the gaussian peak position
+      par[3] is the RMS of the gaussian (and of the pedestal)
+      par[4] is the height of the function
+      par[5] is the fractional height of the charge sharing pedestal (scales with par[4]
+  */
+  Double_t spectrumPixel(Double_t *x, Double_t *par);
+
 
  /** Erf function with charge sharing slope with the correct scan sign
       par[0] is the  pedestal 
@@ -381,6 +411,8 @@ class energyCalibration  {
 
   TF1 *fspectrum; /**< function with which the spectrum will be fitted */
 
+  TF1 *fspixel; /**< function with which the spectrum will be fitted */
+
 #endif
 
   energyCalibrationFunctions *funcs;

slsDetectorCalibration/moench02ModuleData.h

@@ -0,0 +1,139 @@
+#ifndef MOENCH02MODULEDATA_H
+#define  MOENCH02MODULEDATA_H
+#include "slsDetectorData.h"
+
+#include <iostream>
+#include <fstream>
+
+using namespace std;
+
+class moench02ModuleData : public slsDetectorData {
+ 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)
+
+
+  */
+  
+
+  moench02ModuleData(): slsDetectorData(160, 160, 1286*40) {
+    int headerOff=0;
+    int footerOff=1284;
+    int dataOff=4;
+
+    int ix, iy;
+
+    int **dMask;
+    int **dMap;
+
+    dMask=new int*[160];
+    for(int i = 0; i < 160; i++) {
+      dMask[i] = new int[160];
+    }
+    dMap=new int*[160];
+    for(int i = 0; i < 160; i++) {
+      dMap[i] = new int[160];
+    }
+
+    for (int isc=0; isc<4; isc++) {
+      for (int ip=0; ip<10; ip++) {
+
+	for (int ir=0; ir<16; ir++) {
+	  for (int ic=0; ic<40; ic++) {
+
+	    ix=isc*40+ic;
+	    iy=ip*16+ir;
+
+	    dMap[iy][ix]=1280*(isc*10+ip)+2*ir*40+2*ic;
+	    // cout << ix << " " << iy << " " << dMap[ix][iy] << endl;
+	  }
+	}
+      }
+    }
+
+    for (int ix=0; ix<120; ix++) {
+      for (int iy=0; iy<160; iy++)
+	dMask[iy][ix]=0x7fff;
+    }
+    for (int ix=120; ix<160; ix++) {
+      for (int iy=0; iy<160; iy++)
+	dMask[iy][ix]=0x0;
+    }
+
+
+    setDataMap(dMap);
+    setDataMask(dMask);
+  };
+
+  
+  int getFrameNumber(char *buff){
+    return ((*(int*)buff)&(0xffffff00))>>8;;
+  };
+
+  int getPacketNumber(char *buff){
+    return (*(int*)buff)&0xff;
+  };
+
+
+  char *readNextFrame(ifstream &filebin) {
+    char *buff=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(buff+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 [] buff;
+      buff=NULL;
+    }
+    return buff;
+  };
+
+
+};
+
+
+
+#endif

slsDetectorCalibration/moenchReadData.C

@@ -0,0 +1,107 @@
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TPad.h>
+#include <TDirectory.h>
+#include <TEntryList.h>
+#include <TFile.h>
+#include <TMath.h>
+#include <TTree.h>
+#include <TChain.h>
+#include <THStack.h>
+#include <TCanvas.h>
+#include <stdio.h>
+#include <deque>
+#include <list>
+#include <queue>
+#include <fstream>
+#include "RunningStat.h"
+#include "MovingStat.h"
+#include "moench02ModuleData.h"
+
+using namespace std;
+
+
+TH2F *moenchReadData(char *fformat, int runmin, int runmax, int nbins=1500, int hmin=-500, int hmax=1000, int sign=1) {
+
+  moench02ModuleData *decoder=new moench02ModuleData();
+  char *buff;
+  char fname[10000];
+  
+  int nf=0;
+
+  TH2F *h2=NULL;
+  h2=new TH2F("h2",fformat,nbins,hmin-0.5,hmax-0.5,160*160,-0.5,160*160-0.5);
+  
+  int val, dum;
+  double me, sig, tot;
+
+  MovingStat stat[160][160];
+
+  ifstream filebin;
+
+  int nbg=1000;
+
+  int ix, iy, ir, ic;
+
+
+  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);
+    cout << "file name " << fname << endl;
+    filebin.open((const char *)(fname), ios::in | ios::binary);
+   
+    while ((buff=decoder->readNextFrame(filebin))) {
+	
+	for (ix=0; ix<160; ix++)
+	  for (iy=0; iy<160; iy++) {
+
+	    dum=0; //no hit
+	    tot=0;
+
+	    if (nf>1000) {
+	      me=stat[iy][ix].Mean();
+	      sig=stat[iy][ix].StandardDeviation();
+	      val=sign*decoder->getChannelShort(buff,ix,iy)-me;
+	      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;
+      cout << "=" ;
+      nf++;
+    }
+    cout << endl;
+    if (filebin.is_open())
+      filebin.close();	  
+    else
+      cout << "could not open file " << fname << endl;
+  }
+
+  delete decoder;
+  cout << "Read " << nf << " frames" << endl;
+  return h2;
+}
+

slsDetectorCalibration/slsDetectorData.h

@@ -1,8 +1,13 @@
 #ifndef SLSDETECTORDATA_H
 #define  SLSDETECTORDATA_H
 
-class slsDetectorData {
+#include <iostream>
+#include <fstream>
 
+using namespace std;
+
+class slsDetectorData {
+ public:
   /**
 
      Constructor (no error checking if datasize and offsets are compatible!)
@@ -23,26 +28,58 @@ class slsDetectorData {
     else
       dataSize=ds;
 
-    dataMap=new int[nx][ny];
+    dataMask=new int*[ny];
-    dataMask=new int[nx][ny];
+    for(int i = 0; i < ny; i++) {
+      dataMask[i] = new int[nx];
+    }
+    dataMap=new int*[ny];
+    for(int i = 0; i < ny; i++) {
+      dataMap[i] = new int[nx];
+    }
     
+    setDataMap(dMap);
+    setDataMask(dMask);
+
+  };
+
+  ~slsDetectorData() {
+    for(int i = 0; i < ny; i++) {
+    delete [] dataMap[i];
+    delete [] dataMask[i];
+    }
+    delete [] dataMap;
+    delete [] dataMask;
+  }
+
+  void setDataMap(int **dMap=NULL) {
+
+
     if (dMap==NULL) {
       for (int iy=0; iy<ny; iy++)
 	 for (int ix=0; ix<nx; ix++)
-	   dataMap[ix][iy]=iy*nx+ix;
+	   dataMap[iy][ix]=ix*ny+ix;
     } else {
       for (int iy=0; iy<ny; iy++)
 	 for (int ix=0; ix<nx; ix++)
-	   dataMap[ix][iy]=dMap[ix][iy];
+	   dataMap[iy][ix]=dMap[iy][ix];
-    }
-    if (dMap!=NULL) {
-
-      for (int iy=0; iy<ny; iy++)
-	 for (int ix=0; ix<nx; ix++)
-	   dataMask[ix][iy]=dMask[ix][iy];
     }
     
   };
+
+  void setDataMask(int **dMask=NULL){
+
+    if (dMask!=NULL) {
+
+      for (int iy=0; iy<ny; iy++)
+	 for (int ix=0; ix<nx; ix++)
+	   dataMask[iy][ix]=dMask[iy][ix];
+    }
+    
+  };
+
+
+
+
   /**
 
      Returns the value of the selected channel for the given dataset (no error checking if number of pixels are compatible!)
@@ -55,8 +92,8 @@ class slsDetectorData {
   */
 
 
-  int getChannel(int *data, int ix, int iy=1) {
+  int getChannel(char *data, int ix, int iy=1) {
-    return *(data+dataMap[ix][iy])^dataMask[ix][iy];
+    return (*(data+dataMap[iy][ix]))^dataMask[iy][ix];
   };
   
   /**
@@ -70,9 +107,11 @@ class slsDetectorData {
 
   */
 
-
+  u_int16_t getChannelShort(char  *data, int ix, int iy=1) {
-  int getChannel(u_int16_t  *data, int ix, int iy=1) {
+    u_int16_t m=dataMask[iy][ix], d=*(u_int16_t*)(data+dataMap[iy][ix]);
-    return *(data+dataMap[ix][iy])^((u_int16_t)dataMask[ix][iy]);
+    // cout << ix << " " << iy << " " << dataMap[ix][iy] << endl;
+    // return (*(dd+dataMap[ix][iy]))^((u_int16_t)dataMask[ix][iy]);
+    return d^m;
   };