From d1ef6bf94be980cc865b766648be71bc267ec3d2 Mon Sep 17 00:00:00 2001 From: Anna Bergamaschi Date: Tue, 30 Jan 2018 10:24:56 +0100 Subject: [PATCH] some fixes in the cluster finder --- .../gotthardDoubleModuleDataNew.h | 8 +- .../gotthard25umZmqAnalysis.C | 46 ++-- .../etaInterpolationAdaptiveBins.h | 203 +++++++++++------- .../interpolations/etaInterpolationBase.h | 33 +-- .../moenchExecutables/Makefile.cluster_finder | 18 +- .../moenchExecutables/Makefile.tiff_to_th2f | 6 +- .../moench03ClusterFinder.cpp | 170 +++++++++------ .../moench03Interpolation.cpp | 12 +- .../moenchExecutables/readClusters.C | 37 ++-- slsDetectorCalibration/singlePhotonDetector.h | 2 +- 10 files changed, 323 insertions(+), 212 deletions(-) diff --git a/slsDetectorCalibration/dataStructures/gotthardDoubleModuleDataNew.h b/slsDetectorCalibration/dataStructures/gotthardDoubleModuleDataNew.h index 1daa9db82..5aea9c75f 100644 --- a/slsDetectorCalibration/dataStructures/gotthardDoubleModuleDataNew.h +++ b/slsDetectorCalibration/dataStructures/gotthardDoubleModuleDataNew.h @@ -50,9 +50,13 @@ public: for(int ix=0; ixsetROI(0,512,0,1); + singlePhotonDetector *filter=new singlePhotonDetector(det,3, 5, 1, 0, 1000, 100); + filter->setROI(0,2560,0,1); char *buff;//[2*(48+1280*2)]; char *buff0; char *buff1; @@ -82,7 +82,7 @@ int main(int argc, char *argv[]){ mt->StartThreads(); mt->popFree(buff); buff0=buff; - buff1=buff+48+1280*2; + buff1=buff+offset*2+1280*2; int photons[1280*2]; int nf=0; int ok=0; @@ -249,7 +249,6 @@ int main(int argc, char *argv[]){ int runmax=atoi(argv[4]); char *outdir=argv[5]; sprintf(ff,"%s/%s",indir,fformat); - // strcpy(fformat,"/external_pool/gotthard_data/datadir_gotthardI/bchip074075/20170731/Xray/xray_15kV_200uA_5us_d%d_f000000000000_0.raw"); // sprintf(fname0,fformat,0,0); // sprintf(fname1,fformat,1,1); @@ -277,21 +276,37 @@ int main(int argc, char *argv[]){ #ifdef ZMQ + int end_of_acquisition; while(1) { - if ((!zmqsocket0->ReceiveHeader(0, acqIndex0, frameIndex0, subframeIndex0, filename0, fileindex0)) && (!zmqsocket1->ReceiveHeader(0, acqIndex1, frameIndex1, subframeIndex1, filename1, fileindex1))){ + end_of_acquisition=0; + // cout << "Receive header " << nf << endl; + if (!zmqsocket0->ReceiveHeader(0, acqIndex0, frameIndex0, subframeIndex0, filename0, fileindex0)) { + + cout << "************************************************************************** packet0!*****************************"<< endl; + end_of_acquisition++; + } + if (!zmqsocket1->ReceiveHeader(0, acqIndex1, frameIndex1, subframeIndex1, filename1, fileindex1)) { + cout << "************************************************************************** packet1!*****************************"<< endl; + end_of_acquisition++; + } + + // if ((!zmqsocket0->ReceiveHeader(0, acqIndex0, frameIndex0, subframeIndex0, filename0, fileindex0)) && (!zmqsocket1->ReceiveHeader(0, acqIndex1, frameIndex1, subframeIndex1, filename1, fileindex1))){ + if (end_of_acquisition) { + cout << "************************************************************************** END OF FRAME" << end_of_acquisition << " !*****************************"<< endl; + // return 0; while (mt->isBusy()) {;} image=filter->getImage(); if (image) { - fout=fopen(ofname,"w"); - //cout << nf << "*****************" << endl; - for (int i=0; i<512; i++) { - fprintf(fout,"%d %d\n",i,image[i]); - dout[i]=image[i]; + //fout=fopen(ofname,"w"); + cout << nf << "*****************" << endl; + for (int i=0; i<2560; i++) { + // // fprintf(fout,"%d %d\n",i,image[i]); + dout[i]=image[i]; } - fclose(fout);; + // // fclose(fout);; } @@ -330,9 +345,10 @@ int main(int argc, char *argv[]){ cout << "different frame indexes " << frameIndex0 << " and " << frameIndex1 << endl; - + // cout << "Receive data " << nf << endl; length = zmqsocket0->ReceiveData(0, buff0, size/2); - length += zmqsocket0->ReceiveData(0, buff1, size/2); + length += zmqsocket1->ReceiveData(0, buff1, size/2); + irun=fileindex0; sprintf(ofname,"%s_%d.ph",filename0.c_str(),fileindex0); @@ -344,7 +360,7 @@ int main(int argc, char *argv[]){ mt->nextThread(); mt->popFree(buff); buff0=buff; - buff1=buff+48+1280*2; + buff1=buff+offset*2+1280*2; nf++; diff --git a/slsDetectorCalibration/interpolations/etaInterpolationAdaptiveBins.h b/slsDetectorCalibration/interpolations/etaInterpolationAdaptiveBins.h index 7623bc7e9..2100380f5 100644 --- a/slsDetectorCalibration/interpolations/etaInterpolationAdaptiveBins.h +++ b/slsDetectorCalibration/interpolations/etaInterpolationAdaptiveBins.h @@ -3,17 +3,19 @@ #include "tiffIO.h" -#include "etaInterpolationBase.h" +//#include "etaInterpolationBase.h" +#include "etaInterpolationPosXY.h" + +class etaInterpolationAdaptiveBins : public etaInterpolationPosXY { -class etaInterpolationAdaptiveBins : public etaInterpolationBase{ private: - double calcDiff(double avg, double *hx, double *hy) { + double calcDiff(double avg, float *hx, float *hy) { double p_tot=0; double diff=0; double bsize=1./nSubPixels; - for (int ipx=0; ipx=nSubPixels; ipx++) { - for (int ipy=0; ipy=nSubPixels; ipy++) { + for (int ipx=0; ipx=((ipy)*bsize) && hhy[ibx+iby*nbeta]<((ipy+1)*bsize)) { - hx[ibx]+=heta[ibx+iby*nbeta]; - } - if ( hhx[ibx+iby*nbeta]>=((ipy)*bsize) && hhx[ibx+iby*nbeta]<((ipy+1)*bsize)) { - hy[iby]+=heta[ibx+iby*nbeta]; - } - } - } - - - } - - hix[0]=hx[0]; - hiy[0]=hy[0]; - for (int ib=1; ib=((ipy)*bsize) && hhy[ibx+iby*nbeta]<((ipy+1)*bsize)) { - newhhx[ibx+iby*nbeta]=hix[ibx]/((double)tot_eta_x); - } - if ( hhx[ibx+iby*nbeta]>=((ipy)*bsize) && hhx[ibx+iby*nbeta]<((ipy+1)*bsize)) { - newhhy[ibx+iby*nbeta]=hiy[iby]/((double)tot_eta_y); - } - } + hx[ibx]=0; + hy[ibx]=0; } + + tot_eta_x=0; + tot_eta_y=0; + // cout << ipy << " " << ((ipy)*bsize) << " " << ((ipy+1)*bsize) << endl; + for (int ibx=0; ibx=((ipy)*bsize) && hhy[ibx+iby*nbeta]<=((ipy+1)*bsize)) { + hx[ibx]+=heta[ibx+iby*nbeta]; + tot_eta_x+=heta[ibx+iby*nbeta]; + } + + + if (hhx[ibx+iby*nbeta]>=((ipy)*bsize) && hhx[ibx+iby*nbeta]<=((ipy+1)*bsize)) { + hy[iby]+=heta[ibx+iby*nbeta]; + tot_eta_y+=heta[ibx+iby*nbeta]; + } + } + } + + + hix[0]=hx[0]; + hiy[0]=hy[0]; + for (int ib=1; ib=((ipy)*bsize) && hhy[ibx+iby*nbeta]<=((ipy+1)*bsize)) { + newhhx[ibx+iby*nbeta]=hix[ibx]/((double)tot_eta_x); + if (newhhx[ibx+iby*nbeta]>1) cout << "***"<< ibx << " " << iby << newhhx[ibx+iby*nbeta] << endl; + // if (ipy==3 && ibx==10) cout << newhhx[ibx+iby*nbeta] << " " << hix[ibx] << " " << ibx+iby*nbeta << endl; + } + if (hhx[ibx+iby*nbeta]>=((ipy)*bsize) && hhx[ibx+iby*nbeta]<=((ipy+1)*bsize)) { + newhhy[ibx+iby*nbeta]=hiy[iby]/((double)tot_eta_y); + if (newhhy[ibx+iby*nbeta]>1) cout << "***"<< ibx << " " << iby << newhhy[ibx+iby*nbeta] << endl; + // if (ipy==3 && iby==10) cout << newhhy[ibx+iby*nbeta] << " " << hiy[iby] << " " << ibx+iby*nbeta << endl; + } + } + } } - - + + } public: - etaInterpolationAdaptiveBins(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){}; + etaInterpolationAdaptiveBins(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationPosXY(nx,ny,ns, nb, emin,emax){}; - etaInterpolationAdaptiveBins(etaInterpolationAdaptiveBins *orig): etaInterpolationBase(orig){}; + etaInterpolationAdaptiveBins(etaInterpolationAdaptiveBins *orig): etaInterpolationPosXY(orig){}; virtual etaInterpolationAdaptiveBins* Clone() { @@ -103,7 +125,7 @@ class etaInterpolationAdaptiveBins : public etaInterpolationBase{ virtual void prepareInterpolation(int &ok) { ok=1; - + cout << "Adaptive bins" << endl; ///*Eta Distribution Rebinning*/// double bsize=1./nSubPixels; //precision @@ -113,7 +135,6 @@ class etaInterpolationAdaptiveBins : public etaInterpolationBase{ double tot_eta_y=0; for (int ip=0; ip1 || etah[ii]<0 ) cout << "***"<< ii << etah[ii] << endl; + } + sprintf(tit,"/scratch/neweta_hhy_%d.tiff",iint); + WriteToTiff(etah, tit, etabins, etabins); +#endif + // if (new_diff + + //#include "moench03T1ZmqData.h" +#ifdef NEWRECEIVER +#include "moench03T1ReceiverDataNew.h" +#endif +#ifdef CSAXS_FP #include "moench03T1ReceiverData.h" +#endif +#ifdef OLDDATA +#include "moench03Ctb10GbT1Data.h" +#endif // #include "interpolatingDetector.h" //#include "etaInterpolationPosXY.h" @@ -24,12 +34,13 @@ using namespace std; int main(int argc, char *argv[]) { - if (argc<7) { - cout << "Usage is " << argv[0] << "indir outdir fname runmin runmax pedfname" << endl; + if (argc<6) { + cout << "Usage is " << argv[0] << "indir outdir fname runmin runmax " << endl; + return 1; } int p=10000; int fifosize=1000; - int nthreads=20; + int nthreads=8; int nsubpix=25; int etabins=nsubpix*10; double etamin=-1, etamax=2; @@ -42,11 +53,30 @@ int main(int argc, char *argv[]) { int ok; int iprog=0; + + + + +#ifdef NEWRECEIVER moench03T1ReceiverDataNew *decoder=new moench03T1ReceiverDataNew(); + cout << "RECEIVER DATA WITH ONE HEADER!"<setFrameMode(eFrame); - // mt->setFrameMode(ePedestal); - cout << pedfname<< endl; + // filter->setFrameMode(eFrame); + // mt->setFrameMode(ePedestal); + // cout << pedfname<< endl; - filebin.open((const char *)(pedfname), ios::in | ios::binary); - //filebin.open((const char *)(fname), ios::in | ios::binary); + // filebin.open((const char *)(pedfname), ios::in | ios::binary); + // //filebin.open((const char *)(fname), ios::in | ios::binary); - // //open file - if (filebin.is_open()){ - // //while read frame - cout << "pedestal file " << endl; - while (decoder->readNextFrame(filebin, ff, np,data)) { - // cout << ff << " " << np << endl; - // //push - // mt->pushData(buff); - // // //pop - //mt->nextThread(); - // // // cout << " " << (void*)buff; - //mt->popFree(buff); - filter->processData(data); - } - filebin.close(); - // //close file - // //join threads - // while (mt->isBusy()) {;}//wait until all data are processed from the queues - // cout << outfname << endl; - // filter->writePedestals(outfname); - // sprintf(outfname,"%s/%s_pedimg.tiff",outdir,fn); + // // //open file + // if (filebin.is_open()){ + // // //while read frame + // cout << "pedestal file " << endl; + // while (decoder->readNextFrame(filebin, ff, np,data)) { + // // cout << ff << " " << np << endl; + // // //push + // // mt->pushData(buff); + // // // //pop + // //mt->nextThread(); + // // // // cout << " " << (void*)buff; + // //mt->popFree(buff); + // filter->processData(data); + // } + // filebin.close(); + // // //close file + // // //join threads + // // while (mt->isBusy()) {;}//wait until all data are processed from the queues + // // cout << outfname << endl; + // // filter->writePedestals(outfname); + // // sprintf(outfname,"%s/%s_pedimg.tiff",outdir,fn); - // cout << outfname << endl; + // // cout << outfname << endl; - // filter->writeImage(outfname); - // //mt->clearImage(); - } else - cout << "Could not open "<< pedfname << " for reading " << endl; + // // filter->writeImage(outfname); + // // //mt->clearImage(); + // } else + // cout << "Could not open "<< pedfname << " for reading " << endl; - // for (int ix=0; ix<400; ix++) - // for (int iy=0; iy<400; iy++) - // cout << ix << " " << iy << " " << filter->getPedestal(ix,iy) << " " << filter->getPedestalRMS(ix,iy) << endl; + // // for (int ix=0; ix<400; ix++) + // // for (int iy=0; iy<400; iy++) + // // cout << ix << " " << iy << " " << filter->getPedestal(ix,iy) << " " << filter->getPedestalRMS(ix,iy) << endl; @@ -149,11 +179,12 @@ filter->setFrameMode(eFrame); // mt->setFrameMode(eFrame); //need to find a way to switch between flat and frames! // mt->prepareInterpolation(ok); - mt->setFrameMode(eFrame); + mt->setFrameMode(eFrame); mt->StartThreads(); mt->popFree(buff); + cout << "mt " << endl; int ifr=0; // //loop on files @@ -165,37 +196,37 @@ filter->setFrameMode(eFrame); - for (int irun=runmin; irunreadNextFrame(filebin, ff, np,buff)) { - // cout << "*"<pushData(buff); - // // //pop - mt->nextThread(); - // // // cout << " " << (void*)buff; - mt->popFree(buff); + // filebin.open((const char *)(fname), ios::in | ios::binary); + // // //open file + // if (filebin.is_open()){ + // // //while read frame + // while (decoder->readNextFrame(filebin, ff, np,buff)) { + // // cout << "*"<pushData(buff); + // // // //pop + // mt->nextThread(); + // // // // cout << " " << (void*)buff; + // mt->popFree(buff); - } - // cout << "--" << endl; - filebin.close(); + // } + // // cout << "--" << endl; + // filebin.close(); - } + // } - } + // } - while (mt->isBusy()) {;}//wait until all data are processed from the queues - mt->clearImage(); + // while (mt->isBusy()) {;}//wait until all data are processed from the queues + // mt->clearImage(); for (int irun=runmin; irunsetFrameMode(eFrame); return 1; } // //while read frame + ff=-1; while (decoder->readNextFrame(filebin, ff, np,buff)) { // cout << "*"<pushData(buff); + mt->pushData(buff); // // //pop - mt->nextThread(); + mt->nextThread(); // // // cout << " " << (void*)buff; - mt->popFree(buff); - + mt->popFree(buff); + ff=-1; } - // cout << "--" << endl; + cout << "--" << endl; filebin.close(); // //close file // //join threads diff --git a/slsDetectorCalibration/moenchExecutables/moench03Interpolation.cpp b/slsDetectorCalibration/moenchExecutables/moench03Interpolation.cpp index d2932bfa8..87ce84af0 100644 --- a/slsDetectorCalibration/moenchExecutables/moench03Interpolation.cpp +++ b/slsDetectorCalibration/moenchExecutables/moench03Interpolation.cpp @@ -5,8 +5,8 @@ #include "moench03T1ZmqData.h" #include "single_photon_hit.h" - #include "etaInterpolationPosXY.h" - +// #include "etaInterpolationPosXY.h" +#include "etaInterpolationAdaptiveBins.h" using namespace std; #define NC 400 #define NR 400 @@ -43,9 +43,13 @@ int main(int argc, char *argv[]) { FILE *f=NULL; single_photon_hit cl(3,3); - etaInterpolationPosXY *interp=new etaInterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax); + // etaInterpolationPosXY *interp=new etaInterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax); + etaInterpolationAdaptiveBins *interp=new etaInterpolationAdaptiveBins (NC, NR, nsubpix, etabins, etamin, etamax); + //#ifndef FF + cout << "read ff " << argv[2] << endl; interp->readFlatField(argv[2]); interp->prepareInterpolation(ok); + //#endif int *img; float *totimg=new float[NC*NR*nsubpix*nsubpix]; @@ -133,7 +137,7 @@ int main(int argc, char *argv[]) { } } - + cout << "writing eta!" << endl; WriteToTiff(ffimg, outfname,NC*nsubpix,NR*nsubpix); diff --git a/slsDetectorCalibration/moenchExecutables/readClusters.C b/slsDetectorCalibration/moenchExecutables/readClusters.C index d75338508..f28cb95e0 100644 --- a/slsDetectorCalibration/moenchExecutables/readClusters.C +++ b/slsDetectorCalibration/moenchExecutables/readClusters.C @@ -1,5 +1,10 @@ #include "../single_photon_hit.h" //#include "etaVEL/etaInterpolationPosXY.h" +#include +#include +#include +#include +using namespace std; TH2F *readClusters(char *fname, int nx, int ny, TH2F *h2=NULL) { FILE *f=fopen(fname,"r"); @@ -13,21 +18,21 @@ TH2F *readClusters(char *fname, int nx, int ny, TH2F *h2=NULL) { h2=new TH2F("h2",fname,nx, -0.5, nx-0.5, ny, -0.5, ny-0.5); //h2mult=new TH2F("h2mult",fname,nx, -0.5, nx-0.5, ny, -0.5, ny-0.5); // TH2F *hint=new TH2F("hint",fname,nx*ns, -0.5, nx-0.5, ny*ns, -0.5, ny-0.5); - TH1F *hf=new TH1F("hf","hf",1000,0,3000000); + TH1F *hf=new TH1F("hf","hf",1000,0,10E6); //TH2F *hff=new TH2F("hff","hff",ns, -0.5, 0.5, ns, -0.5, +0.5); - TH1F *hsp=new TH1F("hsp",fname,500,0,2000); - TH1F *hsp1=new TH1F("hsp1",fname,500,0,2000); + TH1F *hsp=new TH1F("hsp",fname,500,0,10000); + // TH1F *hsp1=new TH1F("hsp1",fname,500,0,10000); // TH1F *hsp2=new TH1F("hsp2",fname,500,0,1000); // TH1F *hsp3=new TH1F("hsp3",fname,500,0,1000); - hsp1->SetLineColor(2); + // hsp1->SetLineColor(2); // hsp2->SetLineColor(3); // hsp3->SetLineColor(4); TCanvas *c=new TCanvas(); - c->SetLogz(kTRUE); h2->Draw("colz"); TCanvas *c1=new TCanvas(); hsp->Draw(); - hsp1->Draw("same"); + c1->SetLogy(kTRUE); + // hsp1->Draw("same"); // hsp2->Draw("same"); // hsp3->Draw("same"); TCanvas *c2=new TCanvas(); @@ -40,6 +45,8 @@ TH2F *readClusters(char *fname, int nx, int ny, TH2F *h2=NULL) { double phw=340, phs=62; int f0=-1; double tl, bl, tr, br, qt; + int iimage=0; + while (cl.read(f)) { //cl.get_pixel(x1, y1); //cout << cl.iframe << " " << cl.x << " " << cl.y << endl; @@ -64,7 +71,6 @@ TH2F *readClusters(char *fname, int nx, int ny, TH2F *h2=NULL) { // if (iy>0) top+=cl.get_data(ix,iy); } - qt=bl; if (br>qt) qt=br; if (tl>qt) qt=tl; @@ -75,22 +81,22 @@ TH2F *readClusters(char *fname, int nx, int ny, TH2F *h2=NULL) { //max at 340 //if (tot>200) { w=1; - + if (qt>1000) { if (qt/tot>0.8 && qt/tot<1.2){ if (f0<0) f0=cl.iframe; hf->Fill(cl.iframe-f0); - if (qt>540) w++; - if (qt>820) w++; + // if (qt>540) w++; + // if (qt>820) w++; //(tot+3.5*phs)/phw; //} else //w=0; // if (w) { // if (cl.y<350) { // if (cl.y<100 || cl.y>300) { - if (cl.x>150 && cl.x<250 && cl.y>200 && cl.y<250) - hsp1->Fill(qt); - else + // if (cl.x>150 && cl.x<250 && cl.y>200 && cl.y<250) + // hsp1->Fill(qt); + // else hsp->Fill(qt); // hsp2->Fill(cl.get_data(0,0)); @@ -122,10 +128,11 @@ TH2F *readClusters(char *fname, int nx, int ny, TH2F *h2=NULL) { c2->Modified();; c2->Update(); } - + // if (iph>1E7) + // break; //} // if (iph>0.5E7) break; - + } } fclose(f); // hff->Scale(hff->GetNbinsX()*hff->GetNbinsY()/hff->Integral()); diff --git a/slsDetectorCalibration/singlePhotonDetector.h b/slsDetectorCalibration/singlePhotonDetector.h index 38074b581..7bc09ebd4 100644 --- a/slsDetectorCalibration/singlePhotonDetector.h +++ b/slsDetectorCalibration/singlePhotonDetector.h @@ -447,7 +447,7 @@ int *getClusters(char *data, int *ph=NULL) { (clusters+nph)->rms=getPedestalRMS(ix,iy); - + cluster=clusters+nph; for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) {