diff --git a/examples/moench03_T1.config b/examples/moench03_T1.config index 4102605a7..d062ebb58 100644 --- a/examples/moench03_T1.config +++ b/examples/moench03_T1.config @@ -7,28 +7,28 @@ hostname bchip181+ ############################################# ### edit with hostname or 1Gbs IP address of your server ############################################ -rx_hostname mpc2608 +rx_hostname mpc2011 rx_tcpport 1954 ############################################# ### edit with 10 Gbs IP of your server ############################################ -udp_dstip 10.1.2.117 +udp_dstip 10.1.1.102 ############################################# ### edit with any number in the subnet of your server (first 3 numbers as above) ############################################ -udp_srcip 10.1.2.19 +udp_srcip 10.1.1.19 udp_dstport 32411 ############################################# ### edit with 10 Gbs IP of your server ############################################ -rx_zmqip 10.1.2.117 +rx_zmqip 10.1.1.102 rx_zmqport 50003 ############################################# ### edit with 1 Gbs IP of PC where you will run the GUI ############################################ -zmqip 129.129.202.86 +zmqip 129.129.202.136 zmqport 50001 @@ -39,8 +39,8 @@ rx_datastream 1 -frames 1 -period 0.001 +frames 100000 +period 0.0006 ############################################# ### edit with directory you want to write to ############################################ @@ -48,8 +48,8 @@ fpath /mnt/moench_data/scratch/ fwrite 0 rx_datastream 1 -rx_jsonpara frameMode newpedestal -rx_jsonpara detectorMode analog +rx_jsonpara frameMode frame +rx_jsonpara detectorMode counting rx_discardpolicy discardpartial @@ -59,3 +59,5 @@ powerchip 1 vhighvoltage 90 #adcreg 0x14 0x40 + + diff --git a/slsDetectorCalibration/MovingStat.h b/slsDetectorCalibration/MovingStat.h index 46917c98e..31fb59c59 100755 --- a/slsDetectorCalibration/MovingStat.h +++ b/slsDetectorCalibration/MovingStat.h @@ -14,27 +14,27 @@ class MovingStat /** constructor \param nn number of samples parameter to be used */ - MovingStat(int nn=1000) : n(nn), m_n(0) {} - + MovingStat(int nn) : n(nn), m_n(0), m_newM(0),m_newM2(0) {} + // void setPointers(double *me, double *va) {mean=me; var=va;} /** clears the moving average number of samples parameter, mean and standard deviation */ - void Clear() - { + void Clear() + { m_n = 0; m_newM=0; m_newM2=0; - } + } /** clears the moving average number of samples parameter, mean and standard deviation */ void Set(double val, double rms=0, int m=-1) { - if (m>0) m_n = m; else m_n = n; + if (m>=0) m_n = m; else m_n = n; m_newM=val*m_n; - // cout << "set " << val << " " << m << " " << m_n << " " << m_newM << endl; SetRMS(rms); + // cout << "set " << val << " " << m << " " << m_n << " " << m_newM << endl; } /** clears the moving average number of samples parameter, mean and standard deviation @@ -42,12 +42,15 @@ class MovingStat void SetRMS(double rms) { if (rms<=0) { - m_newM2=m_newM*m_newM/n; + if (m_n>0) + m_newM2=m_newM*m_newM/m_n; + else + m_newM2=0; //m_n=0; } else { - if (m_n>0) + if (m_n>0) { m_newM2=(m_n*rms*rms+m_newM*m_newM/m_n); - else { + } else { m_newM2=(m_n*rms*rms+m_newM*m_newM/n); m_n=0; } @@ -102,13 +105,14 @@ class MovingStat m_newM = x; m_newM2 = x*x; m_n++; - } else { - m_newM = m_newM + x - m_newM/m_n; - m_newM2 = m_newM2 + x*x - m_newM2/m_n; - } - + } + else { + m_newM = m_newM + x - m_newM/m_n; + m_newM2 = m_newM2 + x*x - m_newM2/m_n; + } + } - + /** returns the current number of elements of the moving average \returns returns the current number of elements of the moving average */ @@ -122,7 +126,9 @@ class MovingStat inline double Mean() const { // cout << "get " << m_n << " " << m_newM << " " << m_newM/m_n << endl; - return (m_n > 0) ? m_newM/m_n : 0.0; + + + return (m_n > 0) ? m_newM/m_n : 0.0; } /** returns the squared mean, 0 if no elements are inside @@ -138,7 +144,7 @@ class MovingStat */ inline double Variance() const { - return ( (m_n > 1) ? (M2()-Mean()*Mean()) : 0.0 ); + return (m_n > 0) ? m_newM2/m_n-m_newM/m_n*m_newM/m_n : 0.0; } /** returns the standard deviation, 0 if no elements are inside @@ -146,7 +152,8 @@ class MovingStat */ inline double StandardDeviation() const { - return ( (Variance() > 0) ? sqrt( Variance() ) : -1 ); + + return sqrt(Variance());// } private: diff --git a/slsDetectorCalibration/analogDetector.h b/slsDetectorCalibration/analogDetector.h index 984271305..56b3feb29 100644 --- a/slsDetectorCalibration/analogDetector.h +++ b/slsDetectorCalibration/analogDetector.h @@ -68,10 +68,15 @@ template class analogDetector { det->getDetectorSize(nx,ny); stat=new pedestalSubtraction*[ny]; + /* pedMean=new double*[ny]; */ + /* pedVariance=new double*[ny]; */ for (int i=0; i class analogDetector { /** destructor. Deletes the pdestalSubtraction array and the image */ - virtual ~analogDetector() {for (int i=0; i class analogDetector { stat=new pedestalSubtraction*[ny]; + /* pedMean=new double*[ny]; */ + /* pedVariance=new double*[ny]; */ for (int i=0; iSetNPedestals(); //cout << nped << " " << orig->getPedestal(ix,iy) << orig->getPedestalRMS(ix,iy) << endl; - for (int iy=0; iygetPedestal(ix,iy),orig->getPedestalRMS(ix,iy),orig->GetNPedestals(ix,iy)); + /* if (ix==50 && iy==50) */ + /* cout << "clone ped " << " " << ix << " " << iy << " " << getPedestal(ix,iy) << " " << getPedestalRMS(ix,iy)<< " " << GetNPedestals(ix,iy) << endl; */ } } image=new int[nx*ny]; @@ -231,8 +252,8 @@ template class analogDetector { if (gm) { if (gmap) delete [] gmap; gmap=new double[nnx*nny]; - for (int iy=0; iy class analogDetector { void *ret; if (gmap) { gm=new float[nx*ny]; - for (int iy=0; iy class analogDetector { virtual void newDataSet(){ iframe=-1; - for (int iy=0; iy class analogDetector { }; /** resets the commonModeSubtraction and increases the frame index */ - virtual void newFrame(){iframe++; if (cmSub) cmSub->newFrame();}; + virtual void newFrame(){iframe++; if (cmSub) cmSub->newFrame(); det->newFrame();}; /** resets the commonModeSubtraction and increases the frame index */ virtual void newFrame(char *data){ iframe++; if (cmSub) cmSub->newFrame(); + det->newFrame(); + // det->getData(data); calcGhost(data); // cout << getId() << " Calc ghost " << getGhost(15,15) << endl; - }; + }; /** sets the commonModeSubtraction algorithm to be used @@ -359,8 +382,8 @@ template class analogDetector { // cout << "+"<< getId() << endl; if (cmSub) { //cout << "*" << endl; - for (int iy=ymin; iy0) // if (det->isGood(ix,iy)) { addToCommonMode(data, ix, iy); @@ -403,8 +426,10 @@ template class analogDetector { if (ix>=0 && ix=0 && iy0) { return stat[iy][ix].getPedestal()+getCommonMode(ix,iy); - - } else return stat[iy][ix].getPedestal(); + //return pedMean[iy][ix]+getCommonMode(ix,iy); + } + //return pedMean[iy][ix]; + return stat[iy][ix].getPedestal(); } else return -1; }; @@ -422,7 +447,7 @@ template class analogDetector { g=gmap[iy*nx+ix]; if (g==0) g=-1.; } - + // return sqrt(pedVariance[iy][ix])/g; return stat[iy][ix].getPedestalRMS()/g;//divide by gain? } return -1; @@ -444,8 +469,8 @@ template class analogDetector { if (ped==NULL) { ped=new double[nx*ny]; } - for (int iy=0; iy class analogDetector { if (ped==NULL) { ped=new double[nx*ny]; } - for (int iy=0; iy class analogDetector { */ virtual void setPedestal(double *ped, double *rms=NULL, int m=-1){ double rr=0; - for (int iy=ymin; iy class analogDetector { \param rms pointer to array of pedestal rms */ virtual void setPedestalRMS(double *rms){ - for (int iy=ymin; iy class analogDetector { #endif gm=new float[nx*ny]; - for (int iy=0; iySetBinContent(ix+1, iy+1,image[iy*nx+ix]); @@ -622,8 +647,8 @@ template class analogDetector { TH2F *hmap=new TH2F("hmap","hmap",nx, -0.5,nx-0.5, ny, -0.5, ny-0.5); #endif - for (int iy=0; iygetCommonMode(); */ /* else */ @@ -669,8 +694,8 @@ template class analogDetector { if (gm) { - for (int iy=0; iy class analogDetector { if (gm) { - for (int iy=0; iy class analogDetector { float *gm=NULL; void *ret; gm=new float[nx*ny]; - for (int iy=0; iy class analogDetector { if (nnx>nx) nnx=nx; if (nny>ny) nny=ny; if (gm) { - for (int iy=0; iy class analogDetector { //cout << xmin << " " << xmax << endl; // cout << ymin << " " << ymax << endl; - for (int iy=ymin; iyisGood(ix,iy)) { // addToPedestal(data,ix,iy,1); addToPedestal(data,ix,iy,cm); + /* if (ix==50 && iy==50) */ + /* cout<< "*ped* " << id << " " << ix << " " << iy << " " << det->getChannel(data,ix,iy) << " " << stat[iy][ix].getPedestal() << " " << stat[iy][ix].getPedestalRMS() <<" " << stat[iy][ix]. GetNPedestals() << endl; */ //if (ix==10 && iy==10) // cout < class analogDetector { //calcGhost(data); - for (int iy=ymin; iyisGood(ix,iy)) val[iy*nx+ix]+=subtractPedestal(data, ix, iy,cm); } @@ -954,8 +981,8 @@ template class analogDetector { hs->Fill(val,(iy-ymin)*(xmax-xmin)+(ix-xmin)); #ifdef ROOTCLUST double v3=0,v5=0,v7=0,v9=0; - for (int iix=-4; iix<5; iix++) - for (int iiy=-4; iiy<5; iiy++) { + for (int iix=-4; iix<5; i++ix) + for (int iiy=-4; iiy<5; i++iy) { if (det) val= (dataSign*det->getValue(data, ix+iix, iy+iiy)-getPedestal(ix+iix,iy+iiy,cm))/g; else @@ -1049,8 +1076,8 @@ template class analogDetector { //calcGhost(data); addToCommonMode(data); - for (int iy=ymin; iyisGood(ix,iy)) nph[iy*nx+ix]+=getNPhotons(data, ix, iy); } @@ -1064,8 +1091,8 @@ template class analogDetector { */ virtual void clearImage(){ - for (int iy=0; iy class analogDetector { int SetNPedestals(int i=-1) { int ix=0, iy=0; if (i>0) - for (iy=0; iy class analogDetector { newFrame(data); //calcGhost(data); addToCommonMode(data); - for (int iy=ymi; iyisGood(ix,iy)) { if (ix>=0 && ix=0 && iydataSize) dsize=dataSize; + for (int ip=0; ip<(el); ip++) { getPixel(ip,ix,iy); if (ix>=0 && ix=0 && iy -//#include +#include #include // time_t #include using namespace std; -//using namespace std::chrono; +using namespace std::chrono; //#define SLS_DETECTOR_JSON_HEADER_VERSION 0x2 @@ -77,9 +77,10 @@ int main(int argc, char *argv[]) { int ok; - // high_resolution_clock::time_point t1; - // high_resolution_clock::time_point t2 ; - time_t begin,end,finished; + high_resolution_clock::time_point t1; + high_resolution_clock::time_point t2 ; + std::chrono::steady_clock::time_point begin,end,finished; + //time_t begin,end,finished; int rms=0; @@ -141,8 +142,12 @@ int main(int argc, char *argv[]) { int ncol_cm=CM_ROWS; double xt_ghost=C_GHOST; - moench03CommonMode *cm=new moench03CommonMode(ncol_cm); - moench03GhostSummation *gs=new moench03GhostSummation(det, xt_ghost); + moench03CommonMode *cm=NULL; + moench03GhostSummation *gs=NULL; +#ifdef CORR + cm=new moench03CommonMode(ncol_cm); + gs=new moench03GhostSummation(det, xt_ghost); +#endif double *gainmap=NULL; float *gm; double *gmap=NULL; @@ -170,7 +175,7 @@ int main(int argc, char *argv[]) { //analogDetector *filter=new analogDetector(det,1,NULL,1000); #ifndef INTERP - singlePhotonDetector *filter=new singlePhotonDetector(det,3, nSigma, 1, cm, 1000, 10, -1, -1, gainmap, gs); + singlePhotonDetector *filter=new singlePhotonDetector(det,3, nSigma, 1, cm, 1000, 100, -1, -1, gainmap, gs); multiThreadedCountingDetector *mt=new multiThreadedCountingDetector(filter,nthreads,fifosize); @@ -356,10 +361,13 @@ int main(int argc, char *argv[]) { #endif // if (!zmqsocket->ReceiveHeader(0, acqIndex, frameIndex, subframeIndex, filename, fileindex)) { - // cprintf(RED, "Got Dummy\n"); + cprintf(RED, "Got Dummy\n"); // t1=high_resolution_clock::now(); - time(&end); + //time(&end); + //cout << "Measurement lasted " << difftime(end,begin) << endl; + end = std::chrono::steady_clock::now(); + cout << "Measurement lasted " << (end-begin).count()*0.000001 << " ms" << endl; while (mt->isBusy()) {;}//wait until all data are processed from the queues @@ -521,14 +529,16 @@ int main(int argc, char *argv[]) { mt->clearImage(); newFrame=1; - //t2 = high_resolution_clock::now(); - time(&finished); - // auto meas_duration = duration_cast( t2 - t0 ).count(); - // auto real_duration = duration_cast( t2 - t1 ).count(); - cout << "Measurement lasted " << difftime(end,begin) << endl; - cout << "Processing lasted " << difftime(finished,begin) << endl; + //time(&finished); + //cout << "Processing lasted " << difftime(finished,begin) << endl; + + finished = std::chrono::steady_clock::now(); + cout << "Processing lasted " << (finished-begin).count()*0.000001 << " ms" << endl; +#ifdef OPTIMIZE + return 0; +#endif continue; //continue to not get out @@ -536,7 +546,8 @@ int main(int argc, char *argv[]) { #ifdef NEWZMQ if (newFrame) { - time(&begin); + begin = std::chrono::steady_clock::now(); + //time(&begin); // t0 = high_resolution_clock::now(); //cout <<"new frame" << endl; @@ -808,6 +819,7 @@ int main(int argc, char *argv[]) { // timestamp=doc["timestamp"].GetUint(); packetNumber=doc["packetNumber"].GetUint(); // cout << acqIndex << " " << frameIndex << " " << subFrameIndex << " "<< bunchId << " " << timestamp << " " << packetNumber << endl; + //cprintf(GREEN, "frame\n"); if (packetNumber>=40) { //*((int*)buff)=frameIndex; if (insubframe==0) f0=frameIndex; diff --git a/slsDetectorCalibration/pedestalSubtraction.h b/slsDetectorCalibration/pedestalSubtraction.h index 062b0d84e..35c010d60 100644 --- a/slsDetectorCalibration/pedestalSubtraction.h +++ b/slsDetectorCalibration/pedestalSubtraction.h @@ -10,7 +10,7 @@ class pedestalSubtraction { \param nn number of samples to calculate the moving average (defaults to 1000) */ pedestalSubtraction (int nn=1000) : stat(nn) {}; - + /* void setPointers(double *me, double *va) {mean=me; var=va; stat.setPointers(mean, var);}; */ /** virtual destructorr */ virtual ~pedestalSubtraction() {}; diff --git a/slsDetectorCalibration/singlePhotonDetector.h b/slsDetectorCalibration/singlePhotonDetector.h index f863bc28a..a6c9fe955 100644 --- a/slsDetectorCalibration/singlePhotonDetector.h +++ b/slsDetectorCalibration/singlePhotonDetector.h @@ -50,7 +50,6 @@ public analogDetector { */ - singlePhotonDetector(slsDetectorData *d, int csize=3, @@ -58,7 +57,7 @@ public analogDetector { int sign=1, commonModeSubtraction *cm=NULL, int nped=1000, - int nd=100, int nnx=-1, int nny=-1, double *gm=NULL, ghostSummation *gs=NULL) : analogDetector(d, sign, cm, nped, nnx, nny, gm, gs), nDark(nd), eventMask(NULL),nSigma (nsigma), eMin(-1), eMax(-1), clusterSize(csize), clusterSizeY(csize), clusters(NULL), quad(UNDEFINED_QUADRANT), tot(0), quadTot(0) { + int nd=100, int nnx=-1, int nny=-1, double *gm=NULL, ghostSummation *gs=NULL) : analogDetector(d, sign, cm, nped, nnx, nny, gm, gs), nDark(nd), eventMask(NULL),nSigma (nsigma), eMin(-1), eMax(-1), clusterSize(csize), clusterSizeY(csize), c2(1),c3(1), clusters(NULL), quad(UNDEFINED_QUADRANT), tot(0), quadTot(0) { @@ -66,13 +65,18 @@ public analogDetector { fm=new pthread_mutex_t ; eventMask=new eventType*[ny]; + // val=new double*[ny]; for (int i=0; i { myFile=orig->myFile; eventMask=new eventType*[ny]; + // val=new double*[ny]; for (int i=0; ieMin; eMax=orig->eMax; @@ -111,6 +117,10 @@ public analogDetector { clusterSize=orig->clusterSize; clusterSizeY=orig->clusterSizeY; // cluster=new single_photon_hit(clusterSize,clusterSizeY); + + c2=sqrt((clusterSizeY+1)/2* (clusterSize+1)/2); + c3=sqrt(clusterSizeY*clusterSize); + clusters=new single_photon_hit[nx*ny]; // cluster=clusters; @@ -124,7 +134,9 @@ public analogDetector { gmap=orig->gmap; nphTot=0; nphFrame=0; - + nphTot=0; + nphFrame=0; + } @@ -196,7 +208,6 @@ public analogDetector { int nn=0; double max=0, tl=0, tr=0, bl=0,br=0, v; double rms=0; - int cm=0; if (cmSub) cm=1; @@ -218,8 +229,8 @@ public analogDetector { cout << "add to common mode?"<< endl; addToCommonMode(data); } - for (int iy=ymin; iyisGood(ix,iy)) { val=subtractPedestal(data,ix,iy, cm); @@ -236,8 +247,8 @@ public analogDetector { } } - for (int iy=ymin; iyisGood(ix,iy)) { eventMask[iy][ix]=PEDESTAL; @@ -346,18 +357,23 @@ int *getClusters(char *data, int *ph=NULL) { int nph=0; - double val[ny][nx]; - int cy=(clusterSizeY+1)/2; - int cs=(clusterSize+1)/2; + // const int cy=(clusterSizeY+1)/2; + //const int cs=(clusterSize+1)/2; + //int ir, ic; - + eventType ee; double max=0, tl=0, tr=0, bl=0,br=0, *v; int cm=0; int good=1; - if (cmSub) cm=1; + int ir, ic; + // double quadTot; + //quadrant quad; + double rms; + //if (cmSub) cm=1; + double val[ny][nx]; if (ph==NULL) ph=image; - + if (iframeisGood(ix,iy)) { + for (iy=ymin; iyisGood(ix,iy)==0) continue; + max=0; tl=0; tr=0; @@ -384,97 +402,120 @@ int *getClusters(char *data, int *ph=NULL) { - eventMask[iy][ix]=PEDESTAL; + //eventMask[iy][ix] + ee=PEDESTAL; - (clusters+nph)->rms=getPedestalRMS(ix,iy); + rms=getPedestalRMS(ix,iy); // cluster=clusters+nph; - - for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) { - for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) { + // cout << ix << " " << iy << endl; + for (ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) { + for (ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) { if ((iy+ir)>=iy && (iy+ir)=ix && (ix+ic)=0) + br+=*v; + if (ir>=0 && ic<=0) + tl+=*v; + if (ir>=0 && ic>=0) + tr+=*v; + if (*v>max) { + max=*v; + } + + } - v=&(val[iy+ir][ix+ic]); - tot+=*v; - if (ir<=0 && ic<=0) - bl+=*v; - if (ir<=0 && ic>=0) - br+=*v; - if (ir>=0 && ic<=0) - tl+=*v; - if (ir>=0 && ic>=0) - tr+=*v; - if (*v>max) { - max=*v; - } - - - if (ir==0 && ic==0) { - if (*v<-nSigma*(clusters+nph)->rms) - eventMask[iy][ix]=NEGATIVE_PEDESTAL; - else if (*v>nSigma*(clusters+nph)->rms) - eventMask[iy][ix]=PHOTON; - } } } - if (eventMask[iy][ix]==PHOTON && val[iy][ix]=br && bl>=tl && bl>=tr) { - (clusters+nph)->quad=BOTTOM_LEFT; - (clusters+nph)->quadTot=bl; - } else if (br>=bl && br>=tl && br>=tr) { - (clusters+nph)->quad=BOTTOM_RIGHT; - (clusters+nph)->quadTot=br; - } else if (tl>=br && tl>=bl && tl>=tr) { - (clusters+nph)->quad=TOP_LEFT; - (clusters+nph)->quadTot=tl; - } else if (tr>=bl && tr>=tl && tr>=br) { - (clusters+nph)->quad=TOP_RIGHT; - (clusters+nph)->quadTot=tr; + /* if (ix==50 && iy==50) */ + /* cout << id << " " << ix << " " << iy << " " << det->getValue(data,ix,iy)<< " " << val[iy][ix] << " " << getPedestal(ix,iy) << " " << rms << endl; */ + if (val[iy][ix]<-nSigma*rms){ + ee=NEGATIVE_PEDESTAL; + continue; + } + if (max>nSigma*rms){ + // cout << "ph1 " << max << " " << nSigma*rms << endl; + ee=PHOTON; + if (val[iy][ix]nSigma*(clusters+nph)->rms || tot>sqrt(clusterSizeY*clusterSize)*nSigma*(clusters+nph)->rms || ((clusters+nph)->quadTot)>sqrt(cy*cs)*nSigma*(clusters+nph)->rms) { - if (val[iy][ix]>=max) { - eventMask[iy][ix]=PHOTON_MAX; - (clusters+nph)->tot=tot; - (clusters+nph)->x=ix; - (clusters+nph)->y=iy; - // (clusters+nph)->iframe=det->getFrameNumber(data); - // cout << det->getFrameNumber(data) << " " << (clusters+nph)->iframe << endl; - (clusters+nph)->ped=getPedestal(ix,iy,0); - for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) { - for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) { - (clusters+nph)->set_data(val[iy+ir][ix+ic],ic,ir); + else if (tot>c3*nSigma*rms) { + // cout << "ph3 " << tot << " " << c3*nSigma*rms << endl; + ee=PHOTON; + } +#ifndef WRITE_QUAD + else { +#endif + quad=BOTTOM_RIGHT; + quadTot=br; + if (bl>=quadTot) { + quad=BOTTOM_LEFT; + quadTot=bl; + } + if (tl>=quadTot) { + quad=TOP_LEFT; + quadTot=tl; + } + if (tr>=quadTot) { + quad=TOP_RIGHT; + quadTot=tr; + } + + if (quadTot>c2*nSigma*rms) { + // cout << "ph2 " << quadTot << " " << c2*nSigma*rms << endl; + ee=PHOTON; } - } - good=1; - if (eMin>0 && tot0 && tot>eMax) good=0; - if (good) { - nph++; - image[iy*nx+ix]++; - } - - } else { - eventMask[iy][ix]=PHOTON; - } - } else if (eventMask[iy][ix]==PEDESTAL) { - addToPedestal(data,ix,iy,cm); +#ifndef WRITE_QUAD } +#endif + if (ee==PHOTON && val[iy][ix]==max) { + ee=PHOTON_MAX; + // cout << "**" <tot=tot; + (clusters+nph)->x=ix; + (clusters+nph)->y=iy; + (clusters+nph)->quad=quad; + (clusters+nph)->quadTot=quadTot; + //(clusters+nph)->rms=rms; + // (clusters+nph)->iframe=det->getFrameNumber(data); + // cout << det->getFrameNumber(data) << " " << (clusters+nph)->iframe << endl; + // (clusters+nph)->ped=getPedestal(ix,iy,0); + for (ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) { + for (ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) { + if ((iy+ir)>=iy && (iy+ir)=ix && (ix+ic)set_data(val[iy+ir][ix+ic],ic,ir); + } } - } + good=1; + if (eMin>0 && tot0 && tot>eMax) good=0; + if (good) { + nph++; + image[iy*nx+ix]++; + } + + + } else if (ee==PEDESTAL) { + addToPedestal(data,ix,iy,cm); + } /*else { + eventMask[iy][ix]=PHOTON; + }*/ + //eventMask[iy][ix]=ee; } +} + nphFrame=nph; nphTot+=nph; //cout << nphFrame << endl; - // cout <<"**********************************"<< det->getFrameNumber(data) << " " << nphFrame << endl; + //cout <getFrameNumber(data) << " " << nphFrame << endl; writeClusters(det->getFrameNumber(data)); return image; @@ -561,13 +602,20 @@ int *getClusters(char *data, int *ph=NULL) { */ static void writeClusters(FILE *f, single_photon_hit *cl, int nph, int fn=0){ - -/* #ifndef OLDFORMAT */ -/* if (fwrite((void*)&fn, 1, sizeof(int), f)) */ -/* if (fwrite((void*)&nph, 1, sizeof(int), f)) */ -/* #endif */ - for (int i=0; iwrite(f); -}; + if (nph>0) { +#ifndef OLDFORMAT + if (fwrite((void*)&fn, 1, sizeof(int), f)) + if (fwrite((void*)&nph, 1, sizeof(int), f)) +#endif + for (int i=0; iwrite(f); + } + }; + + + + + + void writeClusters(FILE *f, int fn=0){ writeClusters(f,clusters,nphFrame, fn); //for (int i=0; i