diff --git a/slsDetectorCalibration/energyCalibration.cpp b/slsDetectorCalibration/energyCalibration.cpp index 6264182cf..df548f1cf 100644 --- a/slsDetectorCalibration/energyCalibration.cpp +++ b/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]; - f=par[4]*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]; - return f; + if (par[3]!=0) arg=sign*(x[0]-par[2])/par[3]; + f=TMath::Exp(-1*arg*arg/2.); + f=f+par[5]/2.*(TMath::Erfc(arg/(TMath::Sqrt(2.)))); + 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 diff --git a/slsDetectorCalibration/energyCalibration.h b/slsDetectorCalibration/energyCalibration.h index 3e7a9717f..904705aa3 100644 --- a/slsDetectorCalibration/energyCalibration.h +++ b/slsDetectorCalibration/energyCalibration.h @@ -18,6 +18,7 @@ #define MYROOT #endif +#define MYROOT #ifdef MYROOT #include @@ -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; diff --git a/slsDetectorCalibration/moench02ModuleData.h b/slsDetectorCalibration/moench02ModuleData.h new file mode 100644 index 000000000..5009a8686 --- /dev/null +++ b/slsDetectorCalibration/moench02ModuleData.h @@ -0,0 +1,139 @@ +#ifndef MOENCH02MODULEDATA_H +#define MOENCH02MODULEDATA_H +#include "slsDetectorData.h" + +#include +#include + +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 diff --git a/slsDetectorCalibration/moenchReadData.C b/slsDetectorCalibration/moenchReadData.C new file mode 100644 index 000000000..7d9f425e3 --- /dev/null +++ b/slsDetectorCalibration/moenchReadData.C @@ -0,0 +1,107 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#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; irunreadNextFrame(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; +} + diff --git a/slsDetectorCalibration/slsDetectorData.h b/slsDetectorCalibration/slsDetectorData.h index 709afdb5b..519014b7d 100644 --- a/slsDetectorCalibration/slsDetectorData.h +++ b/slsDetectorCalibration/slsDetectorData.h @@ -1,8 +1,13 @@ #ifndef SLSDETECTORDATA_H #define SLSDETECTORDATA_H -class slsDetectorData { +#include +#include +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[nx][ny]; + dataMask=new int*[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