From 19e7ced3325a62c585d41eba6750afc268eb3e84 Mon Sep 17 00:00:00 2001 From: Anna Bergamaschi Date: Thu, 13 Sep 2018 16:07:43 +0200 Subject: [PATCH] moench zmq process implemented - command implementation is missing --- slsDetectorCalibration/analogDetector.h | 30 +- .../interpolatingDetector.h | 219 +++----- .../interpolations/eta2InterpolationBase.h | 24 +- .../interpolations/eta3InterpolationBase.h | 20 +- .../interpolations/etaInterpolationBase.h | 36 +- .../interpolations/etaInterpolationPosXY.h | 4 + .../interpolations/slsInterpolation.h | 20 +- .../moenchExecutables/moenchZmqProcess.cpp | 369 ++++++++----- .../multiThreadedAnalogDetector.h | 209 ++------ .../multiThreadedCountingDetector.h | 55 ++ .../multiThreadedInterpolatingDetector.h | 136 +++++ slsDetectorCalibration/singlePhotonDetector.h | 483 +++++++++--------- 12 files changed, 867 insertions(+), 738 deletions(-) create mode 100644 slsDetectorCalibration/multiThreadedCountingDetector.h create mode 100644 slsDetectorCalibration/multiThreadedInterpolatingDetector.h diff --git a/slsDetectorCalibration/analogDetector.h b/slsDetectorCalibration/analogDetector.h index 5df0b5b1f..029808fd8 100644 --- a/slsDetectorCalibration/analogDetector.h +++ b/slsDetectorCalibration/analogDetector.h @@ -80,7 +80,6 @@ template class analogDetector { fMode=ePedestal; thr=0; myFile=NULL; - fm=new pthread_mutex_t ; #ifdef ROOTSPECTRUM hs=new TH2F("hs","hs",2000,-100,10000,nx*ny,-0.5,nx*ny-0.5); #ifdef ROOTCLUST @@ -128,7 +127,6 @@ template class analogDetector { // nSigma=orig->nSigma; fMode=orig->fMode; myFile=orig->myFile; - fm=orig->fm; stat=new pedestalSubtraction*[ny]; @@ -329,7 +327,8 @@ template class analogDetector { for (int ix=xmin; ix0) - addToCommonMode(data, ix, iy); + if (det->isGood(ix,iy)) + addToCommonMode(data, ix, iy); } } //cout << "cm " << getCommonMode(0,0) << " " << getCommonMode(1,0) << endl; @@ -701,12 +700,14 @@ template class analogDetector { // cout << ymin << " " << ymax << endl; for (int ix=xmin; ixisGood(ix,iy)) { + addToPedestal(data,ix,iy,1); //if (ix==10 && iy==10) // cout < class analogDetector { for (int ix=xmin; ixisGood(ix,iy)) + val[iy*nx+ix]+=subtractPedestal(data, ix, iy,cm); } } return val; @@ -950,7 +952,8 @@ template class analogDetector { for (int ix=xmin; ixisGood(ix,iy)) + nph[iy*nx+ix]+=getNPhotons(data, ix, iy); } } return nph; @@ -1028,8 +1031,11 @@ template class analogDetector { for (int ix=xmi; ix=0 && ix=0 && iyisGood(ix,iy)) { + if (ix>=0 && ix=0 && iy class analogDetector { } }; - virtual char *getInterpolation(){return NULL;}; + // virtual char *getInterpolation(){return NULL;}; /** sets the current frame mode for the detector \param f frame mode to be set @@ -1095,7 +1101,6 @@ FILE *setFilePointer(FILE *f){myFile=f; return myFile;}; \returns current file pointer */ FILE *getFilePointer(){return myFile;}; - void setMutex(pthread_mutex_t *m){fm=m;}; protected: slsDetectorData *det; /**< slsDetectorData to be used */ @@ -1127,7 +1132,6 @@ FILE *getFilePointer(){return myFile;}; TH2F *hs9; #endif #endif - pthread_mutex_t *fm; }; #endif diff --git a/slsDetectorCalibration/interpolatingDetector.h b/slsDetectorCalibration/interpolatingDetector.h index 738427a87..b99e62950 100644 --- a/slsDetectorCalibration/interpolatingDetector.h +++ b/slsDetectorCalibration/interpolatingDetector.h @@ -49,17 +49,21 @@ class interpolatingDetector : public singlePhotonDetector { int nd=100, int nnx=-1, int nny=-1) : singlePhotonDetector(d, 3,nsigma,sign, cm, nped, nd, nnx, nny) , interp(inte), id(0) { //cout << "**"<< xmin << " " << xmax << " " << ymin << " " << ymax << endl; + fi=new pthread_mutex_t ; }; interpolatingDetector(interpolatingDetector *orig) : singlePhotonDetector(orig) { - if (orig->interp) - interp=(orig->interp)->Clone(); - else - interp=orig->interp; + // if (orig->interp) + // interp=(orig->interp)->Clone(); + // else + + interp=orig->interp; + id=orig->id; + fi=orig->fi; } @@ -67,7 +71,11 @@ class interpolatingDetector : public singlePhotonDetector { return new interpolatingDetector(this); } - virtual int setId(int i) {id=i; interp->setId(id); return id;}; + virtual int setId(int i) { + id=i; + // interp->setId(id); + return id; + }; virtual void prepareInterpolation(int &ok) { /* cout << "*"<< endl; */ @@ -81,20 +89,28 @@ class interpolatingDetector : public singlePhotonDetector { /* sprintf(tit,"/scratch/gmap_%d.tiff",id); */ /* writeGainMap(tit); */ /* } */ -/* #endif */ +/* #endif */ + pthread_mutex_lock(fi); if (interp) interp->prepareInterpolation(ok); + pthread_mutex_unlock(fi); } void clearImage() { - if (interp) interp->clearInterpolatedImage(); - else singlePhotonDetector::clearImage(); + if (interp) { + pthread_mutex_lock(fi); + interp->clearInterpolatedImage(); + pthread_mutex_unlock(fi); + } else + singlePhotonDetector::clearImage(); }; int getImageSize(int &nnx, int &nny, int &ns) { - if (interp) return interp->getImageSize(nnx, nny, ns); - else return analogDetector::getImageSize(nnx, nny, ns); + if (interp) + return interp->getImageSize(nnx, nny, ns); + else + return analogDetector::getImageSize(nnx, nny, ns); }; @@ -152,128 +168,7 @@ class interpolatingDetector : public singlePhotonDetector { return NULL; } -/* int addFrame(char *data, int *ph=NULL, int ff=0) { */ - -/* int nph=0; */ -/* double val[ny][nx]; */ -/* int cy=(clusterSizeY+1)/2; */ -/* int cs=(clusterSize+1)/2; */ -/* int ir, ic; */ -/* double rms; */ -/* double int_x,int_y, eta_x, eta_y; */ -/* double max=0, tl=0, tr=0, bl=0,br=0, *v, vv; */ -/* // cout << "********** Add frame "<< iframe << endl; */ -/* if (ph==NULL) */ -/* ph=image; */ - -/* if (iframerms=rms; */ - - - -/* for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) { */ -/* for (int 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; */ -/* } */ - - -/* if (ir==0 && ic==0) { */ -/* if (*v<-nSigma*rms) */ -/* eventMask[iy][ix]=NEGATIVE_PEDESTAL; */ -/* } */ - -/* } */ -/* } */ - -/* if (bl>=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 (max>nSigma*rms || tot>sqrt(clusterSizeY*clusterSize)*nSigma*rms || ((clusters+nph)->quadTot)>sqrt(cy*cs)*nSigma*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); */ -/* (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); */ -/* } */ -/* } */ - - - -/* nph++; */ -/* image[iy*nx+ix]++; */ - -/* } else { */ -/* eventMask[iy][ix]=PHOTON; */ -/* } */ -/* } else if (eventMask[iy][ix]==PEDESTAL) { */ -/* addToPedestal(data,ix,iy); */ -/* } */ - - -/* } */ -/* } */ -/* nphFrame=nph; */ -/* nphTot+=nph; */ -/* writeClusters(); */ - -/* return nphFrame; */ - -/* }; */ int addFrame(char *data, int *ph=NULL, int ff=0) { @@ -283,6 +178,7 @@ int addFrame(char *data, int *ph=NULL, int ff=0) { double int_x, int_y; double eta_x, eta_y; if (interp) { + pthread_mutex_lock(fi); for (nph=0; nphaddToFlatField((clusters+nph)->quadTot,(clusters+nph)->quad,(clusters+nph)->get_cluster(),eta_x, eta_y); @@ -290,38 +186,75 @@ int addFrame(char *data, int *ph=NULL, int ff=0) { interp->getInterpolatedPosition((clusters+nph)->x, (clusters+nph)->y, (clusters+nph)->quadTot,(clusters+nph)->quad,(clusters+nph)->get_cluster(),int_x, int_y); interp->addToImage(int_x, int_y); } - } + } + pthread_mutex_lock(fi); } return nphFrame; }; virtual void processData(char *data, int *val=NULL) { - if (interp){ - switch(fMode) { - case ePedestal: - addToPedestal(data); - break; - case eFlat: - addFrame(data,val,1); + switch (dMode) { + case eAnalog: + analogDetector::processData(data,val); break; + case ePhotonCounting: + singlePhotonDetector::processData(data,val); default: - addFrame(data,val,0); + switch(fMode) { + case ePedestal: + addToPedestal(data); + break; + case eFlat: + if (interp) + addFrame(data,val,1); + else + singlePhotonDetector::processData(data,val); + break; + default: + if (interp) + addFrame(data,val,0); + else + singlePhotonDetector::processData(data,val); + } } - } else - singlePhotonDetector::processData(data,val); + }; virtual char *getInterpolation(){return (char*)interp;}; - virtual char *setInterpolation(char *ii){int ok; interp=(slsInterpolation*)ii; interp->prepareInterpolation(ok); return (char*)interp;}; + virtual char *setInterpolation(char *ii){ + int ok; + interp=(slsInterpolation*)ii; + pthread_mutex_lock(fi); + interp->prepareInterpolation(ok); + pthread_mutex_unlock(fi); + return (char*)interp;}; + virtual void resetFlatField() { if (interp) { + pthread_mutex_lock(fi); + interp->resetFlatField(); + pthread_mutex_unlock(fi); + } + } + + virtual int getNSubPixels(){ if (interp) return interp->getNSubPixels(); else return 1;} + virtual int setNSubPixels(int ns) { + if (interp) { + pthread_mutex_lock(fi); + interp->getNSubPixels(); + pthread_mutex_unlock(fi); + } + return getNSubPixels(); + } + protected: slsInterpolation *interp; int id; + pthread_mutex_t *fi; }; diff --git a/slsDetectorCalibration/interpolations/eta2InterpolationBase.h b/slsDetectorCalibration/interpolations/eta2InterpolationBase.h index 55ed32e2a..6bfade227 100644 --- a/slsDetectorCalibration/interpolations/eta2InterpolationBase.h +++ b/slsDetectorCalibration/interpolations/eta2InterpolationBase.h @@ -16,13 +16,13 @@ class eta2InterpolationBase : public virtual etaInterpolationBase { eta2InterpolationBase(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny, ns, nb, emin, emax) { // cout << "e2ib " << nb << " " << emin << " " << emax << endl; - if (nbeta<=0) { - nbeta=nSubPixels*10; - } + /* if (nbeta<=0) { */ + /* nbeta=nSubPixels*10; */ + /* } */ if (etamin>=etamax) { etamin=-1; etamax=2; - cout << ":" <=etamax) { etamin=-1; etamax=1; @@ -33,13 +33,13 @@ class eta3InterpolationBase : public virtual etaInterpolationBase { hhy=new TH2D("hhy","hhy",nbeta,etamin,etamax,nbeta,etamin,etamax); #endif #ifndef MYROOT1 - delete [] heta; - delete [] hhx; - delete [] hhy; + /* delete [] heta; */ + /* delete [] hhx; */ + /* delete [] hhy; */ - heta=new int[nbeta*nbeta]; - hhx=new float[nbeta*nbeta]; - hhy=new float[nbeta*nbeta]; + /* heta=new int[nbeta*nbeta]; */ + /* hhx=new float[nbeta*nbeta]; */ + /* hhy=new float[nbeta*nbeta]; */ #endif // cout << nbeta << " " << etamin << " " << etamax << endl; @@ -47,7 +47,7 @@ class eta3InterpolationBase : public virtual etaInterpolationBase { eta3InterpolationBase(eta3InterpolationBase *orig): etaInterpolationBase(orig){ }; - virtual eta3InterpolationBase* Clone()=0; + /* virtual eta3InterpolationBase* Clone()=0; */ diff --git a/slsDetectorCalibration/interpolations/etaInterpolationBase.h b/slsDetectorCalibration/interpolations/etaInterpolationBase.h index 1218320f9..9b50b8e46 100644 --- a/slsDetectorCalibration/interpolations/etaInterpolationBase.h +++ b/slsDetectorCalibration/interpolations/etaInterpolationBase.h @@ -19,11 +19,11 @@ class etaInterpolationBase : public slsInterpolation { // cout << "eb " << nb << " " << emin << " " << emax << endl; // cout << nb << " " << etamin << " " << etamax << endl; if (nbeta<=0) { - cout << "aaa:" <=etamax) { - cout << "aaa:" <Reset(); + hhx->Reset(); + hhy->Reset(); +#endif + }; - */ - - + #ifdef MYROOT1 TH2D *setEta(TH2D *h, int nb=-1, double emin=1, double emax=0) @@ -117,6 +136,9 @@ class etaInterpolationBase : public slsInterpolation { { return setEta(h, nb, emin, emax); }; + + + int *getFlatField(){return setEta(NULL);}; int *getFlatField(int &nb, double &emin, double &emax){ @@ -183,6 +205,8 @@ class etaInterpolationBase : public slsInterpolation { }; + + #endif diff --git a/slsDetectorCalibration/interpolations/etaInterpolationPosXY.h b/slsDetectorCalibration/interpolations/etaInterpolationPosXY.h index 7730771cd..49b8ea044 100644 --- a/slsDetectorCalibration/interpolations/etaInterpolationPosXY.h +++ b/slsDetectorCalibration/interpolations/etaInterpolationPosXY.h @@ -154,11 +154,15 @@ class eta2InterpolationPosXY : public virtual eta2InterpolationBase, public virt eta2InterpolationPosXY(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny,ns, nb, emin,emax),eta2InterpolationBase(nx,ny,ns, nb, emin,emax),etaInterpolationPosXY(nx,ny,ns, nb, emin,emax){ // cout << "e2pxy " << nb << " " << emin << " " << emax << endl; }; + eta2InterpolationPosXY(eta2InterpolationPosXY *orig): etaInterpolationBase(orig), etaInterpolationPosXY(orig) {}; virtual eta2InterpolationPosXY* Clone() { return new eta2InterpolationPosXY(this);}; }; + + + class eta3InterpolationPosXY : public virtual eta3InterpolationBase, public virtual etaInterpolationPosXY { public: eta3InterpolationPosXY(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny,ns, nb, emin,emax),eta3InterpolationBase(nx,ny,ns, nb, emin,emax), etaInterpolationPosXY(nx,ny,ns, nb, emin,emax){ diff --git a/slsDetectorCalibration/interpolations/slsInterpolation.h b/slsDetectorCalibration/interpolations/slsInterpolation.h index 8059ddcc9..ad095488d 100644 --- a/slsDetectorCalibration/interpolations/slsInterpolation.h +++ b/slsDetectorCalibration/interpolations/slsInterpolation.h @@ -56,18 +56,30 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny); #endif #ifndef MYROOT1 - hint=new int[nSubPixels*nPixelsX*nSubPixels*nPixelsY]; - memcpy(hint, orig->hint,nSubPixels*nPixelsX*nSubPixels*nPixelsY*sizeof(int)); + hint=new int[nSubPixels*nPixelsX*nSubPixels*nPixelsY]; + memcpy(hint, orig->hint,nSubPixels*nPixelsX*nSubPixels*nPixelsY*sizeof(int)); #endif }; virtual int setId(int i) {id=i; return id;}; - virtual slsInterpolation* Clone() = 0; + virtual slsInterpolation* Clone() =0; /*{ + return new slsInterpolation(this); + }*/ int getNSubPixels() {return nSubPixels;}; + + int setNSubPixels(int ns) { + if (ns>0 && ns!=nSubPixels) { + delete [] hint; + nSubPixels=ns; + hint=new int[nSubPixels*nPixelsX*nSubPixels*nPixelsY]; + } + return nSubPixels; + } + int getImageSize(int &nnx, int &nny, int &ns) { nnx=nSubPixels*nPixelsX; nny=nSubPixels*nPixelsY; @@ -187,6 +199,8 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny); virtual int *getFlatField(int &nb, double &emin, double &emax){nb=0; emin=0; emax=0; return getFlatField();}; #endif + virtual void resetFlatField()=0; + //virtual void Streamer(TBuffer &b); static int calcQuad(int *cl, double &sum, double &totquad, double sDum[2][2]){ diff --git a/slsDetectorCalibration/moenchExecutables/moenchZmqProcess.cpp b/slsDetectorCalibration/moenchExecutables/moenchZmqProcess.cpp index 30d3c795c..d9cd2c1e9 100644 --- a/slsDetectorCalibration/moenchExecutables/moenchZmqProcess.cpp +++ b/slsDetectorCalibration/moenchExecutables/moenchZmqProcess.cpp @@ -17,9 +17,12 @@ #include //#include "analogDetector.h" +//#include "multiThreadedAnalogDetector.h" //#include "singlePhotonDetector.h" -#include "interpolatingDetector.h" -#include "multiThreadedAnalogDetector.h" +//#include "interpolatingDetector.h" +//#include "multiThreadedCountingDetector.h" +#include "multiThreadedInterpolatingDetector.h" +#include "etaInterpolationPosXY.h" #include "ansi.h" #include using namespace std; @@ -36,27 +39,34 @@ int main(int argc, char *argv[]) { */ FILE *of=NULL; int fifosize=1000; - int nthreads=20; + int etabins=1000;//nsubpix*2*100; + double etamin=-1, etamax=2; // help if (argc < 3 ) { - cprintf(RED, "Help: ./trial [receive socket ip] [receive starting port number] [send_socket ip] [send starting port number]\n"); - return EXIT_FAILURE; + cprintf(RED, "Help: ./trial [receive socket ip] [receive starting port number] [send_socket ip] [send starting port number] [nthreads] [nsubpix] [etafile]\n"); + return EXIT_FAILURE; } // receive parameters bool send = false; char* socketip=argv[1]; uint32_t portnum = atoi(argv[2]); - int maxSize = 32*2*8192;//5000;//atoi(argv[3]); - int size= 32*2*5000; - int multisize=size; // send parameters if any char* socketip2 = 0; uint32_t portnum2 = 0; - if (argc > 3) { - send = true; + + + int ok; + + + + + + if (argc > 4) { socketip2 = argv[3]; portnum2 = atoi(argv[4]); + if (portnum2>0) + send = true; } cout << "\nrx socket ip : " << socketip << "\nrx port num : " << portnum ; @@ -64,26 +74,56 @@ int main(int argc, char *argv[]) { cout << "\ntx socket ip : " << socketip2 << "\ntx port num : " << portnum2; } - cout << endl; + int nthreads=5; + if (argc>5) + nthreads=atoi(argv[5]); + + cout << "Number of threads is: " << nthreads << endl; + int nSubPixels=2; + if (argc>6) + nSubPixels=atoi(argv[6]); + cout << "Number of subpixels is: " << nSubPixels << endl; + char *etafname=NULL; + if (argc>7) { + etafname=argv[7]; + cout << "Eta file name is: " << etafname << endl; + } + + //slsDetectorData *det=new moench03T1ZmqDataNew(); moench03T1ZmqDataNew *det=new moench03T1ZmqDataNew(); cout << endl << " det" <getDetectorSize(npx, npy); + + + + + int maxSize = npx*npy*2;//32*2*8192;//5000;//atoi(argv[3]); + int size= maxSize;//32*2*5000; + int multisize=size; + int dataSize=size; + + + + //analogDetector *filter=new analogDetector(det,1,NULL,1000); - singlePhotonDetector *filter=new singlePhotonDetector(det,3, 5, 1, 0, 1000, 10); - // interpolatingDetector *filter=new interpolatingDetector(det,NULL, 5, 1, 0, 1000, 10); + //singlePhotonDetector *filter=new singlePhotonDetector(det,3, 5, 1, 0, 1000, 10); + eta2InterpolationPosXY *interp=new eta2InterpolationPosXY(npx, npy, nSubPixels, etabins, etamin, etamax); + + if (etafname) interp->readFlatField(etafname); + + interpolatingDetector *filter=new interpolatingDetector(det,interp, 5, 1, 0, 1000, 10); cout << " filter" <setFrameMode(eFrame); mt->StartThreads(); - cout << " start" <popFree(buff); - cout << " pop" <IsError()) { - cprintf(RED, "Error: Could not create Zmq socket server on port %d and ip %s\n", portnum2, socketip2); - delete zmqsocket2; - delete zmqsocket; - return EXIT_FAILURE; + cprintf(RED, "AAA Error: Could not create Zmq socket server on port %d and ip %s\n", portnum2, socketip2); + // delete zmqsocket2; + //delete zmqsocket; + // return EXIT_FAILURE; + send = false; } #endif if (zmqsocket2->Connect()) { - cprintf(RED, "Error: Could not connect to socket %s\n", + cprintf(RED, "BBB Error: Could not connect to socket %s\n", zmqsocket2->GetZmqServerAddress()); - delete zmqsocket2; - return EXIT_FAILURE; + // delete zmqsocket2; + send = false; + // return EXIT_FAILURE; } else printf("Zmq Client at %s\n", zmqsocket2->GetZmqServerAddress()); @@ -176,7 +220,6 @@ int main(int argc, char *argv[]) { int *detimage; int nnx, nny,nns; uint32_t imageSize = 0, nPixelsX = 0, nPixelsY = 0, dynamicRange = 0; - filter->getImageSize(nnx, nny,nns); // infinite loop uint32_t packetNumber = 0; uint64_t bunchId = 0; @@ -187,7 +230,7 @@ int main(int argc, char *argv[]) { uint16_t zCoord = 0; uint32_t debug = 0; uint32_t dr = 16; - int16_t *dout=new int16_t [nnx*nny]; + int16_t *dout;//=new int16_t [nnx*nny]; //uint32_t dr = 32; //int32_t *dout=new int32_t [nnx*nny]; uint32_t nSigma=5; @@ -214,6 +257,13 @@ int main(int argc, char *argv[]) { frameMode fMode; double *ped; + filter->getImageSize(nnx, nny,nns); + + + + + + while(1) { @@ -237,48 +287,76 @@ int main(int argc, char *argv[]) { while (mt->isBusy()) {;}//wait until all data are processed from the queues - detimage=mt->getImage(nnx,nny,nns); - if (fMode==ePedestal) { - nns=1; - // cout << "get pedestal " << endl; - ped=mt->getPedestal(); - // cout << "got pedestal " << endl; - for (int ix=0; ixwritePedestal(ofname); + } else if (fMode==eFlat) { + mt->prepareInterpolation(ok); + sprintf(ofname,"eta_%s_%d.tiff",fname,fileindex); + mt->writeFlatField(ofname); + } else { + sprintf(ofname,"%s_%d.tiff",fname,fileindex); + mt->writeImage(ofname); } // cout << nns*nnx*nny*nns*dr/8 << " " << length << endl; if (send) { - strcpy(fname,filename.c_str()); -#ifdef NEWZMQ - //zmqsocket2->SendHeaderData (0, false, SLS_DETECTOR_JSON_HEADER_VERSION, dynamicRange, fileindex, - // nnx, nny, nns*dynamicRange/8,acqIndex, frameIndex, fname, acqIndex, subFrameIndex, packetNumber,bunchId, timestamp, modId, xCoord, yCoord, zCoord,debug, roundRNumber, detType, version, flippedData, additionalJsonHeader); + + if (fMode==ePedestal) { + nns=1; + nnx=npx; + nny=npy; + dout= new int16_t[nnx*nny*nns*nns]; + // cout << "get pedestal " << endl; + ped=mt->getPedestal(); + // cout << "got pedestal " << endl; + for (int ix=0; ixSendHeaderData (0, false, SLS_DETECTOR_JSON_HEADER_VERSION, dr, fileindex, nnx, nny, nns*dr/8,acqIndex, frameIndex, fname, acqIndex, subFrameIndex, packetNumber,bunchId, timestamp, modId, xCoord, yCoord, zCoord,debug, roundRNumber, detType, version, flippedData, additionalJsonHeader); + dout[ix]=ped[ix]; + } + + } else if (fMode==eFlat) { + int nb; + double emi, ema; + int *ff=mt->getFlatField(nb, emi, ema); + nnx=nb; + nny=nb; + dout= new int16_t[nb*nb]; + for (int ix=0; ixgetImage(nnx,nny,nns); + // nns=1; + // nnx=npx; + // nny=npy; + nnx=nnx*nns; + nny=nny*nns; + dout= new int16_t[nnx*nny]; + for (int ix=0; ixSendHeaderData (0, false, SLS_DETECTOR_JSON_HEADER_VERSION, dr, fileindex, nnx, nny, nnx*nny*dr/8,acqIndex, frameIndex, fname, acqIndex, subFrameIndex, packetNumber,bunchId, timestamp, modId, xCoord, yCoord, zCoord,debug, roundRNumber, detType, version, flippedData, additionalJsonHeader); @@ -289,27 +367,27 @@ int main(int argc, char *argv[]) { #endif zmqsocket2->SendData((char*)dout,length);//nns*dr/8); - cprintf(GREEN, "Sent Data\n"); - } - - - // stream dummy to socket2 to signal end of acquisition - if (send) { - zmqsocket2->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION); - cprintf(RED, "Sent Dummy\n"); + cprintf(GREEN, "Sent Data\n"); + + zmqsocket2->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION); + cprintf(RED, "Sent Dummy\n"); + delete [] dout; + } + + mt->clearImage(); - if (of) { - fclose(of); - of=NULL; - } newFrame=1; continue; //continue to not get out + + } #ifdef NEWZMQ if (newFrame) { + + // acqIndex, frameIndex, subframeIndex, filename, fileindex size = doc["size"].GetUint(); multisize = size;// * zmqsocket->size(); @@ -334,36 +412,38 @@ int main(int argc, char *argv[]) { detType=doc["detType"].GetUint(); version=doc["version"].GetUint(); + dataSize=size; + strcpy(fname,filename.c_str()); - cprintf(BLUE, "Header Info:\n" - "size: %u\n" - "multisize: %u\n" - "dynamicRange: %u\n" - "nPixelsX: %u\n" - "nPixelsY: %u\n" - "currentFileName: %s\n" - "currentAcquisitionIndex: %lu\n" - "currentFrameIndex: %lu\n" - "currentFileIndex: %lu\n" - "currentSubFrameIndex: %u\n" - "xCoordX: %u\n" - "yCoordY: %u\n" - "zCoordZ: %u\n" - "flippedDataX: %u\n" - "packetNumber: %u\n" - "bunchId: %u\n" - "timestamp: %u\n" - "modId: %u\n" - "debug: %u\n" - "roundRNumber: %u\n" - "detType: %u\n" - "version: %u\n", - size, multisize, dynamicRange, nPixelsX, nPixelsY, - filename.c_str(), acqIndex, - frameIndex, fileindex, subFrameIndex, - xCoord, yCoord,zCoord, - flippedDataX, packetNumber, bunchId, timestamp, modId, debug, roundRNumber, detType, version); + // cprintf(BLUE, "Header Info:\n" + // "size: %u\n" + // "multisize: %u\n" + // "dynamicRange: %u\n" + // "nPixelsX: %u\n" + // "nPixelsY: %u\n" + // "currentFileName: %s\n" + // "currentAcquisitionIndex: %lu\n" + // "currentFrameIndex: %lu\n" + // "currentFileIndex: %lu\n" + // "currentSubFrameIndex: %u\n" + // "xCoordX: %u\n" + // "yCoordY: %u\n" + // "zCoordZ: %u\n" + // "flippedDataX: %u\n" + // "packetNumber: %u\n" + // "bunchId: %u\n" + // "timestamp: %u\n" + // "modId: %u\n" + // "debug: %u\n" + // "roundRNumber: %u\n" + // "detType: %u\n" + // "version: %u\n", + // size, multisize, dynamicRange, nPixelsX, nPixelsY, + // filename.c_str(), acqIndex, + // frameIndex, fileindex, subFrameIndex, + // xCoord, yCoord,zCoord, + // flippedDataX, packetNumber, bunchId, timestamp, modId, debug, roundRNumber, detType, version); /* Analog detector commands */ isPedestal=0; @@ -386,7 +466,7 @@ int main(int argc, char *argv[]) { fMode=eFlat; isFlat=1; } else if (frameMode_s == "newFlatfield") { - //mt->resetFlatfield(); + mt->resetFlatField(); isFlat=1; //cprintf(MAGENTA, "Resetting flatfield\n"); fMode=eFlat; @@ -397,16 +477,15 @@ int main(int argc, char *argv[]) { cprintf(MAGENTA, "%s\n" , frameMode_s.c_str()); mt->setFrameMode(fMode); - threshold=0; + // threshold=0; cprintf(MAGENTA, "Threshold: "); if (doc.HasMember("threshold")) { - if (doc["threshold"].IsInt()) + if (doc["threshold"].IsInt()) { threshold=doc["threshold"].GetInt(); - } - mt->setThreshold(threshold); + mt->setThreshold(threshold); + } + } cprintf(MAGENTA, "%d\n", threshold); - - xmin=0; xmax=npx; @@ -433,12 +512,13 @@ int main(int argc, char *argv[]) { } } - cprintf(MAGENTA, "%d %d %d %d\n", xmin, xmax, ymin, ymax); + cprintf(MAGENTA, "ROI: %d %d %d %d\n", xmin, xmax, ymin, ymax); mt->setROI(xmin, xmax, ymin, ymax); - // if (doc.HasMember("dynamicRange")) { - // dr=doc["dynamicRange"].GetUint(); - // } + if (doc.HasMember("dynamicRange")) { + dr=doc["dynamicRange"].GetUint(); + dr=16; + } dMode=eAnalog; detectorMode_s="analog"; @@ -463,28 +543,44 @@ int main(int argc, char *argv[]) { /* Single Photon Detector commands */ - - // if (doc.HasMember("nSigma")) { - // nSigma=doc["nSigma"].GetUint(); - - // } - + nSigma=5; + if (doc.HasMember("nSigma")) { + if (doc["nSigma"].IsInt()) + nSigma=doc["nSigma"].GetInt(); + mt->setNSigma(nSigma); + } + + emin=-1; + emax=-1; + if (doc.HasMember("energyRange")) { + if (doc["energyRange"].IsArray()) { + if (doc["energyRange"].Size() > 0 ) + if (doc["energyRange"][0].IsInt()) + emin=doc["energyRange"][0].GetInt(); + + if (doc["energyRange"].Size() > 1 ) + if (doc["energyRange"][1].IsInt()) + emax=doc["energyRange"][1].GetUint(); + } + } + if (doc.HasMember("eMin")) { + if (doc["eMin"][1].IsInt()) + emin=doc["eMin"].GetInt(); + } + if (doc.HasMember("eMax")) { + if (doc["eMax"][1].IsInt()) + emin=doc["eMax"].GetInt(); + } + mt->setEnergyRange(emin,emax); - // if (doc.HasMember("energyRange")) { - // emin=doc["energyRange"][0].GetUint(); - // emax=doc["energyRange"][1].GetUint(); - - // } - /* interpolating detector commands */ - // if (doc.HasMember("nSubPixels")) { - // nsubPixels=doc["nSubPixels"].GetUint(); - // } + if (doc.HasMember("nSubPixels")) { + if (doc["nSubPixels"].IsUint()) + nSubPixels=doc["nSubPixels"].GetUint(); + mt->setNSubPixels(nSubPixels); + } - //if (doc.HasMember("interpolation")) { - // intMode_s=doc["intepolation"].GetString(); - //} newFrame=0; zmqsocket->CloseHeaderMessage(); @@ -514,6 +610,7 @@ int main(int argc, char *argv[]) { iframe++; + } // exiting infinite loop diff --git a/slsDetectorCalibration/multiThreadedAnalogDetector.h b/slsDetectorCalibration/multiThreadedAnalogDetector.h index 14f7d130a..6657f6720 100644 --- a/slsDetectorCalibration/multiThreadedAnalogDetector.h +++ b/slsDetectorCalibration/multiThreadedAnalogDetector.h @@ -1,6 +1,8 @@ #ifndef MULTITHREADED_ANALOG_DETECTOR_H #define MULTITHREADED_ANALOG_DETECTOR_H +#define MAXTHREADS 1000 + #include #include #include @@ -14,11 +16,8 @@ #include #include -//#include "analogDetector.h" -#include "singlePhotonDetector.h" -//#include "interpolatingDetector.h" +#include "analogDetector.h" #include "circularFifo.h" -#include "slsInterpolation.h" #include #include #include @@ -57,6 +56,9 @@ public: return fMode; }; virtual double setThreshold(double th) {return det->setThreshold(th);}; + + + virtual void setROI(int xmin, int xmax, int ymin, int ymax) {det->setROI(xmin,xmax,ymin,ymax);}; virtual int setDetectorMode(int dm) { if (dm>=0) { @@ -111,43 +113,7 @@ public: virtual void clearImage(){det->clearImage();}; - virtual void prepareInterpolation(int &ok){ - slsInterpolation *interp=(slsInterpolation*)det->getInterpolation(); - if (interp) - interp->prepareInterpolation(ok); - } - - virtual int *getFlatField(){ - slsInterpolation *interp=(slsInterpolation*)det->getInterpolation(); - if (interp) - return interp->getFlatField(); - else - return NULL; - } - - virtual int *setFlatField(int *ff, int nb, double emin, double emax){ - slsInterpolation *interp=(slsInterpolation*)det->getInterpolation(); - if (interp) - return interp->setFlatField(ff, nb, emin, emax); - else - return NULL; - } - - virtual int *getFlatField(int &nb, double emi, double ema){ - slsInterpolation *interp=(slsInterpolation*)det->getInterpolation(); - int *ff; - if (interp) { - ff=interp->getFlatField(nb,emi,ema); - // cout << "tdgff* ff has " << nb << " bins " << endl; - return ff; - } else - return NULL; - } - - virtual char *getInterpolation() { - return det->getInterpolation(); - } - + virtual void setPedestal(double *ped, double *rms=NULL, int m=-1){det->setPedestal(ped,rms,m);}; @@ -173,9 +139,8 @@ public: FILE *getFilePointer(){return det->getFilePointer();}; - void setMutex(pthread_mutex_t *fmutex){det->setMutex(fmutex);}; -private: +protected: analogDetector *det; int fMode; int dMode; @@ -219,7 +184,7 @@ private: class multiThreadedAnalogDetector { public: - multiThreadedAnalogDetector(analogDetector *d, int n, int fs=1000) : nThreads(n), ithread(0) { + multiThreadedAnalogDetector(analogDetector *d, int n, int fs=1000) : stop(0), nThreads(n), ithread(0) { dd[0]=d; if (nThreads==1) dd[0]->setId(100); @@ -239,6 +204,7 @@ public: image=new int[nn]; ff=NULL; ped=NULL; + cout << "Ithread is " << ithread << endl; } ~multiThreadedAnalogDetector() { @@ -251,13 +217,15 @@ public: } - int setFrameMode(int fm) { int ret; for (int i=0; isetFrameMode(fm); return ret;}; - double setThreshold(int fm) { double ret; for (int i=0; isetThreshold(fm); return ret;}; - int setDetectorMode(int dm) { int ret; for (int i=0; isetDetectorMode(dm); return ret;}; - void setROI(int xmin, int xmax, int ymin, int ymax) { for (int i=0; isetROI(xmin, xmax,ymin,ymax);}; + virtual int setFrameMode(int fm) { int ret; for (int i=0; isetFrameMode(fm);} return ret;}; + virtual double setThreshold(int fm) { double ret; for (int i=0; isetThreshold(fm); return ret;}; + virtual int setDetectorMode(int dm) { int ret; for (int i=0; isetDetectorMode(dm); return ret;}; + virtual void setROI(int xmin, int xmax, int ymin, int ymax) { for (int i=0; isetROI(xmin, xmax,ymin,ymax);}; - void newDataSet(){for (int i=0; inewDataSet();}; - int *getImage(int &nnx, int &nny, int &ns) { + + virtual void newDataSet(){for (int i=0; inewDataSet();}; + + virtual int *getImage(int &nnx, int &nny, int &ns) { int *img; // int nnx, nny, ns; int nn=dets[0]->getImageSize(nnx, nny, ns); @@ -280,7 +248,7 @@ public: } - void clearImage() { + virtual void clearImage() { for (int ii=0; iiclearImage(); @@ -315,10 +283,11 @@ public: virtual void StartThreads() { - for (int i=0; iStartThread(); - } + + } @@ -329,7 +298,7 @@ public: } - int isBusy() { + virtual int isBusy() { int ret=0, ret1; for (int i=0; iisBusy(); @@ -340,125 +309,21 @@ public: } - bool pushData(char* &ptr) { + virtual bool pushData(char* &ptr) { dets[ithread]->pushData(ptr); } - bool popFree(char* &ptr) { + virtual bool popFree(char* &ptr) { + cout << ithread << endl; dets[ithread]->popFree(ptr); } - int nextThread() { + virtual int nextThread() { ithread++; if (ithread==nThreads) ithread=0; return ithread; } - virtual void prepareInterpolation(int &ok){ - getFlatField(); //sum up all etas - setFlatField(); //set etas to all detectors - for (int i=0; iprepareInterpolation(ok); - } - } - - virtual int *getFlatField(){ - int nb=0; - double emi, ema; - int *f0; - slsInterpolation* inte=(slsInterpolation*)dets[0]->getInterpolation(); - if (inte) { - if (inte->getFlatField(nb,emi,ema)) { - if (ff) delete [] ff; - ff=new int[nb*nb]; - for (int i=0; igetInterpolation(); - f0=inte->getFlatField(); - if (f0) { - // cout << "ff " << i << endl; - for (int ib=0; ib0) */ - /* cout << i << " " << ib << " " << f0[ib] << " " << ff[ib] << endl; */ - - } - } - } - return ff; - } - } - return NULL; - } - - - virtual int *setFlatField(int *h=NULL, int nb=-1, double emin=1, double emax=0){ - //int nb=0; - double emi, ema; - slsInterpolation* inte=(slsInterpolation*)dets[0]->getFlatField(nb,emi,ema); - if (inte) { - if (h==NULL) h=ff; - for (int i=0; isetFlatField(h, nb, emin, emax); - } - } - return NULL; - }; - - void *writeFlatField(const char * imgname){ - - int nb=0; - double emi, ema; - slsInterpolation* inte=(slsInterpolation*)dets[0]->getFlatField(nb,emi,ema); - if (inte) { - if (getFlatField()) { - // cout << "mtwff* ff has " << nb << " bins " << endl; - float *gm=new float[nb*nb]; - if (gm) { - for (int ix=0; ixgetFlatField(nb,emin,emax); - if (inte) { - uint32 nnx; - uint32 nny; - float *gm=ReadFromTiff(imgname, nnx, nny); - if (ff) delete [] ff; - if (nnx>nb) nb=nnx; - if (nny>nb) nb=nny; - ff=new int[nb*nb]; - - for (int ix=0; ixgetDetectorSize(nx,ny); @@ -520,7 +385,7 @@ public: }; - void *readPedestal(const char * imgname, int nb=-1, double emin=1, double emax=0){ + virtual void *readPedestal(const char * imgname, int nb=-1, double emin=1, double emax=0){ int nx, ny; dets[0]->getDetectorSize(nx,ny); @@ -549,7 +414,7 @@ public: \param f file pointer \returns current file pointer */ - FILE *setFilePointer(FILE *f){ + virtual FILE *setFilePointer(FILE *f){ for (int i=0; isetFilePointer(f); //dets[i]->setMutex(&fmutex); @@ -557,18 +422,18 @@ public: return dets[0]->getFilePointer(); }; -/** gets file pointer where to write the clusters to - \returns current file pointer -*/ -FILE *getFilePointer(){return dets[0]->getFilePointer();}; + /** gets file pointer where to write the clusters to + \returns current file pointer + */ + virtual FILE *getFilePointer(){return dets[0]->getFilePointer();}; - private: + protected: bool stop; const int nThreads; - threadedAnalogDetector *dets[20]; - analogDetector *dd[20]; + threadedAnalogDetector *dets[MAXTHREADS]; + analogDetector *dd[MAXTHREADS]; int ithread; int *image; int *ff; diff --git a/slsDetectorCalibration/multiThreadedCountingDetector.h b/slsDetectorCalibration/multiThreadedCountingDetector.h new file mode 100644 index 000000000..9564151ca --- /dev/null +++ b/slsDetectorCalibration/multiThreadedCountingDetector.h @@ -0,0 +1,55 @@ +#ifndef MULTITHREADED_COUNTING_DETECTOR_H +#define MULTITHREADED_COUNTING_DETECTOR_H + + +#include "singlePhotonDetector.h" +#include "multiThreadedAnalogDetector.h" +//#include + + +using namespace std; + + +class threadedCountingDetector : public threadedAnalogDetector +{ +public: + threadedCountingDetector(singlePhotonDetector *d, int fs=10000) : threadedAnalogDetector(d,fs) {}; + + + virtual double setNSigma(double n) {return ((singlePhotonDetector*)det)->setNSigma(n);}; + virtual void setEnergyRange(double emi, double ema) {return ((singlePhotonDetector*)det)->setEnergyRange(emi,ema);}; + + +}; + + + +class multiThreadedCountingDetector : public multiThreadedAnalogDetector +{ +public: + multiThreadedCountingDetector(singlePhotonDetector *d, int n, int fs=1000) : multiThreadedAnalogDetector(d,n,fs) { }; + + virtual double setNSigma(double n) {double ret; for (int i=0; isetNSigma(n); return ret;}; + virtual void setEnergyRange(double emi, double ema) {for (int i=0; isetEnergyRange(emi,ema);}; + +}; + +#endif + + + + + + + + + + + + + + + + + + diff --git a/slsDetectorCalibration/multiThreadedInterpolatingDetector.h b/slsDetectorCalibration/multiThreadedInterpolatingDetector.h new file mode 100644 index 000000000..02d4cb415 --- /dev/null +++ b/slsDetectorCalibration/multiThreadedInterpolatingDetector.h @@ -0,0 +1,136 @@ +#ifndef MULTITHREADED_INTERPOLATING_DETECTOR_H +#define MULTITHREADED_INTERPOLATING_DETECTOR_H + + +#include "interpolatingDetector.h" +#include "multiThreadedCountingDetector.h" +//#include + + +using namespace std; + + +class threadedInterpolatingDetector : public threadedCountingDetector +{ +public: + threadedInterpolatingDetector(interpolatingDetector *d, int fs=10000) : threadedCountingDetector(d,fs) {}; + + virtual void prepareInterpolation(int &ok){ + slsInterpolation *interp=(slsInterpolation*)((interpolatingDetector*)det)->getInterpolation(); + if (interp) + interp->prepareInterpolation(ok); + } + + virtual int *getFlatField(){ + slsInterpolation *interp=(slsInterpolation*)((interpolatingDetector*)det)->getInterpolation(); + if (interp) + return interp->getFlatField(); + else + return NULL; + } + + virtual int *setFlatField(int *ff, int nb, double emin, double emax){ + slsInterpolation *interp=(slsInterpolation*)((interpolatingDetector*)det)->getInterpolation(); + if (interp) + return interp->setFlatField(ff, nb, emin, emax); + else + return NULL; + } + + void *writeFlatField(const char * imgname) { + + slsInterpolation *interp=(slsInterpolation*)((interpolatingDetector*)det)->getInterpolation(); + if (interp) + interp->writeFlatField(imgname); + + } + void *readFlatField(const char * imgname, int nb=-1, double emin=1, double emax=0){ + slsInterpolation *interp=(slsInterpolation*)((interpolatingDetector*)det)->getInterpolation(); + if (interp) + interp->readFlatField(imgname, nb, emin, emax); + } + + virtual int *getFlatField(int &nb, double emi, double ema){ + slsInterpolation *interp=(slsInterpolation*)((interpolatingDetector*)det)->getInterpolation(); + int *ff=NULL; + if (interp) { + ff=interp->getFlatField(nb,emi,ema); + } + return ff; + } + + virtual char *getInterpolation() { + return ((interpolatingDetector*)det)->getInterpolation(); + } + + virtual void resetFlatField() { ((interpolatingDetector*)det)->resetFlatField();}; + + + virtual int setNSubPixels(int ns) { return ((interpolatingDetector*)det)->setNSubPixels(ns);}; + + + +}; + + + +class multiThreadedInterpolatingDetector : public multiThreadedCountingDetector +{ +public: + multiThreadedInterpolatingDetector(interpolatingDetector *d, int n, int fs=1000) : multiThreadedCountingDetector(d,n,fs) { }; + + virtual void prepareInterpolation(int &ok){ + /* getFlatField(); //sum up all etas */ + /* setFlatField(); //set etas to all detectors */ + /* for (int i=0; iprepareInterpolation(ok); + // } + } + + virtual int *getFlatField(){ + return ((threadedInterpolatingDetector*)dets[0])->getFlatField(); + } + + virtual int *getFlatField(int &nb, double emi, double ema){ + return ((threadedInterpolatingDetector*)dets[0])->getFlatField(nb,emi,ema); + } + + virtual int *setFlatField(int *h=NULL, int nb=-1, double emin=1, double emax=0){ + return ((threadedInterpolatingDetector*)dets[0])->setFlatField(h,nb,emin,emax); + }; + + void *writeFlatField(const char * imgname){ + ((threadedInterpolatingDetector*)dets[0])->writeFlatField(imgname); + }; + + + void *readFlatField(const char * imgname, int nb=-1, double emin=1, double emax=0){ + ((threadedInterpolatingDetector*)dets[0])->readFlatField(imgname, nb, emin, emax); + }; + + + virtual int setNSubPixels(int ns) { return ((threadedInterpolatingDetector*)dets[0])->setNSubPixels(ns);}; + virtual void resetFlatField() {((threadedInterpolatingDetector*)dets[0])->resetFlatField();}; + + +}; + +#endif + + + + + + + + + + + + + + + + + + diff --git a/slsDetectorCalibration/singlePhotonDetector.h b/slsDetectorCalibration/singlePhotonDetector.h index 8bae7ff62..b35682150 100644 --- a/slsDetectorCalibration/singlePhotonDetector.h +++ b/slsDetectorCalibration/singlePhotonDetector.h @@ -63,6 +63,8 @@ public analogDetector { + fm=new pthread_mutex_t ; + eventMask=new eventType*[ny]; for (int i=0; i { // cluster=clusters; setClusterSize(clusterSize); + fm=orig->fm; quad=UNDEFINED_QUADRANT; tot=0; @@ -138,7 +141,7 @@ public analogDetector { \param n number of sigma to be set (0 or negative gets) \returns actual number of sigma parameter */ - double setNSigma(double n=-1){if (n>0) nSigma=n; return nSigma;} + double setNSigma(double n=-1){if (n>=0) nSigma=n; return nSigma;} /** sets/gets cluster size \param n cluster size to be set, (0 or negative gets). If even is incremented by 1. @@ -217,230 +220,224 @@ public analogDetector { } for (int ix=xmin; ix::getNPhotons(data,ix,iy);//val/thr;// - if (nn>0) { - /* if (nph[ix+nx*iy]>0) { */ - /* cout << nph[ix+nx*iy] << " => "; */ - /* cout << ix << " " << iy << " " << val << " "<< nn << " " << nph[ix+nx*iy]+nn << endl; */ - /* } */ - nph[ix+nx*iy]+=nn; - /* if (nph[ix+nx*iy]>1) { */ - /* cout << nph[ix+nx*iy] << " => "; */ - /* } */ - rest[iy][ix]=(val-nn*thr+0.5*thr); - // if (ix==150 && iy==100) - // cout << ix << " " << iy << " " << val << " "<< nn << " " << rest[iy][ix] << " " << nph[ix+nx*iy] << endl; - nphFrame+=nn; - nphTot+=nn; - } else - rest[iy][ix]=val; - + if (det->isGood(ix,iy)) { + val=subtractPedestal(data,ix,iy, cm); + + nn=analogDetector::getNPhotons(data,ix,iy);//val/thr;// + if (nn>0) { + nph[ix+nx*iy]+=nn; + rest[iy][ix]=(val-nn*thr);//?+0.5*thr + nphFrame+=nn; + nphTot+=nn; + } else + rest[iy][ix]=val; + + } } } - // } for (int ix=xmin; ix0.25*thr) { - eventMask[iy][ix]=NEIGHBOUR; - for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) { - for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) { - if ((iy+ir)>=0 && (iy+ir)=0 && (ix+ic)set_data(rest[iy+ir][ix+ic], ic, ir); + if (det->isGood(ix,iy)) { + eventMask[iy][ix]=PEDESTAL; + max=0; + tl=0; + tr=0; + bl=0; + br=0; + tot=0; + quadTot=0; + + if (rest[iy][ix]>0.25*thr) { + eventMask[iy][ix]=NEIGHBOUR; + for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) { + for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) { + if ((iy+ir)>=0 && (iy+ir)=0 && (ix+ic)set_data(rest[iy+ir][ix+ic], ic, ir); - - v=rest[iy+ir][ix+ic];//clusters->get_data(ic,ir); - 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; + + v=rest[iy+ir][ix+ic];//clusters->get_data(ic,ir); + 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 (ir==0 && ic==0) { - //} } } - } - - if (rest[iy][ix]>=max) { - if (bl>=br && bl>=tl && bl>=tr) { - quad=BOTTOM_LEFT; - quadTot=bl; - } else if (br>=bl && br>=tl && br>=tr) { - quad=BOTTOM_RIGHT; - quadTot=br; - } else if (tl>=br && tl>=bl && tl>=tr) { - quad=TOP_LEFT; - quadTot=tl; - } else if (tr>=bl && tr>=tl && tr>=br) { - quad=TOP_RIGHT; + + if (rest[iy][ix]>=max) { + if (bl>=br && bl>=tl && bl>=tr) { + quad=BOTTOM_LEFT; + quadTot=bl; + } else if (br>=bl && br>=tl && br>=tr) { + quad=BOTTOM_RIGHT; + quadTot=br; + } else if (tl>=br && tl>=bl && tl>=tr) { + quad=TOP_LEFT; + quadTot=tl; + } else if (tr>=bl && tr>=tl && tr>=br) { + quad=TOP_RIGHT; quadTot=tr; - } - - rms=getPedestalRMS(ix,iy); - tthr=nSigma*rms; + } + + if (nSigma==0) { + tthr=thr; + tthr1=thr; + tthr2=thr; + } else { + + rms=getPedestalRMS(ix,iy); + tthr=nSigma*rms; - tthr1=nSigma*sqrt(clusterSize*clusterSizeY)*rms; - tthr2=nSigma*sqrt((clusterSize+1)/2.*((clusterSizeY+1)/2.))*rms; - - - if (thr>tthr) tthr=thr-tthr; - if (thr>tthr1) tthr1=tthr-tthr1; - if (thr>tthr2) tthr2=tthr-tthr2; - - if (tot>tthr1 || quadTot>tthr2 || max>tthr) { - eventMask[iy][ix]=PHOTON; - nph[ix+nx*iy]++; - // if (ix==150 && iy==100) - // if (iy<399) - // cout << "** " << ix << " " << iy << " " << quadTot << endl; - rest[iy][ix]-=thr; - nphFrame++; - nphTot++; + tthr1=nSigma*sqrt(clusterSize*clusterSizeY)*rms; + tthr2=nSigma*sqrt((clusterSize+1)/2.*((clusterSizeY+1)/2.))*rms; + + + if (thr>2*tthr) tthr=thr-tthr; + if (thr>2*tthr1) tthr1=tthr-tthr1; + if (thr>2*tthr2) tthr2=tthr-tthr2; + + } - } + if (tot>tthr1 || quadTot>tthr2 || max>tthr) { + eventMask[iy][ix]=PHOTON; + nph[ix+nx*iy]++; + rest[iy][ix]-=thr; + nphFrame++; + nphTot++; + + } + } } } } } - // cout << iframe << " " << nphFrame << " " << nphTot << endl; - //cout << iframe << " " << nph << endl; } else return getClusters(data, nph); } return NULL; }; - /** finds event type for pixel and fills cluster structure. The algorithm loops only if the evenMask for this pixel is still undefined. - if pixel or cluster around it are above threshold (nsigma*pedestalRMS) cluster is filled and pixel mask is PHOTON_MAX (if maximum in cluster) or NEIGHBOUR; If PHOTON_MAX, the elements of the cluster are also set as NEIGHBOURs in order to speed up the looping - if below threshold the pixel is either marked as PEDESTAL (and added to the pedestal calculator) or NEGATIVE_PEDESTAL is case it's lower than -threshold, otherwise the pedestal average would drift to negative values while it should be 0. + /* /\** finds event type for pixel and fills cluster structure. The algorithm loops only if the evenMask for this pixel is still undefined. */ + /* if pixel or cluster around it are above threshold (nsigma*pedestalRMS) cluster is filled and pixel mask is PHOTON_MAX (if maximum in cluster) or NEIGHBOUR; If PHOTON_MAX, the elements of the cluster are also set as NEIGHBOURs in order to speed up the looping */ + /* if below threshold the pixel is either marked as PEDESTAL (and added to the pedestal calculator) or NEGATIVE_PEDESTAL is case it's lower than -threshold, otherwise the pedestal average would drift to negative values while it should be 0. */ - /param data pointer to the data - /param ix pixel x coordinate - /param iy pixel y coordinate - /param cm enable(1)/disable(0) common mode subtraction (if defined). - /returns event type for the given pixel - */ - eventType getEventType(char *data, int ix, int iy, int cm=0) { + /* /param data pointer to the data */ + /* /param ix pixel x coordinate */ + /* /param iy pixel y coordinate */ + /* /param cm enable(1)/disable(0) common mode subtraction (if defined). */ + /* /returns event type for the given pixel */ + /* *\/ */ + /* eventType getEventType(char *data, int ix, int iy, int cm=0) { */ - // eventType ret=PEDESTAL; - double max=0, tl=0, tr=0, bl=0,br=0, v; - // cout << iframe << endl; + /* // eventType ret=PEDESTAL; */ + /* double max=0, tl=0, tr=0, bl=0,br=0, v; */ + /* // cout << iframe << endl; */ - int cy=(clusterSizeY+1)/2; - int cs=(clusterSize+1)/2; - double val; - tot=0; - quadTot=0; - quad=UNDEFINED_QUADRANT; + /* int cy=(clusterSizeY+1)/2; */ + /* int cs=(clusterSize+1)/2; */ + /* double val; */ + /* tot=0; */ + /* quadTot=0; */ + /* quad=UNDEFINED_QUADRANT; */ - if (iframex=ix; - clusters->y=iy; - clusters->rms=getPedestalRMS(ix,iy); - clusters->ped=getPedestal(ix,iy, cm); + /* clusters->x=ix; */ + /* clusters->y=iy; */ + /* clusters->rms=getPedestalRMS(ix,iy); */ + /* clusters->ped=getPedestal(ix,iy, cm); */ - for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) { - for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) { - if ((iy+ir)>=0 && (iy+ir)=0 && (ix+ic)=0 && (iy+ir)=0 && (ix+ic)set_data(v, ic, ir); - // v=clusters->get_data(ic,ir); - 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; + /* clusters->set_data(v, ic, ir); */ + /* // v=clusters->get_data(ic,ir); */ + /* 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->rms) - eventMask[iy][ix]=NEGATIVE_PEDESTAL; - } - } - } - } + /* if (v>max) { */ + /* max=v; */ + /* } */ + /* if (ir==0 && ic==0) { */ + /* if (v<-nSigma*clusters->rms) */ + /* eventMask[iy][ix]=NEGATIVE_PEDESTAL; */ + /* } */ + /* } */ + /* } */ + /* } */ - if (bl>=br && bl>=tl && bl>=tr) { - quad=BOTTOM_LEFT; - quadTot=bl; - } else if (br>=bl && br>=tl && br>=tr) { - quad=BOTTOM_RIGHT; - quadTot=br; - } else if (tl>=br && tl>=bl && tl>=tr) { - quad=TOP_LEFT; - quadTot=tl; - } else if (tr>=bl && tr>=tl && tr>=br) { - quad=TOP_RIGHT; - quadTot=tr; - } + /* if (bl>=br && bl>=tl && bl>=tr) { */ + /* quad=BOTTOM_LEFT; */ + /* quadTot=bl; */ + /* } else if (br>=bl && br>=tl && br>=tr) { */ + /* quad=BOTTOM_RIGHT; */ + /* quadTot=br; */ + /* } else if (tl>=br && tl>=bl && tl>=tr) { */ + /* quad=TOP_LEFT; */ + /* quadTot=tl; */ + /* } else if (tr>=bl && tr>=tl && tr>=br) { */ + /* quad=TOP_RIGHT; */ + /* quadTot=tr; */ + /* } */ - if (max>nSigma*clusters->rms || tot>sqrt(clusterSizeY*clusterSize)*nSigma*clusters->rms || quadTot>cy*cs*nSigma*clusters->rms) { - if (clusters->get_data(0,0)>=max) { - eventMask[iy][ix]=PHOTON_MAX; - } else { - eventMask[iy][ix]=PHOTON; - } - } else if (eventMask[iy][ix]==PEDESTAL) { - if (cm==0) { - if (det) - val=dataSign*det->getValue(data, ix, iy); - else - val=((double**)data)[iy][ix]; - addToPedestal(val,ix,iy); - } - } + /* if (max>nSigma*clusters->rms || tot>sqrt(clusterSizeY*clusterSize)*nSigma*clusters->rms || quadTot>cy*cs*nSigma*clusters->rms) { */ + /* if (clusters->get_data(0,0)>=max) { */ + /* eventMask[iy][ix]=PHOTON_MAX; */ + /* } else { */ + /* eventMask[iy][ix]=PHOTON; */ + /* } */ + /* } else if (eventMask[iy][ix]==PEDESTAL) { */ + /* if (cm==0) { */ + /* if (det) */ + /* val=dataSign*det->getValue(data, ix, iy); */ + /* else */ + /* val=((double**)data)[iy][ix]; */ + /* addToPedestal(val,ix,iy); */ + /* } */ + /* } */ - return eventMask[iy][ix]; + /* return eventMask[iy][ix]; */ - }; + /* }; */ @@ -482,71 +479,71 @@ int *getClusters(char *data, int *ph=NULL) { for (int ix=xmin; ixisGood(ix,iy)) { + max=0; + tl=0; + tr=0; + bl=0; + br=0; + tot=0; + quadTot=0; + quad=UNDEFINED_QUADRANT; + - (clusters+nph)->rms=getPedestalRMS(ix,iy); - // cluster=clusters+nph; + eventMask[iy][ix]=PEDESTAL; + + + (clusters+nph)->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++) { - - 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; - } - - - if (ir==0 && ic==0) { - if (*v<-nSigma*(clusters+nph)->rms) - eventMask[iy][ix]=NEGATIVE_PEDESTAL; - } - + for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) { + for (int 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; + } + + + if (ir==0 && ic==0) { + if (*v<-nSigma*(clusters+nph)->rms) + eventMask[iy][ix]=NEGATIVE_PEDESTAL; + } + } - } - - if (bl>=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; + } + + if (bl>=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; - } + } 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 (max>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) { + if (max>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; @@ -558,7 +555,6 @@ int *getClusters(char *data, int *ph=NULL) { (clusters+nph)->set_data(val[iy+ir][ix+ic],ic,ir); } } - // cout << (clusters+nph)->iframe << " " << ix << " " << nph << " " << tot << " " << (clusters+nph)->quadTot << endl; good=1; if (eMin>0 && tot0 && tot>eMax) good=0; @@ -566,13 +562,14 @@ int *getClusters(char *data, int *ph=NULL) { nph++; image[iy*nx+ix]++; } - + } else { eventMask[iy][ix]=PHOTON; } - } else if (eventMask[iy][ix]==PEDESTAL) { - addToPedestal(data,ix,iy,cm); + } else if (eventMask[iy][ix]==PEDESTAL) { + addToPedestal(data,ix,iy,cm); + } } } } @@ -704,6 +701,7 @@ void writeClusters(FILE *f){for (int i=0; iwrite(f void setEnergyRange(double emi, double ema){eMin=emi; eMax=ema;}; void getEnergyRange(double &emi, double &ema){emi=eMin; ema=eMax;}; + void setMutex(pthread_mutex_t *m){fm=m;}; protected: @@ -721,6 +719,7 @@ void writeClusters(FILE *f){for (int i=0; iwrite(f int nphTot; int nphFrame; + pthread_mutex_t *fm; };