diff --git a/slsDetectorCalibration/moench02ModuleData.h b/slsDetectorCalibration/moench02ModuleData.h index 5009a8686..7dc2d1bdf 100644 --- a/slsDetectorCalibration/moench02ModuleData.h +++ b/slsDetectorCalibration/moench02ModuleData.h @@ -2,13 +2,31 @@ #define MOENCH02MODULEDATA_H #include "slsDetectorData.h" +#include "RunningStat.h" +#include "MovingStat.h" + #include #include 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 (int ip=0; ip<10; ip++) { + for (isc=0; isc<4; isc++) { + for (ip=0; ip<10; ip++) { - for (int ir=0; ir<16; ir++) { - for (int ic=0; ic<40; ic++) { + for (ir=0; ir<16; ir++) { + 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 (int iy=0; iy<160; iy++) - dMask[iy][ix]=0x7fff; + for (ix=0; ix<120; ix++) { + for (iy=0; iy<160; iy++) + dMask[iy][ix]=0x3fff; } - for (int ix=120; ix<160; ix++) { - for (int iy=0; iy<160; iy++) + for (ix=120; ix<160; ix++) { + 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; + + }; diff --git a/slsDetectorCalibration/moenchReadData.C b/slsDetectorCalibration/moenchReadData.C index 7d9f425e3..078d5931a 100644 --- a/slsDetectorCalibration/moenchReadData.C +++ b/slsDetectorCalibration/moenchReadData.C @@ -14,26 +14,33 @@ #include #include #include -#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; irunreadNextFrame(filebin))) { - for (ix=0; ix<160; ix++) - for (iy=0; iy<160; iy++) { + for (ix=xmin-1; ix1000) { - 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)); - + + if (nf>100) { + thisEvent= decoder->getEventType(ix, iy, hc, tc, 1); + } + + if (thisEvent==moench02ModuleData::PEDESTAL) { + decoder->addToPedestal(decoder->getClusterElement(0,0), 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;