diff --git a/slsDetectorCalibration/analogDetector.h b/slsDetectorCalibration/analogDetector.h index 548a5e181..abd019948 100644 --- a/slsDetectorCalibration/analogDetector.h +++ b/slsDetectorCalibration/analogDetector.h @@ -253,16 +253,54 @@ template class analogDetector { \param iy pixel y coordinate \param cm 1 adds the value to common mod, 0 skips it. Defaults to 0. - not properly implemented */ - virtual void addToPedestal(double val, int ix, int iy=0){ - if (ix>=0 && ix=0 && iy=0 && ix=0 && iy0) { + val-= getCommonMode(ix, iy); + } stat[iy][ix].addToPedestal(val); - if (cmSub) { - if (det) if (det->isGood(ix, iy)==0) return; - cmSub->addToCommonMode(val, ix, iy); - }; + /* if (cmSub && cm>0) { */ + /* if (det) if (det->isGood(ix, iy)==0) return; */ + /* cmSub->addToCommonMode(val, ix, iy); */ + /* }; */ }; } - + + double getCommonMode(int ix, int iy) { + if (cmSub) return cmSub->getCommonMode(ix, iy); + else return 0; + } + + virtual void addToCommonMode(double val, int ix, int iy=0){ + + if (ix>=0 && ix=0 && iyisGood(ix, iy)==0) return; + if (getNumpedestals(ix,iy)>0) + cmSub->addToCommonMode(val-getPedestal(ix,iy), ix, iy); + }; + + } + } + + virtual void addToCommonMode(char *data){ + if (cmSub) { + for (int ix=xmin; ix0) + addToCommonMode(data, ix, iy); + } + } + cout << "cm " << getCommonMode(0,0) << " " << getCommonMode(1,0) << endl; + } + } + virtual void addToCommonMode(char *data, int ix, int iy=0){ + if (cmSub) { + if (det) if (det->isGood(ix, iy)==0) return; + if (getNumpedestals(ix,iy)>0) + cmSub->addToCommonMode(subtractPedestal(data,ix,iy,0), ix, iy); + } + } /** gets pedestal (and common mode) \param ix pixel x coordinate @@ -270,7 +308,13 @@ template class analogDetector { \param cm 0 (default) without common mode subtraction, 1 with common mode subtraction (if defined) \returns pedestal value */ - virtual double getPedestal(int ix, int iy, int cm=0){if (ix>=0 && ix=0 && iy0) return stat[iy][ix].getPedestal()-cmSub->getCommonMode(); else return stat[iy][ix].getPedestal(); else return -1;}; + virtual double getPedestal(int ix, int iy, int cm=0){ + if (ix>=0 && ix=0 && iy0) + return stat[iy][ix].getPedestal()+getCommonMode(ix,iy); + else return stat[iy][ix].getPedestal(); + else return -1; + }; /** gets pedestal rms (i.e. noise) @@ -278,9 +322,17 @@ template class analogDetector { \param iy pixel y coordinate \returns pedestal rms */ - virtual double getPedestalRMS(int ix, int iy){if (ix>=0 && ix=0 && iy=0 && ix=0 && iy=0 && ix=0 && iy class analogDetector { gm=new float[nx*ny]; for (int ix=0; ixgetCommonMode(); - else gm[iy*nx+ix]=stat[iy][ix].getPedestal(); } } @@ -538,17 +587,18 @@ template class analogDetector { \param data pointer to the data */ - virtual void addToPedestal(char *data) { + virtual void addToPedestal(char *data, int cm=0) { newFrame(); - - + if (cmSub) { + addToCommonMode(data); + } + for (int ix=xmin; ix class analogDetector { */ - virtual void addToPedestal(char *data, int ix, int iy=0) { + virtual void addToPedestal(char *data, int ix, int iy=0, int cm=0) { double val; @@ -618,7 +668,8 @@ template class analogDetector { val=dataSign*det->getValue(data, ix, iy); else val=((double*)data)[iy*nx+ix]; - + if (cm && cmSub) + val-=getCommonMode(ix,iy); addToPedestal(val,ix,iy); } return ; @@ -634,7 +685,7 @@ template class analogDetector { */ - virtual double *subtractPedestal(char *data, double *val=NULL) { + virtual double *subtractPedestal(char *data, double *val=NULL, int cm=0) { newFrame(); @@ -643,7 +694,7 @@ template class analogDetector { for (int ix=xmin; ix class analogDetector { - virtual double subtractPedestal(char *data, int ix, int iy=0) { + virtual double subtractPedestal(char *data, int ix, int iy=0, int cm=0) { double g=1.; if (ix>=0 && ix=0 && iygetValue(data, ix, iy)-getPedestal(ix,iy))/g; + + if (det) { + // cout << det->getValue(data, ix, iy) << " " << getPedestal(ix,iy) << " " << (dataSign*det->getValue(data, ix, iy)-getPedestal(ix,iy))/g << endl; + return (dataSign*det->getValue(data, ix, iy)-getPedestal(ix,iy, cm))/g; + } else - return (((double*)data)[iy*nx+ix]-getPedestal(ix,iy))/g; + return (((double*)data)[iy*nx+ix]-getPedestal(ix,iy,cm))/g; } }; @@ -700,15 +753,22 @@ template class analogDetector { virtual int getNPhotons(char *data, int ix, int iy=0) { int nph=0; - double v; + double v; + int cm=0; + if (cmSub) cm=1; if (ix>=0 && ix=0 && iy0) { v+=0.5*thr; nph=v/thr; return nph; - } else - return v; + } else { + // cout << v << endl; + //if (v>0) + return v; + //else + //return 0; + } } return 0; }; @@ -721,10 +781,14 @@ template class analogDetector { */ int *getNPhotons(char *data, int *nph=NULL) { - double val; + // double val; if (nph==NULL) nph=image; newFrame(); + + + addToCommonMode(data); + for (int ix=xmin; ix *filter=new analogDetector(det, 1, cm, 1000); filter->setROI(0,2560,0,1); char *buff;//[2*(48+1280*2)]; char *buff0; char *buff1; multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize); - mt->setFrameMode(eFrame); + // mt->setFrameMode(eFrame); + mt->setFrameMode(ePedestal); mt->StartThreads(); mt->popFree(buff); buff0=buff; @@ -297,39 +289,44 @@ int main(int argc, char *argv[]){ 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<2560; i++) { - // // fprintf(fout,"%d %d\n",i,image[i]); - dout[i]=image[i]; - } - // // fclose(fout);; - } + // while (mt->isBusy()) {;} + // image=filter->getImage(); + // if (image) { + // 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]; + // if (image[i]<0) + // dout[i]=0; + // } + // fclose(fout); + // } - if (send) { - strcpy(fname0,filename0.c_str()); - strcpy(fname1,filename1.c_str()); - // zmqsocket2->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,16,fileindex,400,400,400*400, acqIndex,frameIndex,fname, acqIndex, 0,0,0,0,0,0,0,0,0,0,0,1); - zmqsocket2->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,0,0,0,0,0, 0,0,fname0, 0, 0,0,0,0,0,0,0,0,0,0,0,1); - zmqsocket3->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,0,0,0,0,0, 0,0,fname1, 0, 0,0,0,0,0,0,0,0,0,0,0,1); + if (send) { + // strcpy(fname0,filename0.c_str()); + // strcpy(fname1,filename1.c_str()); + // // zmqsocket2->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,16,fileindex,400,400,400*400, acqIndex,frameIndex,fname, acqIndex, 0,0,0,0,0,0,0,0,0,0,0,1); + // zmqsocket2->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,0,0,0,0,0, 0,0,fname0, 0, 0,0,0,0,0,0,0,0,0,0,0,1); + // zmqsocket3->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,0,0,0,0,0, 0,0,fname1, 0, 0,0,0,0,0,0,0,0,0,0,0,1); - zmqsocket2->SendData((char*)dout,size/2); - zmqsocket3->SendData(((char*)dout)+size/2,size/2); - // cprintf(GREEN, "Sent Data\n"); + // zmqsocket2->SendData((char*)dout,size/2); + // zmqsocket3->SendData(((char*)dout)+size/2,size/2); + // // cprintf(GREEN, "Sent Data\n"); - zmqsocket2->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION); - zmqsocket3->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION); + // zmqsocket2->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION); + // zmqsocket3->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION); - cprintf(RED, "Received %d frames\n", nf); + // cprintf(RED, "Received %d frames\n", nf); - } + zmqsocket2->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION); + zmqsocket3->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION); + } + mt->setFrameMode(eFrame); filter->clearImage(); // std::time(&end_time); // cout << std::ctime(&end_time) << " " << nf << endl; @@ -352,16 +349,61 @@ int main(int argc, char *argv[]){ irun=fileindex0; sprintf(ofname,"%s_%d.ph",filename0.c_str(),fileindex0); + + while (mt->isBusy()) {;} + image=filter->getImage(); + if (image) { + for (int i=0; i<2560; i++) { + // fprintf(fout,"%d %d\n",i,image[i]); + dout[i]=det->getChannel(buff,i,0);//image[i]+1000;//filter->getPedestal(i,0);// + // if (image[i]<0) + // dout[i]=0; + // cout << i << " " << image[i] << " " << dout[i] << endl; + } + } + + + + + + if (send) { + strcpy(fname0,filename0.c_str()); + strcpy(fname1,filename1.c_str()); + // zmqsocket2->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,16,fileindex,400,400,400*400, acqIndex,frameIndex,fname, acqIndex, 0,0,0,0,0,0,0,0,0,0,0,1); + zmqsocket2->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,0,0,0,0,0, 0,0,fname0, 0, 0,0,0,0,0,0,0,0,0,0,0,1); + zmqsocket3->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,0,0,0,0,0, 0,0,fname1, 0, 0,0,0,0,0,0,0,0,0,0,0,1); + + zmqsocket2->SendData((char*)dout,size/2); + zmqsocket3->SendData(((char*)dout)+size/2,size/2); + // cprintf(GREEN, "Sent Data\n"); + + + // zmqsocket2->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION); + // zmqsocket3->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION); + + // cprintf(RED, "Received %d frames\n", nf); + + } + + + if (nf>100) + mt->setFrameMode(eFrame); + filter->clearImage(); + + #endif - - + mt->pushData(buff); mt->nextThread(); + + mt->popFree(buff); buff0=buff; buff1=buff+offset*2+1280*2; + + nf++; diff --git a/slsDetectorCalibration/interpolations/etaInterpolationAdaptiveBins.h b/slsDetectorCalibration/interpolations/etaInterpolationAdaptiveBins.h index 2100380f5..3a785f8f2 100644 --- a/slsDetectorCalibration/interpolations/etaInterpolationAdaptiveBins.h +++ b/slsDetectorCalibration/interpolations/etaInterpolationAdaptiveBins.h @@ -181,53 +181,73 @@ class etaInterpolationAdaptiveBins : public etaInterpolationPosXY { double avg=tot_eta/((double)(nSubPixels*nSubPixels)); cout << "total eta entries is :"<< tot_eta << " avg: "<< avg << endl; cout << "Start " << endl; - double old_diff=calcDiff(avg, hhx, hhy), new_diff=old_diff+1; - cout << " diff= " << new_diff << endl; + double old_diff=calcDiff(avg, hhx, hhy), new_diff=old_diff+1, best_diff=old_diff; + cout << " diff= " << old_diff << endl; int iint=0; float *newhhx=new float[nbeta*nbeta]; //profile x float *newhhy=new float[nbeta*nbeta]; //profile y + float *besthhx=hhx; //profile x + float *besthhy=hhy; //profile y while (iint1 || etah[ii]<0 ) cout << "***"<< ii << etah[ii] << endl; */ - } - sprintf(tit,"/scratch/neweta_hhx_%d.tiff",iint); - WriteToTiff(etah, tit, etabins, etabins); +/* } */ +/* sprintf(tit,"/scratch/neweta_hhx_%d.tiff",iint); */ +/* WriteToTiff(etah, tit, etabins, etabins); */ - for (int ii=0; ii1 || etah[ii]<0 ) cout << "***"<< ii << etah[ii] << endl; +/* for (int ii=0; ii1 || 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=etamax) { + cout << "aaa:" <2) - calcEta(totquad, sDum, etax, etay); - getInterpolatedPosition(x,y,etax,etay,corner,int_x,int_y); + /* virtual void getInterpolatedPosition(int x, int y, int *data, double &int_x, double &int_y)=0; */ + /* virtual void getInterpolatedPosition(int x, int y, double *data, double &int_x, double &int_y)=0; */ + /* virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cl,double &int_x, double &int_y)=0; */ + /* virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,int *cl,double &int_x, double &int_y)=0; */ + /* virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int corner, double &int_x, double &int_y)=0; */ - return; - }; + + /* virtual int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay)=0; */ + /* virtual int addToFlatField(double totquad,int quad,double *cl,double &etax, double &etay)=0; */ + /* virtual int addToFlatField(double *cluster, double &etax, double &etay)=0; */ + /* virtual int addToFlatField(int *cluster, double &etax, double &etay)=0; */ - virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cl,double &int_x, double &int_y) { - - double cc[2][2]; - double *cluster[3]; - int xoff, yoff; - cluster[0]=cl; - cluster[1]=cl+3; - cluster[2]=cl+6; - - switch (quad) { - case BOTTOM_LEFT: - xoff=0; - yoff=0; - break; - case BOTTOM_RIGHT: - xoff=1; - yoff=0; - break; - case TOP_LEFT: - xoff=0; - yoff=1; - break; - case TOP_RIGHT: - xoff=1; - yoff=1; - break; - default: - ; - } - double etax, etay; - if (nSubPixels>2) { - cc[0][0]=cluster[yoff][xoff]; - cc[1][0]=cluster[yoff+1][xoff]; - cc[0][1]=cluster[yoff][xoff+1]; - cc[1][1]=cluster[yoff+1][xoff+1]; - calcEta(totquad,cc,etax,etay); - } - return getInterpolatedPosition(x,y,etax, etay,quad,int_x,int_y); - - } - - - - - - - virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int corner, double &int_x, double &int_y) - { - - - double xpos_eta=0,ypos_eta=0; - double dX,dY; - int ex,ey; - switch (corner) - { - case TOP_LEFT: - dX=-1.; - dY=0; - break; - case TOP_RIGHT: - ; - dX=0; - dY=0; - break; - case BOTTOM_LEFT: - dX=-1.; - dY=-1.; - break; - case BOTTOM_RIGHT: - dX=0; - dY=-1.; - break; - default: - cout << "bad quadrant" << endl; - dX=0.; - dY=0.; - } - - - if (nSubPixels>2) { - -#ifdef MYROOT1 - xpos_eta=(hhx->GetBinContent(hhx->GetXaxis()->FindBin(etax),hhy->GetYaxis()->FindBin(etay)))/((double)nSubPixels); - ypos_eta=(hhy->GetBinContent(hhx->GetXaxis()->FindBin(etax),hhy->GetYaxis()->FindBin(etay)))/((double)nSubPixels); -#endif -#ifndef MYROOT1 - ex=(etax-etamin)/etastep; - ey=(etay-etamin)/etastep; - if (ex<0) { - ex=0; - cout << "x*"<< ex << endl; - } - if (ex>=nbeta) { - ex=nbeta-1; - cout << "x?"<< ex << endl; - - } - if (ey<0) { - ey=0; - cout << "y*"<< ey << endl; - } - if (ey>=nbeta) { - ey=nbeta-1; - cout << "y?"<< ey << endl; - - } - - - - xpos_eta=(((double)hhx[(ey*nbeta+ex)]))+dX ;///((double)nSubPixels); - ypos_eta=(((double)hhy[(ey*nbeta+ex)]))+dY ;///((double)nSubPixels); - //else - //return 0; - -#endif - } else { - xpos_eta=0.5*dX+0.25; - ypos_eta=0.5*dY+0.25; - } - - int_x=((double)x) + xpos_eta+0.5; - int_y=((double)y) + ypos_eta+0.5; - - - } - - - - - - - ///////////////////////////////////////////////////////////////////////////////////////////////// - virtual void getPositionETA3(int x, int y, double *data, double &int_x, double &int_y) - { - double sDum[2][2]; - double tot, totquad; - double eta3x,eta3y; - double ex,ey; - - calcQuad(data, tot, totquad, sDum); - calcEta3(data,eta3x, eta3y,tot); - - double xpos_eta,ypos_eta; - -#ifdef MYROOT1 - xpos_eta=((hhx->GetBinContent(hhx->GetXaxis()->FindBin(eta3x),hhy->GetYaxis()->FindBin(eta3y))))/((double)nSubPixels); - ypos_eta=((hhy->GetBinContent(hhx->GetXaxis()->FindBin(eta3x),hhy->GetYaxis()->FindBin(eta3y))))/((double)nSubPixels); - -#endif -#ifndef MYROOT1 - ex=(eta3x-etamin)/etastep; - ey=(eta3y-etamin)/etastep; - - if (ex<0) ex=0; - if (ex>=nbeta) ex=nbeta-1; - if (ey<0) ey=0; - if (ey>=nbeta) ey=nbeta-1; - - xpos_eta=(((double)hhx[(int)(ey*nbeta+ex)]))/((double)nSubPixels); - ypos_eta=(((double)hhy[(int)(ey*nbeta+ex)]))/((double)nSubPixels); -#endif - - int_x=((double)x) + xpos_eta; - int_y=((double)y) + ypos_eta; - - return; - }; - - virtual int addToFlatField(double totquad,int quad,double *cl,double &etax, double &etay) { - double cc[2][2]; - double *cluster[3]; - int xoff, yoff; - cluster[0]=cl; - cluster[1]=cl+3; - cluster[2]=cl+6; - - switch (quad) { - case BOTTOM_LEFT: - xoff=0; - yoff=0; - break; - case BOTTOM_RIGHT: - xoff=1; - yoff=0; - break; - case TOP_LEFT: - xoff=0; - yoff=1; - break; - case TOP_RIGHT: - xoff=1; - yoff=1; - break; - default: - ; - } - cc[0][0]=cluster[yoff][xoff]; - cc[1][0]=cluster[yoff+1][xoff]; - cc[0][1]=cluster[yoff][xoff+1]; - cc[1][1]=cluster[yoff+1][xoff+1]; - - /* cout << cl[0] << " " << cl[1] << " " << cl[2] << endl; */ - /* cout << cl[3] << " " << cl[4] << " " << cl[5] << endl; */ - /* cout << cl[6] << " " << cl[7] << " " << cl[8] << endl; */ - /* cout <<"******"<Fill(etax,etay); #endif #ifndef MYROOT1 + int ex,ey; ex=(etax-etamin)/etastep; ey=(etay-etamin)/etastep; - // cout << etax << " " << ex << " " << etay << " " << ey << " " << ey*nbeta+ex << endl; if (ey=0 && ey>=0) heta[ey*nbeta+ex]++; - // cout << "*"<< etax << " " << etay << endl; - /* cout << etax << " " << etay << " " << ex << " " << ey << " " << ey*nbeta+ex << endl; */ - /* cout <<"********"<< endl << endl ; */ #endif return 0; -}; + }; + + // virtual void prepareInterpolation(int &ok)=0; + protected: #ifdef MYROOT1 diff --git a/slsDetectorCalibration/interpolations/etaInterpolationPosXY.h b/slsDetectorCalibration/interpolations/etaInterpolationPosXY.h index 5417f165c..e730d2df6 100644 --- a/slsDetectorCalibration/interpolations/etaInterpolationPosXY.h +++ b/slsDetectorCalibration/interpolations/etaInterpolationPosXY.h @@ -2,24 +2,27 @@ #define ETA_INTERPOLATION_POSXY_H -#include "tiffIO.h" +//#include "tiffIO.h" #include "etaInterpolationBase.h" +#include "eta2InterpolationBase.h" +#include "eta3InterpolationBase.h" -class etaInterpolationPosXY : public etaInterpolationBase{ +class etaInterpolationPosXY : public virtual etaInterpolationBase{ public: - etaInterpolationPosXY(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){}; + etaInterpolationPosXY(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 << "epxy " << nb << " " << emin << " " << emax << endl; cout << nbeta << " " << etamin << " " << etamax << endl; + }; - etaInterpolationPosXY(etaInterpolationPosXY *orig): etaInterpolationBase(orig){}; + etaInterpolationPosXY(etaInterpolationPosXY *orig): etaInterpolationBase(orig) {}; - virtual etaInterpolationPosXY* Clone() { + virtual etaInterpolationPosXY* Clone()=0;/** { return new etaInterpolationPosXY(this); - }; + };*/ virtual void prepareInterpolation(int &ok) { - cout <<"?"<< endl; ok=1; #ifdef MYROOT1 if (hhx) delete hhx; @@ -114,10 +117,6 @@ class etaInterpolationPosXY : public etaInterpolationBase{ // cout << "y " << nbeta << " " << (ii+1)*tot_eta_x*bsize << " " << ii << endl; - - - - } #ifdef SAVE_ALL @@ -150,4 +149,25 @@ class etaInterpolationPosXY : public etaInterpolationBase{ }; +class eta2InterpolationPosXY : public virtual eta2InterpolationBase, public virtual etaInterpolationPosXY { + public: + 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){ + cout << "e3pxy " << nbeta << " " << etamin << " " << etamax << " " << nSubPixels<< endl; + }; + + eta3InterpolationPosXY(eta3InterpolationPosXY *orig): etaInterpolationBase(orig), etaInterpolationPosXY(orig) {}; + + virtual eta3InterpolationPosXY* Clone() { return new eta3InterpolationPosXY(this);}; +}; + #endif diff --git a/slsDetectorCalibration/interpolations/noInterpolation.h b/slsDetectorCalibration/interpolations/noInterpolation.h index c5036a9f4..6c199399f 100644 --- a/slsDetectorCalibration/interpolations/noInterpolation.h +++ b/slsDetectorCalibration/interpolations/noInterpolation.h @@ -38,13 +38,31 @@ class noInterpolation : public slsInterpolation{ return ; }; + + virtual void getInterpolatedPosition(int x, int y, int *data, double &int_x, double &int_y) + { + return getInterpolatedPosition(x, y, (double*)NULL, int_x, int_y); + } + + virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int corner, double &int_x, double &int_y) { - getInterpolatedPosition(x, y, NULL, int_x, int_y); + getInterpolatedPosition(x, y, (double*)NULL, int_x, int_y); }; + + + + virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cl,double &etax, double &etay){ - getInterpolatedPosition(x, y, NULL, etax, etay); + getInterpolatedPosition(x, y, (double*)NULL, etax, etay); }; + + + virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,int *cl,double &etax, double &etay){ + getInterpolatedPosition(x, y, (double*)NULL, etax, etay); + }; + + ////////////////////////////////////////////////////////////////////////////////////// virtual void getPositionETA3(int x, int y, double *data, double &int_x, double &int_y) @@ -55,13 +73,27 @@ class noInterpolation : public slsInterpolation{ return ; }; + + virtual void getPositionETA3(int x, int y, int *data, double &int_x, double &int_y) + { + return getPositionETA3(x, y, (double*)NULL, int_x, int_y); + }; + + + ////////////////////////////////////////////////////////////////////////////////////// virtual int addToFlatField(double *cluster, double &etax, double &etay){return 0;}; + + virtual int addToFlatField(int *cluster, double &etax, double &etay){return 0;}; + virtual int addToFlatField(double etax, double etay){return 0;}; virtual int addToFlatField(double totquad,int quad,double *cl,double &etax, double &etay){return 0;}; + virtual int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay){return 0;}; + + protected: ; // TRandom *eventGenerator; diff --git a/slsDetectorCalibration/interpolations/slsInterpolation.h b/slsDetectorCalibration/interpolations/slsInterpolation.h index 559f33597..ea79eb639 100644 --- a/slsDetectorCalibration/interpolations/slsInterpolation.h +++ b/slsDetectorCalibration/interpolations/slsInterpolation.h @@ -8,7 +8,10 @@ #endif #include +#ifndef MY_TIFF_IO_H #include "tiffIO.h" +#endif + #ifndef DEF_QUAD #define DEF_QUAD enum quadrant { @@ -88,7 +91,7 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny); void *writeInterpolatedImage(const char * imgname) { - cout << "!" <=0 && ix<(nPixelsX*nSubPixels) && iy<(nSubPixels*nPixelsY) && iy>=0 )(*(hint+ix+iy*nPixelsX*nSubPixels))+=1; return hint; @@ -147,9 +154,10 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny); virtual int addToFlatField(double *cluster, double &etax, double &etay)=0; - virtual int addToFlatField(double etax, double etay)=0; - + virtual int addToFlatField(int *cluster, double &etax, double &etay)=0; + virtual int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay)=0; virtual int addToFlatField(double totquad,int quad,double *cluster,double &etax, double &etay)=0; + virtual int addToFlatField(double etax, double etay)=0; #ifdef MYROOT1 virtual TH2D *getFlatField(){return NULL;}; @@ -167,6 +175,14 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny); //virtual void Streamer(TBuffer &b); + static int calcQuad(int *cl, double &sum, double &totquad, double sDum[2][2]){ + double cli[3*3];//=new int[3*3]; + for (int i=0; i<9; i++) + cli[i]=cl[i]; + return calcQuad(cli, sum, totquad, sDum); + + } + static int calcQuad(double *cl, double &sum, double &totquad, double sDum[2][2]){ @@ -175,48 +191,64 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny); cluster[0]=cl; cluster[1]=cl+3; cluster[2]=cl+6; + + sum=0; + double sumBL=0; + double sumTL=0; + double sumBR=0; + double sumTR=0; + int xoff=0, yoff=0; + for (int ix=0; ix<3; ix++) { + for (int iy=0; iy<3; iy++) { + sum+=cluster[iy][ix]; + if (ix<=1 && iy<=1) sumBL+=cluster[iy][ix]; + if (ix<=1 && iy>=1) sumTL+=cluster[iy][ix]; + if (ix>=1 && iy<=1) sumBR+=cluster[iy][ix]; + if (ix>=1 && iy>=1) sumTR+=cluster[iy][ix]; + } + } - sum = cluster[0][0] + cluster[1][0] + cluster[2][0] + cluster[0][1] + cluster[1][1] + cluster[2][1] + cluster[0][2] + cluster[1][2] + cluster[2][2]; - - double sumBL = cluster[0][0] + cluster[1][0] + cluster[0][1] + cluster[1][1]; //2 ->BL - double sumTL = cluster[1][0] + cluster[2][0] + cluster[2][1] + cluster[1][1]; //0 ->TL - double sumBR = cluster[0][1] + cluster[0][2] + cluster[1][2] + cluster[1][1]; //3 ->BR - double sumTR = cluster[1][2] + cluster[2][1] + cluster[2][2] + cluster[1][1]; //1 ->TR - double sumMax = 0; - double t, r; - - // if(sumTL >= sumMax){ - sDum[0][0] = cluster[0][0]; sDum[1][0] = cluster[1][0]; - sDum[0][1] = cluster[0][1]; sDum[1][1] = cluster[1][1]; + /* sDum[0][0] = cluster[0][0]; sDum[1][0] = cluster[1][0]; */ + /* sDum[0][1] = cluster[0][1]; sDum[1][1] = cluster[1][1]; */ corner = BOTTOM_LEFT; - sumMax=sumBL; - // } - - if(sumTL >= sumMax){ - sDum[0][0] = cluster[1][0]; sDum[1][0] = cluster[2][0]; - sDum[0][1] = cluster[1][1]; sDum[1][1] = cluster[2][1]; + totquad=sumBL; + + if(sumTL >= totquad){ + /* sDum[0][0] = cluster[1][0]; sDum[1][0] = cluster[2][0]; */ + /* sDum[0][1] = cluster[1][1]; sDum[1][1] = cluster[2][1]; */ + corner = TOP_LEFT; - sumMax=sumTL; + totquad=sumTL; + xoff=0; + yoff=1; } - if(sumBR >= sumMax){ - sDum[0][0] = cluster[0][1]; sDum[1][0] = cluster[1][1]; - sDum[0][1] = cluster[0][2]; sDum[1][1] = cluster[1][2]; + if(sumBR >= totquad){ + /* sDum[0][0] = cluster[0][1]; sDum[1][0] = cluster[1][1]; */ + /* sDum[0][1] = cluster[0][2]; sDum[1][1] = cluster[1][2]; */ + xoff=1; + yoff=0; corner = BOTTOM_RIGHT; - sumMax=sumBR; + totquad=sumBR; } - if(sumTR >= sumMax){ - sDum[0][0] = cluster[1][1]; sDum[1][0] = cluster[2][1]; - sDum[0][1] = cluster[1][2]; sDum[1][1] = cluster[2][2]; - + if(sumTR >= totquad){ + xoff=1; + yoff=1; + /* sDum[0][0] = cluster[1][1]; sDum[1][0] = cluster[2][1]; */ + /* sDum[0][1] = cluster[1][2]; sDum[1][1] = cluster[2][2]; */ corner = TOP_RIGHT; - sumMax=sumTR; + totquad=sumTR; } - totquad=sumMax; + + for (int ix=0; ix<2; ix++) { + for (int iy=0; iy<2; iy++) { + sDum[iy][ix] = cluster[iy+yoff][ix+xoff]; + } + } return corner; @@ -235,7 +267,6 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny); } - static int calcEta(double *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]) { int corner = calcQuad(cl,sum,totquad,sDum); calcEta(totquad, sDum, etax, etay); @@ -244,6 +275,14 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny); } + static int calcEta(int *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]) { + int corner = calcQuad(cl,sum,totquad,sDum); + calcEta(totquad, sDum, etax, etay); + + return corner; + } + + static int calcEtaL(double totquad, int corner, double sDum[2][2], double &etax, double &etay){ double t,r, toth, totv; if (totquad>0) { @@ -294,6 +333,13 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny); return corner; } + static int calcEtaL(int *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]) { + int corner = calcQuad(cl,sum,totquad,sDum); + calcEtaL(totquad, corner, sDum, etax, etay); + + return corner; + } + static int calcEtaC3(double *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]){ @@ -306,23 +352,95 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny); + static int calcEtaC3(int *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]){ + + int corner = calcQuad(cl,sum,totquad,sDum); + calcEta(sum, sDum, etax, etay); + return corner; + + } + + + static int calcEta3(double *cl, double &etax, double &etay, double &sum) { - double l,r,t,b; - sum=cl[0]+cl[1]+cl[2]+cl[3]+cl[4]+cl[5]+cl[6]+cl[7]+cl[8]; + double l=0,r=0,t=0,b=0, val; + sum=0; + // int quad; + for (int ix=0; ix<3; ix++) { + for (int iy=0; iy<3; iy++) { + val=cl[iy+3*ix]; + sum+=val; + if (iy==0) l+=val; + if (iy==2) r+=val; + if (ix==0) b+=val; + if (ix==2) t+=val; + } + } if (sum>0) { - l=cl[0]+cl[3]+cl[6]; - r=cl[2]+cl[5]+cl[8]; - b=cl[0]+cl[1]+cl[2]; - t=cl[6]+cl[7]+cl[8]; etax=(-l+r)/sum; etay=(-b+t)/sum; } + /* if (etax<-1 || etax>1 || etay<-1 || etay>1) { */ + /* cout << "**********" << etax << " " << etay << endl; */ + /* for (int ix=0; ix<3; ix++) { */ + /* for (int iy=0; iy<3; iy++) { */ + /* cout << cl[iy+3*ix] << "\t" ; */ + + /* } */ + /* cout << endl; */ + /* } */ + /* cout << sum << " " << l << " " << r << " " << t << " " << b << endl; */ + + /* } */ + + + if (etax>=0 && etay>=0) + return TOP_RIGHT; + if (etax<0 && etay>=0) + return TOP_LEFT; + if (etax<0 && etay<0) + return BOTTOM_LEFT; + return BOTTOM_RIGHT; + } + + + static int calcEta3(int *cl, double &etax, double &etay, double &sum) { + double cli[9]; + for (int ix=0; ix<9; ix++) cli[ix]=cl[ix]; + + return calcEta3(cli, etax, etay, sum); + } + + + static int calcMyEta(double totquad, int quad, double *cl, double &etax, double &etay) { + double l,r,t,b, sum; + int yoff; + switch (quad) { + case BOTTOM_LEFT: + case BOTTOM_RIGHT: + yoff=0; + break; + case TOP_LEFT: + case TOP_RIGHT: + yoff=1; + break; + default: + ; + } + l=cl[0+yoff*3]+cl[0+yoff*3+3]; + r=cl[2+yoff*3]+cl[2+yoff*3+3]; + b=cl[0+yoff*3]+cl[1+yoff*3]*cl[2+yoff*3]; + t=cl[0+yoff*3+3]+cl[1+yoff*3+3]*cl[0+yoff*3+3]; + sum=t+b; + if (sum>0) { + etax=(-l+r)/sum; + etay=(+t)/sum; + } return -1; } - - static int calcMyEta(double totquad, int quad, double *cl, double &etax, double &etay) { + static int calcMyEta(double totquad, int quad, int *cl, double &etax, double &etay) { double l,r,t,b, sum; int yoff; switch (quad) { @@ -367,6 +485,21 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny); } + static int calcEta3X(int *cl, double &etax, double &etay, double &sum) { + double l,r,t,b; + sum=cl[0]+cl[1]+cl[2]+cl[3]+cl[4]+cl[5]+cl[6]+cl[7]+cl[8]; + if (sum>0) { + l=cl[3]; + r=cl[5]; + b=cl[1]; + t=cl[7]; + etax=(-l+r)/sum; + etay=(-b+t)/sum; + } + return -1; + } + + diff --git a/slsDetectorCalibration/moenchExecutables/Makefile.cluster_finder b/slsDetectorCalibration/moenchExecutables/Makefile.cluster_finder index 9b4c27588..b8414b77f 100644 --- a/slsDetectorCalibration/moenchExecutables/Makefile.cluster_finder +++ b/slsDetectorCalibration/moenchExecutables/Makefile.cluster_finder @@ -11,23 +11,29 @@ MAIN=moench03ClusterFinder.cpp #-lhdf5 #DESTDIR?=../bin -all: moenchClusterFinder moenchMakeEta moenchInterpolation moenchNoInterpolation +all: moenchClusterFinder moenchMakeEta moenchInterpolation moenchNoInterpolation moenchPhotonCounter moenchAnalog moenchClusterFinder: $(MAIN) $(INCS) clean - g++ -o moenchClusterFinder $(MAIN) $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL -DOLDDATA + g++ -o moenchClusterFinder $(MAIN) $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL -DOLDDATA #-DNEWRECEIVER moenchMakeEta: moench03MakeEta.cpp $(INCS) clean g++ -o moenchMakeEta moench03MakeEta.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL moenchInterpolation: moench03Interpolation.cpp $(INCS) clean - g++ -o moenchInterpolation moench03Interpolation.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL + g++ -o moenchInterpolation moench03Interpolation.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL moenchNoInterpolation: moench03NoInterpolation.cpp $(INCS) clean g++ -o moenchNoInterpolation moench03NoInterpolation.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL +moenchPhotonCounter: moenchPhotonCounter.cpp $(INCS) clean + g++ -o moenchPhotonCounter moenchPhotonCounter.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL -DNEWRECEIVER + +moenchAnalog: moenchPhotonCounter.cpp $(INCS) clean + g++ -o moenchAnalog moenchPhotonCounter.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL -DNEWRECEIVER -DANALOG + clean: - rm -f moenchClusterFinder moenchMakeEta moenchInterpolation moenchNoInterpolation + rm -f moenchClusterFinder moenchMakeEta moenchInterpolation moenchNoInterpolation moenchPhotonCounter diff --git a/slsDetectorCalibration/moenchExecutables/moench03ClusterFinder.cpp b/slsDetectorCalibration/moenchExecutables/moench03ClusterFinder.cpp index f6f8d470c..8a58566db 100644 --- a/slsDetectorCalibration/moenchExecutables/moench03ClusterFinder.cpp +++ b/slsDetectorCalibration/moenchExecutables/moench03ClusterFinder.cpp @@ -119,53 +119,11 @@ int main(int argc, char *argv[]) { cout << "input directory is " << indir << endl; cout << "output directory is " << outdir << endl; cout << "fileformat is " << fformat << endl; - // cout << "pedestal file is " << fformat << endl; std::time(&end_time); cout << std::ctime(&end_time) << 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); - - // // //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; - - // // 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; - + @@ -176,9 +134,7 @@ int main(int argc, char *argv[]) { char* buff; multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize); - - // mt->setFrameMode(eFrame); //need to find a way to switch between flat and frames! - // mt->prepareInterpolation(ok); + mt->setFrameMode(eFrame); mt->StartThreads(); mt->popFree(buff); @@ -187,46 +143,7 @@ int main(int argc, char *argv[]) { cout << "mt " << endl; int ifr=0; - // //loop on files - // mt->setFrameMode(eFrame); - //mt->setFrameMode(eFlat); - - - - - - - // for (int irun=runmin; irunreadNextFrame(filebin, ff, np,buff)) { - // // cout << "*"<pushData(buff); - // // // //pop - // mt->nextThread(); - // // // // cout << " " << (void*)buff; - // mt->popFree(buff); - - // } - // // cout << "--" << endl; - // filebin.close(); - - - // } - - // } - - // while (mt->isBusy()) {;}//wait until all data are processed from the queues - // mt->clearImage(); + for (int irun=runmin; irunnextThread(); // // // cout << " " << (void*)buff; mt->popFree(buff); + ifr++; + if (ifr%10000==0) cout << ifr << " " << ff << endl; ff=-1; } cout << "--" << endl; diff --git a/slsDetectorCalibration/moenchExecutables/moench03Interpolation.cpp b/slsDetectorCalibration/moenchExecutables/moench03Interpolation.cpp index 87ce84af0..2d852a164 100644 --- a/slsDetectorCalibration/moenchExecutables/moench03Interpolation.cpp +++ b/slsDetectorCalibration/moenchExecutables/moench03Interpolation.cpp @@ -2,53 +2,101 @@ #include "ansi.h" #include -#include "moench03T1ZmqData.h" -#include "single_photon_hit.h" +//#include "moench03T1ZmqData.h" +//#define DOUBLE_SPH +//#define MANYFILES -// #include "etaInterpolationPosXY.h" -#include "etaInterpolationAdaptiveBins.h" +#ifdef DOUBLE_SPH +#include "single_photon_hit_double.h" +#endif + +#ifndef DOUBLE_SPH +#include "single_photon_hit.h" +#endif + +#include "etaInterpolationPosXY.h" +#include "noInterpolation.h" +//#include "etaInterpolationAdaptiveBins.h" +//#include "etaInterpolationRandomBins.h" using namespace std; #define NC 400 #define NR 400 +#define XTALK + int main(int argc, char *argv[]) { /** * trial.o [socket ip] [starting port number] [outfname] * */ - if (argc<7) { - cout << "Wrong usage! Should be: "<< argv[0] << " infile " << " etafile outfile runmin runmax ns" << endl; + if (argc<9) { + cout << "Wrong usage! Should be: "<< argv[0] << " infile " << " etafile outfile runmin runmax ns cmin cmax" << endl; return 1; } char infname[10000]; + char fname[10000]; char outfname[10000]; int runmin=atoi(argv[4]); int runmax=atoi(argv[5]); int nsubpix=atoi(argv[6]); + float cmin=atof(argv[7]); + float cmax=atof(argv[8]); int etabins=1000;//nsubpix*2*100; double etamin=-1, etamax=2; + double eta3min=-2, eta3max=2; int quad; double sum, totquad; double sDum[2][2]; double etax, etay, int_x, int_y; + double eta3x, eta3y, int3_x, int3_y, noint_x, noint_y; int ok; - + int f0=-1; int ix, iy, isx, isy; - - + int nframes=0, lastframe=-1; + double d_x, d_y, res=5, xx, yy; +#ifdef MANYFILES + int ff=1000; +#endif + int nph=0, badph=0, totph=0; FILE *f=NULL; +#ifdef DOUBLE_SPH + single_photon_hit_double cl(3,3); +#endif + +#ifndef DOUBLE_SPH single_photon_hit cl(3,3); - // etaInterpolationPosXY *interp=new etaInterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax); - etaInterpolationAdaptiveBins *interp=new etaInterpolationAdaptiveBins (NC, NR, nsubpix, etabins, etamin, etamax); +#endif + +#ifdef XTALK + int old_val[3][3]; + int new_val[3][3]; + double xcorr=0.04; + + // int ix=0; +#endif + + + eta2InterpolationPosXY *interp=new eta2InterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax); + eta3InterpolationPosXY *interp3=new eta3InterpolationPosXY(NC, NR, nsubpix, etabins, eta3min, eta3max); + noInterpolation *dummy=new noInterpolation(NC, NR, nsubpix); + noInterpolation *nointerp=new noInterpolation(NC, NR, nsubpix); + noInterpolation *mult=new noInterpolation(NC, NR, nsubpix); + //etaInterpolationAdaptiveBins *interp=new etaInterpolationAdaptiveBins (NC, NR, nsubpix, etabins, etamin, etamax); + //etaInterpolationRandomBins *interp=new etaInterpolationRandomBins (NC, NR, nsubpix, etabins, etamin, etamax); //#ifndef FF cout << "read ff " << argv[2] << endl; - interp->readFlatField(argv[2]); + sprintf(fname,"%s_eta2.tiff",argv[2]); + interp->readFlatField(fname); interp->prepareInterpolation(ok); + + sprintf(fname,"%s_eta3.tiff",argv[2]); + interp3->readFlatField(fname); + interp3->prepareInterpolation(ok); //#endif int *img; @@ -70,20 +118,99 @@ int main(int argc, char *argv[]) { #endif - - for (int irun=runmin; iruncalcQuad(cl.get_cluster(), sum, totquad, sDum); - if (sum>200 && sum<580) { - interp->getInterpolatedPosition(cl.x,cl.y, totquad,quad,cl.get_cluster(),int_x, int_y); - interp->addToImage(int_x, int_y); + totph++; + if (lastframe!=cl.iframe) { + lastframe=cl.iframe; + // cout << cl.iframe << endl; + // f0=cl.iframe; + if (nframes==0) f0=lastframe; + nframes++; +#ifdef MANYFILES + if (nframes%ff==0) { + sprintf(outfname,argv[3],irun,nframes-ff); + cout << outfname << endl; + interp->writeInterpolatedImage(outfname); + interp->clearInterpolatedImage(); } +#endif + + } +#ifdef XTALK + if ((cl.x+1)%25!=0) { + for (int ix=-1; ix<2; ix++) { + for (int iy=-1; iy<2; iy++) { + old_val[iy+1][ix+1]=cl.get_data(ix,iy); + if (ix>=0) { + new_val[iy+1][ix+1]=old_val[iy+1][ix+1]-old_val[iy+1][ix]*xcorr; + cl.set_data(new_val[iy+1][ix+1],ix,iy); + } + } + } + } +#endif + quad=interp->calcQuad(cl.get_cluster(), sum, totquad, sDum); + if (sum>cmin && totquad/sum>0.8 && totquad/sum<1.2 && totquad200 && sum<580) { + // interp->getInterpolatedPosition(cl.x,cl.y, totquad,quad,cl.get_cluster(),int_x, int_y); + interp->getInterpolatedPosition(cl.x,cl.y, cl.get_cluster(),int_x, int_y); + interp->addToImage(int_x, int_y); + interp3->getInterpolatedPosition(cl.x,cl.y, cl.get_cluster(),int3_x, int3_y); + interp3->addToImage(int3_x, int3_y); + nointerp->getInterpolatedPosition(cl.x,cl.y, cl.get_cluster(),noint_x, noint_y); + nointerp->addToImage(noint_x, noint_y); + + + d_x= (int_x-int3_x)*25.; + d_y= (int_y-int3_y)*25.; + dummy->calcEta(totquad, sDum, etax, etay); + xx=int_x; + yy=int_y; + if (etax<0.1 || etax>0.9) xx=int3_x; + if (etay<0.1 || etay>0.9) yy=int3_y; + dummy->addToImage(xx,yy); + + if (d_x>res || d_x<-res || d_y>res || d_y<-res) { + badph++; + // cout << "delta (um): "<< d_x << " " << d_y << " " << cl.x << " " << cl.y << endl; + // cout << sum << " " << totquad << " " << etax << " "<< etay << endl; + // //cout<< int_x << " " << int_y << " " << int3_x << " " << int3_y << endl; + } + mult->addToImage(noint_x, noint_y); + + if (nph%1000000==0) cout << nph << endl; + if (nph%100000000==0) { + sprintf(outfname,"%s_inteta2.tiff", argv[3]); + interp->writeInterpolatedImage(outfname); + sprintf(outfname,"%s_inteta3.tiff", argv[3]); + interp3->writeInterpolatedImage(outfname); + sprintf(outfname,"%s_mix.tiff", argv[3]); + dummy->writeInterpolatedImage(outfname); + sprintf(outfname,"%s_noint.tiff", argv[3]); + nointerp->writeInterpolatedImage(outfname); + sprintf(outfname,"%s_mult.tiff", argv[3]); + mult->writeInterpolatedImage(outfname); + } + + } else { + mult->getInterpolatedPosition(cl.x,cl.y, cl.get_cluster(),int_x, int_y); + for (int imult=0; imult<2.*sum/(cmax+cmin); imult++) mult->addToImage(int_x, int_y); + } + + } fclose(f); #ifdef FF @@ -147,8 +274,22 @@ int main(int argc, char *argv[]) { } #endif #ifndef FF + +#ifdef MANYFILES + sprintf(outfname,argv[3],irun,nframes-ff); + cout << outfname << endl; +#endif + sprintf(outfname,"%s_inteta2.tiff", argv[3]); interp->writeInterpolatedImage(outfname); - img=interp->getInterpolatedImage(); + sprintf(outfname,"%s_inteta3.tiff", argv[3]); + interp3->writeInterpolatedImage(outfname); + sprintf(outfname,"%s_mix.tiff", argv[3]); + dummy->writeInterpolatedImage(outfname); + sprintf(outfname,"%s_mult.tiff", argv[3]); + mult->writeInterpolatedImage(outfname); + +#ifndef MANYFILES + img=interp->getInterpolatedImage(); for (ix=0; ixclearInterpolatedImage(); + interp3->clearInterpolatedImage(); + dummy->clearInterpolatedImage(); + mult->clearInterpolatedImage(); - } + } else + cout << "could not open file " << infname << endl; } - - + cout << irun << " " << runmax << endl; +#ifndef MANYFILES sprintf(outfname,argv[3],11111); WriteToTiff(totimg, outfname,NC*nsubpix,NR*nsubpix); - +#endif + cout << "Filled " << nph << " (/"<< totph <<") of which " << badph << " badly interpolated " << endl; return 0; } diff --git a/slsDetectorCalibration/moenchExecutables/moench03MakeEta.cpp b/slsDetectorCalibration/moenchExecutables/moench03MakeEta.cpp index 4f1119471..c2f2f92cd 100644 --- a/slsDetectorCalibration/moenchExecutables/moench03MakeEta.cpp +++ b/slsDetectorCalibration/moenchExecutables/moench03MakeEta.cpp @@ -11,6 +11,7 @@ using namespace std; #define NC 400 #define NR 400 +#define XTALK int main(int argc, char *argv[]) { /** @@ -20,6 +21,7 @@ int main(int argc, char *argv[]) { int nsubpix=10; int etabins=nsubpix*100; double etamin=-1, etamax=2; + double eta3min=-2, eta3max=2; int quad; double sum, totquad; double sDum[2][2]; @@ -27,16 +29,29 @@ int main(int argc, char *argv[]) { double etax, etay; int runmin, runmax; single_photon_hit cl(3,3); - - if (argc<5) { - cout << "Wrong usage! Should be: "<< argv[0] << " infile " << " outfile runmin runmax" << endl; +int iph=0; + + if (argc<7) { + cout << "Wrong usage! Should be: "<< argv[0] << " infile " << " outfile runmin runmax cmin cmax" << endl; return 1; } - etaInterpolationPosXY *interp=new etaInterpolationPosXY(NR, NC, nsubpix, etabins, etamin, etamax); + eta2InterpolationPosXY *interp2=new eta2InterpolationPosXY(NR, NC, nsubpix, etabins, etamin, etamax); + cout << "###########"<< endl; + eta3InterpolationPosXY *interp3=new eta3InterpolationPosXY(NR, NC, nsubpix, etabins, eta3min, eta3max); + // cout << eta3min << " " << eta3max << endl; runmin=atoi(argv[3]); runmax=atoi(argv[4]); + double cmin=atof(argv[5]); //200 + double cmax=atof(argv[6]); //3000 +#ifdef XTALK + int old_val[3][3]; + int new_val[3][3]; + double xcorr=0.04; + + // int ix=0; +#endif FILE *f; for (int i=runmin; icalcQuad(cl.get_cluster(), sum, totquad, sDum); - if (sum>200 && sum<580 && cl.y<350) - interp->addToFlatField(cl.get_cluster(),etax, etay); +#ifdef XTALK + if ((cl.x+1)%25!=0) { + for (int ix=-1; ix<2; ix++) { + for (int iy=-1; iy<2; iy++) { + old_val[iy+1][ix+1]=cl.get_data(ix,iy); + if (ix>=0) { + new_val[iy+1][ix+1]=old_val[iy+1][ix+1]-old_val[iy+1][ix]*xcorr; + cl.set_data(new_val[iy+1][ix+1],ix,iy); + } + } + } + } +#endif + quad=interp2->calcQuad(cl.get_cluster(), sum, totquad, sDum); + + if (sum>cmin && totquad/sum>0.8 && totquad/sum<1.2 && totquadaddToFlatField(cl.get_cluster(),etax, etay); + // if (etax>0.49 && etax<0.51 && etay>0.49 && etay<0.51 ) { + // cout << cl.y << " " << cl.x << " " << quad << " "<< totquad << " " <addToFlatField(cl.get_cluster(),etax, etay); + iph++; + if (iph%1000000==0) cout << iph << endl; + if (iph%100000000==0) { + sprintf(fname,"%s_eta2.tiff",argv[2]); + interp2->writeFlatField(fname); + sprintf(fname,"%s_eta3.tiff",argv[2]); + interp3->writeFlatField(fname); + } + // if (iph>1E8) break; + } + // } + } fclose(f); - interp->writeFlatField(argv[2]); } else cout << "could not open file " << fname << endl; } - interp->writeFlatField(argv[2]); + sprintf(fname,"%s_eta2.tiff",argv[2]); + interp2->writeFlatField(fname); + sprintf(fname,"%s_eta3.tiff",argv[2]); + interp3->writeFlatField(fname); return 0; } diff --git a/slsDetectorCalibration/moenchExecutables/moench03OnTheFlyAnalysis.C b/slsDetectorCalibration/moenchExecutables/moench03OnTheFlyAnalysis.C deleted file mode 100644 index b247a9892..000000000 --- a/slsDetectorCalibration/moenchExecutables/moench03OnTheFlyAnalysis.C +++ /dev/null @@ -1,160 +0,0 @@ -#include -#include -#include -#include -#include -#include -//#include -//#include -//#include -#include - -#include "moench03Ctb10GbT1Data.h" - -#include "interpolatingDetector.h" -#include "etaInterpolationPosXY.h" -#include "linearInterpolation.h" -#include "noInterpolation.h" -#include "multiThreadedDetector.h" - -#include - -#define NC 400 -#define NR 400 - -#include "tiffIO.h" - - -void *moenchProcessFrame() { - char fname[10000]; - strcpy(fname,"/mnt/moench_data/m03-15_mufocustube/plant_40kV_10uA/m03-15_100V_g4hg_300us_dtf_0.raw"); - - int nthreads=3; - - int nph, nph1; - single_photon_hit clusters[NR*NC]; - // cout << "hits "<< endl; - int etabins=550; - double etamin=-1, etamax=2; - int nsubpix=4; - float *etah=new float[etabins*etabins]; - // cout << "etah "<< endl; - cout << "image size "<< nsubpix*nsubpix*NC*NR << endl; - float *image=new float[nsubpix*nsubpix*NC*NR]; - int *heta, *himage; - - moench03Ctb10GbT1Data *decoder=new moench03Ctb10GbT1Data(); - // cout << "decoder "<< endl; - etaInterpolationPosXY *interp=new etaInterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax); - // cout << "interp "<< endl; - //linearInterpolation *interp=new linearInterpolation(NC, NR, nsubpix); - //noInterpolation *interp=new noInterpolation(NC, NR, nsubpix); - interp->readFlatField("/scratch/eta_100.tiff",etamin,etamax); - interpolatingDetector *filter=new interpolatingDetector(decoder,interp, 5, 1, 0, 1000, 10); - filter->readPedestals("/scratch/ped_100.tiff"); - cout << "filter "<< endl; - - - - char *buff; - int nf=0; - int ok=0; - ifstream filebin; - std::time_t end_time; - - int iFrame=-1; - int np=-1; - - - filter->newDataSet(); - - - multiThreadedDetector *mt=new multiThreadedDetector(filter,nthreads,100); - nph=0; - nph1=0; - //int np; - int iph; - - cout << "file name " << fname << endl; - filebin.open((const char *)(fname), ios::in | ios::binary); - if (filebin.is_open()) - cout << "Opened file " << fname<< endl; - else - cout << "Could not open file " << fname<< endl; - mt->setFrameMode(eFrame); - mt->prepareInterpolation(ok); - mt->StartThreads(); - mt->popFree(buff); - - while ((decoder->readNextFrame(filebin, iFrame, np, buff)) && nf<1.5E4) { - if (nf<9E3) - ; - else { - - // if (nf>1.1E4 && ok==0) { - // mt->prepareInterpolation(ok); - // mt->setFrameMode(eFrame); - // //ok=1; - // } - - mt->pushData(buff); - mt->nextThread(); - // cout << " " << (void*)buff; - mt->popFree(buff); - - // if (ok==0) { - // cout << "**************************************************************************"<< endl; - // heta=interp->getFlatField(); - // // for (int ii=0; iiprepareInterpolation(ok); - // cout << "**************************************************************************"<< endl; - // std::time(&end_time); - // cout << std::ctime(&end_time) << " " << nf << endl; - // } - // filter->processData(buff,eFrame); - // } - - // nph+=nph1; - } - if (nf%1000==0) { - std::time(&end_time); - cout << std::ctime(&end_time) << " " << nf << endl; - } - - nf++; - //delete [] buff; - iFrame=-1; - } - - if (filebin.is_open()) - filebin.close(); - else - cout << "could not open file " << fname << endl; - - mt->StopThreads(); - - char tit[10000]; - sprintf(tit,"/scratch/int_image_mt%d.tiff",nthreads); - - mt->writeInterpolatedImage(tit); - // delete [] etah; - - // delete interp; - //delete decoder; - //cout << "Read " << nf << " frames" << endl; - return NULL; -} - -int main(int argc, char *argv[]){ - - moenchProcessFrame(); - -} diff --git a/slsDetectorCalibration/moenchExecutables/readClusters.C b/slsDetectorCalibration/moenchExecutables/readClusters.C deleted file mode 100644 index f28cb95e0..000000000 --- a/slsDetectorCalibration/moenchExecutables/readClusters.C +++ /dev/null @@ -1,267 +0,0 @@ -#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"); - int iph=0; - int ns=4; - double px, py; - double left, right, top, bottom; - if (f) { - int x1,y1; - if (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,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,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); - // hsp2->SetLineColor(3); - // hsp3->SetLineColor(4); - TCanvas *c=new TCanvas(); - h2->Draw("colz"); - TCanvas *c1=new TCanvas(); - hsp->Draw(); - c1->SetLogy(kTRUE); - // hsp1->Draw("same"); - // hsp2->Draw("same"); - // hsp3->Draw("same"); - TCanvas *c2=new TCanvas(); - hf->Draw(); - // hint->Draw("colz"); - //c2->SetLogz(kTRUE); - single_photon_hit cl(3,3); - double tot; - int w; - 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; - //if (cl.x>80 && cl.x<280 && cl.y>80 && cl.y<300) { - tot=0; /* - left=0; - right=0; - top=0; - bottom=0;*/ - tl=0; tr=0; bl=0; br=0; - for (int ix=-1; ix<2; ix++) - for (int iy=-1; iy<2; iy++){ - tot+=cl.get_data(ix,iy); - if (ix<=0 && iy<=0) bl+=cl.get_data(ix,iy); - if (ix<=0 && iy>=0) tl+=cl.get_data(ix,iy); - if (ix>=0 && iy<=0) br+=cl.get_data(ix,iy); - if (ix>=0 && iy>=0) tr+=cl.get_data(ix,iy); - - //if (ix<0) left+=cl.get_data(ix,iy); - // if (ix>0) right+=cl.get_data(ix,iy); - // if (iy<0) bottom+=cl.get_data(ix,iy); - // if (iy>0) top+=cl.get_data(ix,iy); - - } - qt=bl; - if (br>qt) qt=br; - if (tl>qt) qt=tl; - if (tr>qt) qt=tr; - /* - px=(-left+right)/tot; - py=(-bottom+top)/tot;*/ - //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++; - //(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 - hsp->Fill(qt); - - // hsp2->Fill(cl.get_data(0,0)); - // } else { - // hsp1->Fill(qt); - // hsp3->Fill(cl.get_data(0,0)); - // } - //if (cl.x>160 && cl.x<260 && cl.y>30 && cl.y<80 && tot>0) - //if (w==1) { - // if (w==2) - h2->Fill(cl.x,cl.y,w); - // } - //} - // hint->Fill(px+cl.x,py+cl.y,w); - // if (cl.y<350) - // hff->Fill(px,py,w); - - // } - - //h2mult->Fill(cl.x,cl.y,w); - - } - iph+=w; - if (iph%100000==0) { - c->Modified(); - c->Update(); - c1->Modified();; - c1->Update(); - c2->Modified();; - c2->Update(); - } - // if (iph>1E7) - // break; - //} - // if (iph>0.5E7) break; - } - } - fclose(f); - // hff->Scale(hff->GetNbinsX()*hff->GetNbinsY()/hff->Integral()); - // TH2F *hint2=(TH2F*)hint->Clone("hint2"); - // double ff; - // for (int ibx=0; ibxGetNbinsX(); ibx++) { - // for (int iby=0; ibyGetNbinsY(); iby++) { - // ff=hff->GetBinContent(ibx%ns+1, iby%ns+1); - // // cout << ibx << " " << iby << " " << ibx%ns << " " << iby%ns << " " << ff << endl; - // if (ff>0) - // hint2->SetBinContent(ibx+1, iby+1,hint->GetBinContent(ibx+1,iby+1)/ff); - // } - // } - - - - - - - return h2; - -} else - cout << "could not open file " << fname << endl; - return NULL; - - -} - -// TH2F *getEta(char *fname, int nx, int ny, TH2F *h2=NULL) { -// int ns=10; -// slsInterpolation *inte=new etaInterpolationPosXY(nx,ny,ns, - -// FILE *f=fopen(fname,"r"); -// int iph=0; -// int ns=4; -// double px, py; -// double left, right, top, bottom; -// if (f) { -// int x1,y1; -// if (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); -// // 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); -// TCanvas *c=new TCanvas(); -// c->SetLogz(kTRUE); -// h2->Draw("colz"); -// TCanvas *c1=new TCanvas(); -// hsp->Draw(); -// TCanvas *c2=new TCanvas(); -// hint->Draw("colz"); -// c2->SetLogz(kTRUE); -// single_photon_hit cl(3,3); -// double tot; -// int w; -// double phw=340, phs=62; -// while (cl.read(f)) { -// //cl.get_pixel(x1, y1); -// //cout << cl.iframe << " " << cl.x << " " << cl.y << endl; -// //if (cl.x>80 && cl.x<280 && cl.y>80 && cl.y<300) { -// tot=0; -// left=0; -// right=0; -// top=0; -// bottom=0; -// for (int ix=-1; ix<2; ix++) -// for (int iy=-1; iy<2; iy++){ -// tot+=cl.get_data(ix,iy); -// if (ix<0) left+=cl.get_data(ix,iy); -// if (ix>0) right+=cl.get_data(ix,iy); -// if (iy<0) bottom+=cl.get_data(ix,iy); -// if (iy>0) top+=cl.get_data(ix,iy); - -// } -// px=(-left+right)/tot; -// py=(-bottom+top)/tot; -// //max at 340 -// if (tot>200) { -// w=(tot+3.5*phs)/phw; -// } else -// w=0; -// if (w) { -// hsp->Fill(tot); -// if (w==1) { -// h2->Fill(cl.x,cl.y,w); -// hint->Fill(px+cl.x,py+cl.y,w); -// } - -// h2mult->Fill(cl.x,cl.y,w); - -// // if (cl.y<350) -// // hff->Fill(px,py,w); -// //} -// iph+=w; -// if (iph%100000==0) { -// // c->Modified(); -// // c->Update(); -// c1->Modified();; -// c1->Update(); -// // c2->Modified();; -// // c2->Update(); -// } -// } -// // if (iph>0.5E7) break; -// } -// fclose(f); -// // hff->Scale(hff->GetNbinsX()*hff->GetNbinsY()/hff->Integral()); -// // TH2F *hint2=(TH2F*)hint->Clone("hint2"); -// // double ff; -// // for (int ibx=0; ibxGetNbinsX(); ibx++) { -// // for (int iby=0; ibyGetNbinsY(); iby++) { -// // ff=hff->GetBinContent(ibx%ns+1, iby%ns+1); -// // // cout << ibx << " " << iby << " " << ibx%ns << " " << iby%ns << " " << ff << endl; -// // if (ff>0) -// // hint2->SetBinContent(ibx+1, iby+1,hint->GetBinContent(ibx+1,iby+1)/ff); -// // } -// // } - - - - - - -// return h2; - -// } else -// cout << "could not open file " << fname << endl; -// return NULL; - - -// } diff --git a/slsDetectorCalibration/pedestalSubtraction.h b/slsDetectorCalibration/pedestalSubtraction.h index d05f9cfca..062b0d84e 100644 --- a/slsDetectorCalibration/pedestalSubtraction.h +++ b/slsDetectorCalibration/pedestalSubtraction.h @@ -51,8 +51,11 @@ class pedestalSubtraction { /** sets the moving average */ virtual void setPedestalRMS(double rms) {stat.SetRMS(rms);} - - + /**sets/gets the number of samples for the moving average + \returns actual number of samples for the moving average + */ + virtual int getNumpedestals() {return stat.NumDataValues();}; + private: MovingStat stat; /**< approximated moving average struct */ diff --git a/slsDetectorCalibration/singlePhotonDetector.h b/slsDetectorCalibration/singlePhotonDetector.h index 7bc09ebd4..f634ebdbd 100644 --- a/slsDetectorCalibration/singlePhotonDetector.h +++ b/slsDetectorCalibration/singlePhotonDetector.h @@ -58,7 +58,7 @@ public analogDetector { int sign=1, commonModeSubtraction *cm=NULL, int nped=1000, - int nd=100, int nnx=-1, int nny=-1, double *gm=NULL) : analogDetector(d, sign, cm, nped, nnx, nny, gm), nDark(nd), eventMask(NULL),nSigma (nsigma), clusterSize(csize), clusterSizeY(csize), cluster(NULL), quad(UNDEFINED_QUADRANT), tot(0), quadTot(0) { + int nd=100, int nnx=-1, int nny=-1, double *gm=NULL) : analogDetector(d, sign, cm, nped, nnx, nny, gm), nDark(nd), eventMask(NULL),nSigma (nsigma), clusterSize(csize), clusterSizeY(csize), clusters(NULL), quad(UNDEFINED_QUADRANT), tot(0), quadTot(0) { @@ -74,7 +74,7 @@ public analogDetector { // cluster=new single_photon_hit(clusterSize,clusterSizeY); clusters=new single_photon_hit[nx*ny]; - cluster=clusters; + // cluster=clusters; setClusterSize(csize); nphTot=0; nphFrame=0; @@ -82,7 +82,7 @@ public analogDetector { /** destructor. Deletes the cluster structure, the pdestalSubtraction and the image array */ - virtual ~singlePhotonDetector() {delete cluster; for (int i=0; i { // cluster=new single_photon_hit(clusterSize,clusterSizeY); clusters=new single_photon_hit[nx*ny]; - cluster=clusters; + // cluster=clusters; setClusterSize(clusterSize); @@ -147,8 +147,8 @@ public analogDetector { if (n%2==0) n+=1; clusterSize=n; - if (cluster) - delete cluster; + // if (clusters) + // delete [] clusters; if (ny>clusterSize) clusterSizeY=clusterSize; else @@ -192,29 +192,37 @@ public analogDetector { double max=0, tl=0, tr=0, bl=0,br=0, v; - if (thr>=0) { + int cm=0; + if (cmSub) cm=1; + + if (thr>0) { cy=1; cs=1; ccs=1; ccy=1; } - if (iframe0) { - + newFrame(); + if (cmSub) { + addToCommonMode(data); + } for (int ix=xmin; ix::getNPhotons(data,ix,iy); + nn=val/tthr;//analogDetector::getNPhotons(data,ix,iy); nph[ix+nx*iy]+=nn; rest[iy][ix]=(val-nn*tthr); + nphFrame+=nn; + nphTot+=nn; } } // } @@ -235,10 +243,10 @@ public analogDetector { 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); + //clusters->set_data(rest[iy+ir][ix+ic], ic, ir); - v=rest[iy+ir][ix+ic];//cluster->get_data(ic,ir); + v=rest[iy+ir][ix+ic];//clusters->get_data(ic,ir); tot+=v; if (ir<=0 && ic<=0) @@ -278,14 +286,15 @@ public analogDetector { if (max>tthr || tot>sqrt(ccy*ccs)*tthr || quadTot>sqrt(cy*cs)*tthr) { eventMask[iy][ix]=PHOTON; nph[ix+nx*iy]++; - nphFrame+=nph[ix+nx*iy]; - nphTot+=nph[ix+nx*iy]; + nphFrame++; + nphTot++; } } } } - cout << iframe << " " << nph << endl; + // cout << iframe << " " << nphFrame << " " << nphTot << endl; + //cout << iframe << " " << nph << endl; } else return getClusters(data, nph); } return NULL; @@ -329,10 +338,10 @@ public analogDetector { eventMask[iy][ix]=PEDESTAL; - cluster->x=ix; - cluster->y=iy; - cluster->rms=getPedestalRMS(ix,iy); - cluster->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++) { @@ -341,8 +350,8 @@ public analogDetector { v=subtractPedestal(data, ix+ic, iy+ir); - cluster->set_data(v, ic, ir); - // v=cluster->get_data(ic,ir); + clusters->set_data(v, ic, ir); + // v=clusters->get_data(ic,ir); tot+=v; if (ir<=0 && ic<=0) bl+=v; @@ -357,7 +366,7 @@ public analogDetector { max=v; } if (ir==0 && ic==0) { - if (v<-nSigma*cluster->rms) + if (v<-nSigma*clusters->rms) eventMask[iy][ix]=NEGATIVE_PEDESTAL; } } @@ -378,8 +387,8 @@ public analogDetector { quadTot=tr; } - if (max>nSigma*cluster->rms || tot>sqrt(clusterSizeY*clusterSize)*nSigma*cluster->rms || quadTot>cy*cs*nSigma*cluster->rms) { - if (cluster->get_data(0,0)>=max) { + 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; @@ -447,7 +456,7 @@ int *getClusters(char *data, int *ph=NULL) { (clusters+nph)->rms=getPedestalRMS(ix,iy); - cluster=clusters+nph; + // cluster=clusters+nph; for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) { @@ -473,7 +482,7 @@ int *getClusters(char *data, int *ph=NULL) { if (ir==0 && ic==0) { - if (*v<-nSigma*cluster->rms) + if (*v<-nSigma*(clusters+nph)->rms) eventMask[iy][ix]=NEGATIVE_PEDESTAL; } @@ -494,7 +503,7 @@ int *getClusters(char *data, int *ph=NULL) { (clusters+nph)->quadTot=tr; } - if (max>nSigma*cluster->rms || tot>sqrt(clusterSizeY*clusterSize)*nSigma*cluster->rms || ((clusters+nph)->quadTot)>sqrt(cy*cs)*nSigma*cluster->rms) { + 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; @@ -559,7 +568,7 @@ int *getClusters(char *data, int *ph=NULL) { \param ir y coordinate (center is (0,0)) \returns cluster element */ - double getClusterElement(int ic, int ir=0){return cluster->get_data(ic,ir);}; + double getClusterElement(int ic, int ir=0){return clusters->get_data(ic,ir);}; /** returns event mask for the given pixel \param ic x coordinate (center is (0,0)) @@ -581,25 +590,25 @@ int *getClusters(char *data, int *ph=NULL) { if (iFrame) tall->Branch("iFrame",iFrame,"iframe/I"); else - tall->Branch("iFrame",&(cluster->iframe),"iframe/I"); + tall->Branch("iFrame",&(clusters->iframe),"iframe/I"); - tall->Branch("x",&(cluster->x),"x/I"); - tall->Branch("y",&(cluster->y),"y/I"); + tall->Branch("x",&(clusters->x),"x/I"); + tall->Branch("y",&(clusters->y),"y/I"); char tit[100]; sprintf(tit,"data[%d]/D",clusterSize*clusterSizeY); - tall->Branch("data",cluster->data,tit); - tall->Branch("pedestal",&(cluster->ped),"pedestal/D"); - tall->Branch("rms",&(cluster->rms),"rms/D"); - tall->Branch("tot",&(cluster->tot),"tot/D"); - tall->Branch("quadTot",&(cluster->quadTot),"quadTot/D"); - tall->Branch("quad",&(cluster->quad),"quad/I"); + tall->Branch("data",clusters->data,tit); + tall->Branch("pedestal",&(clusters->ped),"pedestal/D"); + tall->Branch("rms",&(clusters->rms),"rms/D"); + tall->Branch("tot",&(clusters->tot),"tot/D"); + tall->Branch("quadTot",&(clusters->quadTot),"quadTot/D"); + tall->Branch("quad",&(clusters->quad),"quad/I"); return tall; }; #else /** write cluster to filer \param f file pointer */ - void writeCluster(FILE* f){cluster->write(f);}; + void writeCluster(FILE* f){clusters->write(f);}; /** write clusters to file @@ -643,7 +652,7 @@ void writeClusters(FILE *f){for (int i=0; iwrite(f double nSigma; /**< number of sigma parameter for photon discrimination */ int clusterSize; /**< cluster size in the x direction */ int clusterSizeY; /**< cluster size in the y direction i.e. 1 for strips, clusterSize for pixels */ - single_photon_hit *cluster; /**< single photon hit data structure */ + // single_photon_hit *cluster; /**< single photon hit data structure */ single_photon_hit *clusters; /**< single photon hit data structure */ quadrant quad; /**< quadrant where the photon is located */ double tot; /**< sum of the 3x3 cluster */ diff --git a/slsDetectorCalibration/single_photon_hit.h b/slsDetectorCalibration/single_photon_hit.h index c992123d8..696868b87 100644 --- a/slsDetectorCalibration/single_photon_hit.h +++ b/slsDetectorCalibration/single_photon_hit.h @@ -1,5 +1,7 @@ #ifndef SINGLE_PHOTON_HIT_H #define SINGLE_PHOTON_HIT_H +#include +#include typedef double double32_t; typedef float float32_t; @@ -26,7 +28,7 @@ class single_photon_hit { \param nx cluster size in x direction \param ny cluster size in y direction (defaults to 1 for 1D detectors) */ - single_photon_hit(int nx=3, int ny=3): dx(nx), dy(ny) {data=new double[dx*dy];}; + single_photon_hit(int nx=3, int ny=3): dx(nx), dy(ny) {data=new int[dx*dy];}; ~single_photon_hit(){delete [] data;}; /**< destructor, deletes the data array */ @@ -35,8 +37,9 @@ class single_photon_hit { */ size_t write(FILE *myFile) { //fwrite((void*)this, 1, 3*sizeof(int)+4*sizeof(double)+sizeof(quad), myFile); - if (fwrite((void*)this, 1, 3*sizeof(int), myFile)) - return fwrite((void*)data, 1, dx*dy*sizeof(double), myFile); + + if (fwrite((void*)this, 1, sizeof(int)+2*sizeof(int16_t), myFile)) + return fwrite((void*)data, 1, dx*dy*sizeof(int), myFile); return 0; }; @@ -47,8 +50,8 @@ class single_photon_hit { size_t read(FILE *myFile) { //fread((void*)this, 1, 3*sizeof(int)+4*sizeof(double)+sizeof(quad), myFile); - if (fread((void*)this, 1, 3*sizeof(int), myFile)) - return fread((void*)data, 1, dx*dy*sizeof(double), myFile); + if (fread((void*)this, 1, sizeof(int)+2*sizeof(int16_t), myFile)) + return fread((void*)data, 1, dx*dy*sizeof(int), myFile); return 0; }; @@ -60,7 +63,7 @@ class single_photon_hit { */ void set_data(double v, int ix, int iy=0){data[(iy+dy/2)*dx+ix+dx/2]=v;}; - void set_cluster_size(int nx, int ny) {if (nx>0) dx=nx; if (ny>0) dy=ny; delete [] data; data=new double[dx*dy];}; + void set_cluster_size(int nx, int ny) {if (nx>0) dx=nx; if (ny>0) dy=ny; delete [] data; data=new int[dx*dy];}; void get_cluster_size(int &nx, int &ny) {nx=dx; ny=dy;}; void get_pixel(int &x1, int &y1) {x1=x; y1=y;}; @@ -71,11 +74,11 @@ class single_photon_hit { \returns value of the cluster element */ double get_data(int ix, int iy=0){return data[(iy+dy/2)*dx+ix+dx/2];}; - double *get_cluster() {return data;}; + int *get_cluster() {return data;}; int iframe; /**< frame number */ - int x; /**< x-coordinate of the center of hit */ - int y; /**< x-coordinate of the center of hit */ + int16_t x; /**< x-coordinate of the center of hit */ + int16_t y; /**< x-coordinate of the center of hit */ double rms; /**< noise of central pixel l -- at some point it can be removed*/ double ped; /**< pedestal of the central pixel -- at some point it can be removed*/ double tot; /**< sum of the 3x3 cluster */ @@ -83,7 +86,7 @@ class single_photon_hit { double quadTot; /**< sum of the maximum 2x2cluster */ int dx; /**< size of data cluster in x */ int dy; /**< size of data cluster in y */ - double *data; /**< pointer to data */ + int *data; /**< pointer to data */ }; diff --git a/slsDetectorCalibration/tiffIO.cpp b/slsDetectorCalibration/tiffIO.cpp index b0dd03b71..6521df159 100644 --- a/slsDetectorCalibration/tiffIO.cpp +++ b/slsDetectorCalibration/tiffIO.cpp @@ -1,5 +1,6 @@ - +#ifndef MY_TIFF_IO_H #include "tiffIO.h" +#endif #include using namespace std; // #undef cbf_failnez diff --git a/slsDetectorSoftware/mythenDetectorServer/Makefile b/slsDetectorSoftware/mythenDetectorServer/Makefile index da0c0f291..b72e33681 100755 --- a/slsDetectorSoftware/mythenDetectorServer/Makefile +++ b/slsDetectorSoftware/mythenDetectorServer/Makefile @@ -13,7 +13,7 @@ INSTMODE= 0777 SRCS= server.c server_funcs.c communication_funcs.c firmware_funcs.c mcb_funcs.c trimming_funcs.c sharedmemory.c OBJS= $(SRCS:%.c=%.o) -VFLAGS= +VFLAGS= #-DVERBOSE #-DVERYVERBOSE CFLAGS+= -Wall -DC_ONLY -DMCB_FUNCS -DDACS_INT $(VFLAGS) @@ -32,7 +32,7 @@ versioning: $(PROGS): $(OBJS) echo $(OBJS) - $(CC) $(LDFLAGS) $^ $(LDLIBS) $(CFLAGS) -o $@ + $(CC) $(LDFLAGS) $^ $(LDLIBS) $(CFLAGS) -o $@ -DVERBOSE install: $(PROGS) $(INSTALL) -d $(INSTDIR) diff --git a/slsDetectorSoftware/mythenDetectorServer/firmware_funcs.c b/slsDetectorSoftware/mythenDetectorServer/firmware_funcs.c index 49cbdb982..6e90105cb 100755 --- a/slsDetectorSoftware/mythenDetectorServer/firmware_funcs.c +++ b/slsDetectorSoftware/mythenDetectorServer/firmware_funcs.c @@ -17,6 +17,7 @@ #include #include +//#define VERBOSE //for memory mapping u_int64_t CSP0BASE; @@ -702,9 +703,9 @@ int setNMod(int n) { shiftfifo=SHIFTFIFO; -#ifdef VERBOSE + //#ifdef VERBOSE printf("SetNMod called arg %d -- dr %d shiftfifo %d\n",n,dynamicRange,shiftfifo); -#endif + //#endif if (n>=0 && n<=ntot) { nModX=n; @@ -713,12 +714,16 @@ int setNMod(int n) { reg=bus_r(FIFO_COUNTR_REG_OFF+(ififo<>FIFO_NM_OFF, (reg&FIFO_NC_MASK)>>FIFO_NC_OFF, FIFO_COUNTR_REG_OFF+(ififo<>FIFO_NM_OFF, (reg&FIFO_NC_MASK)>>FIFO_NC_OFF, FIFO_COUNTR_REG_OFF+(ififo<