detectors/slsDetectorPackage
25
0
Fork 0
mirror of https://github.com/slsdetectorgroup/slsDetectorPackage.git synced 2025-04-20 02:40:03 +02:00
Code Issues Actions Packages Releases Activity

moench zmq process implemented - command implementation is missing

Browse Source
This commit is contained in:
bergamaschi 2018-09-13 16:07:43 +02:00
parent f288390255
commit 19e7ced332
12 changed files with 867 additions and 738 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
slsDetectorCalibration/analogDetector.h
View File

@ -80,7 +80,6 @@ template <class dataType> class analogDetector {
fMode=ePedestal; fMode=ePedestal;
thr=0; thr=0;
myFile=NULL; myFile=NULL;
fm=new pthread_mutex_t ;
#ifdef ROOTSPECTRUM #ifdef ROOTSPECTRUM
hs=new TH2F("hs","hs",2000,-100,10000,nx*ny,-0.5,nx*ny-0.5); hs=new TH2F("hs","hs",2000,-100,10000,nx*ny,-0.5,nx*ny-0.5);
#ifdef ROOTCLUST #ifdef ROOTCLUST
@ -128,7 +127,6 @@ template <class dataType> class analogDetector {
// nSigma=orig->nSigma; // nSigma=orig->nSigma;
fMode=orig->fMode; fMode=orig->fMode;
myFile=orig->myFile; myFile=orig->myFile;
fm=orig->fm;
stat=new pedestalSubtraction*[ny]; stat=new pedestalSubtraction*[ny];
@ -329,7 +327,8 @@ template <class dataType> class analogDetector {
for (int ix=xmin; ix<xmax; ix++) { for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) { for (int iy=ymin; iy<ymax; iy++) {
// if (getNumpedestals(ix,iy)>0) // if (getNumpedestals(ix,iy)>0)
addToCommonMode(data, ix, iy); if (det->isGood(ix,iy))
addToCommonMode(data, ix, iy);
} }
} }
//cout << "cm " << getCommonMode(0,0) << " " << getCommonMode(1,0) << endl; //cout << "cm " << getCommonMode(0,0) << " " << getCommonMode(1,0) << endl;
@ -701,12 +700,14 @@ template <class dataType> class analogDetector {
// cout << ymin << " " << ymax << endl; // cout << ymin << " " << ymax << endl;
for (int ix=xmin; ix<xmax; ix++) { for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) { for (int iy=ymin; iy<ymax; iy++) {
addToPedestal(data,ix,iy,1); if (det->isGood(ix,iy)) {
addToPedestal(data,ix,iy,1);
//if (ix==10 && iy==10) //if (ix==10 && iy==10)
// cout <<ix << " " << iy << " " << getPedestal(ix,iy)<< endl; // cout <<ix << " " << iy << " " << getPedestal(ix,iy)<< endl;
#ifdef ROOTSPECTRUM #ifdef ROOTSPECTRUM
subtractPedestal(data,ix,iy,cm); subtractPedestal(data,ix,iy,cm);
#endif #endif
}
} }
} }
@ -823,7 +824,8 @@ template <class dataType> class analogDetector {
for (int ix=xmin; ix<xmax; ix++) { for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) { for (int iy=ymin; iy<ymax; iy++) {
val[iy*nx+ix]+=subtractPedestal(data, ix, iy,cm); if (det->isGood(ix,iy))
val[iy*nx+ix]+=subtractPedestal(data, ix, iy,cm);
} }
} }
return val; return val;
@ -950,7 +952,8 @@ template <class dataType> class analogDetector {
for (int ix=xmin; ix<xmax; ix++) { for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) { for (int iy=ymin; iy<ymax; iy++) {
nph[iy*nx+ix]+=getNPhotons(data, ix, iy); if (det->isGood(ix,iy))
nph[iy*nx+ix]+=getNPhotons(data, ix, iy);
} }
} }
return nph; return nph;
@ -1028,8 +1031,11 @@ template <class dataType> class analogDetector {
for (int ix=xmi; ix<xma; ix++) for (int ix=xmi; ix<xma; ix++)
for (int iy=ymi; iy<yma; iy++) for (int iy=ymi; iy<yma; iy++)
if (ix>=0 && ix<nx && iy>=0 && iy<ny) if (det->isGood(ix,iy)) {
val+=getNPhotons(data, ix, iy); if (ix>=0 && ix<nx && iy>=0 && iy<ny)
val+=getNPhotons(data, ix, iy);
}
return val; return val;
}; };
@ -1056,7 +1062,7 @@ template <class dataType> class analogDetector {
} }
}; };
virtual char *getInterpolation(){return NULL;}; // virtual char *getInterpolation(){return NULL;};
/** sets the current frame mode for the detector /** sets the current frame mode for the detector
\param f frame mode to be set \param f frame mode to be set
@ -1095,7 +1101,6 @@ FILE *setFilePointer(FILE *f){myFile=f; return myFile;};
\returns current file pointer \returns current file pointer
*/ */
FILE *getFilePointer(){return myFile;}; FILE *getFilePointer(){return myFile;};
void setMutex(pthread_mutex_t *m){fm=m;};
protected: protected:
slsDetectorData<dataType> *det; /**< slsDetectorData to be used */ slsDetectorData<dataType> *det; /**< slsDetectorData to be used */
@ -1127,7 +1132,6 @@ FILE *getFilePointer(){return myFile;};
TH2F *hs9; TH2F *hs9;
#endif #endif
#endif #endif
pthread_mutex_t *fm;
}; };
#endif #endif

219
slsDetectorCalibration/interpolatingDetector.h
View File

@ -49,17 +49,21 @@ class interpolatingDetector : public singlePhotonDetector {
int nd=100, int nnx=-1, int nny=-1) : int nd=100, int nnx=-1, int nny=-1) :
singlePhotonDetector(d, 3,nsigma,sign, cm, nped, nd, nnx, nny) , interp(inte), id(0) { singlePhotonDetector(d, 3,nsigma,sign, cm, nped, nd, nnx, nny) , interp(inte), id(0) {
//cout << "**"<< xmin << " " << xmax << " " << ymin << " " << ymax << endl; //cout << "**"<< xmin << " " << xmax << " " << ymin << " " << ymax << endl;
fi=new pthread_mutex_t ;
}; };
interpolatingDetector(interpolatingDetector *orig) : singlePhotonDetector(orig) { interpolatingDetector(interpolatingDetector *orig) : singlePhotonDetector(orig) {
if (orig->interp) // if (orig->interp)
interp=(orig->interp)->Clone(); // interp=(orig->interp)->Clone();
else // else
interp=orig->interp;
interp=orig->interp;
id=orig->id; id=orig->id;
fi=orig->fi;
} }
@ -67,7 +71,11 @@ class interpolatingDetector : public singlePhotonDetector {
return new interpolatingDetector(this); 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) { virtual void prepareInterpolation(int &ok) {
/* cout << "*"<< endl; */ /* cout << "*"<< endl; */
@ -81,20 +89,28 @@ class interpolatingDetector : public singlePhotonDetector {
/* sprintf(tit,"/scratch/gmap_%d.tiff",id); */ /* sprintf(tit,"/scratch/gmap_%d.tiff",id); */
/* writeGainMap(tit); */ /* writeGainMap(tit); */
/* } */ /* } */
/* #endif */ /* #endif */
pthread_mutex_lock(fi);
if (interp) if (interp)
interp->prepareInterpolation(ok); interp->prepareInterpolation(ok);
pthread_mutex_unlock(fi);
} }
void clearImage() { void clearImage() {
if (interp) interp->clearInterpolatedImage(); if (interp) {
else singlePhotonDetector::clearImage(); pthread_mutex_lock(fi);
interp->clearInterpolatedImage();
pthread_mutex_unlock(fi);
} else
singlePhotonDetector::clearImage();
}; };
int getImageSize(int &nnx, int &nny, int &ns) { int getImageSize(int &nnx, int &nny, int &ns) {
if (interp) return interp->getImageSize(nnx, nny, ns); if (interp)
else return analogDetector<uint16_t>::getImageSize(nnx, nny, ns); return interp->getImageSize(nnx, nny, ns);
else
return analogDetector<uint16_t>::getImageSize(nnx, nny, ns);
}; };
@ -152,128 +168,7 @@ class interpolatingDetector : public singlePhotonDetector {
return NULL; 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 (iframe<nDark) { */
/* addToPedestal(data); */
/* return 0; */
/* } */
/* newFrame(); */
/* // cout << "********** Data "<< endl; */
/* for (int ix=xmin; ix<xmax; ix++) { */
/* for (int iy=ymin; iy<ymax; iy++) { */
/* max=0; */
/* tl=0; */
/* tr=0; */
/* bl=0; */
/* br=0; */
/* tot=0; */
/* quadTot=0; */
/* quad=UNDEFINED_QUADRANT; */
/* eventMask[iy][ix]=PEDESTAL; */
/* rms=getPedestalRMS(ix,iy); */
/* (clusters+nph)->rms=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)<ny && (ix+ic)>=ix && (ix+ic)<nx) { */
/* val[iy+ir][ix+ic]=subtractPedestal(data,ix+ic,iy+ir); */
/* } */
/* v=&(val[iy+ir][ix+ic]); */
/* tot+=*v; */
/* if (ir<=0 && ic<=0) */
/* bl+=*v; */
/* if (ir<=0 && ic>=0) */
/* br+=*v; */
/* if (ir>=0 && ic<=0) */
/* tl+=*v; */
/* if (ir>=0 && ic>=0) */
/* tr+=*v; */
/* if (*v>max) { */
/* max=*v; */
/* } */
/* if (ir==0 && ic==0) { */
/* if (*v<-nSigma*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) { 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 int_x, int_y;
double eta_x, eta_y; double eta_x, eta_y;
if (interp) { if (interp) {
pthread_mutex_lock(fi);
for (nph=0; nph<nphFrame; nph++) { for (nph=0; nph<nphFrame; nph++) {
if (ff) { if (ff) {
interp->addToFlatField((clusters+nph)->quadTot,(clusters+nph)->quad,(clusters+nph)->get_cluster(),eta_x, eta_y); interp->addToFlatField((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->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); interp->addToImage(int_x, int_y);
} }
} }
pthread_mutex_lock(fi);
} }
return nphFrame; return nphFrame;
}; };
virtual void processData(char *data, int *val=NULL) { virtual void processData(char *data, int *val=NULL) {
if (interp){ switch (dMode) {
switch(fMode) { case eAnalog:
case ePedestal: analogDetector<uint16_t>::processData(data,val);
addToPedestal(data);
break;
case eFlat:
addFrame(data,val,1);
break; break;
case ePhotonCounting:
singlePhotonDetector::processData(data,val);
default: 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 *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: protected:
slsInterpolation *interp; slsInterpolation *interp;
int id; int id;
pthread_mutex_t *fi;
}; };

24
slsDetectorCalibration/interpolations/eta2InterpolationBase.h
View File

@ -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) { 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; // cout << "e2ib " << nb << " " << emin << " " << emax << endl;
if (nbeta<=0) { /* if (nbeta<=0) { */
nbeta=nSubPixels*10; /* nbeta=nSubPixels*10; */
} /* } */
if (etamin>=etamax) { if (etamin>=etamax) {
etamin=-1; etamin=-1;
etamax=2; etamax=2;
cout << ":" <<endl; // cout << ":" <<endl;
} }
etastep=(etamax-etamin)/nbeta; etastep=(etamax-etamin)/nbeta;
#ifdef MYROOT1 #ifdef MYROOT1
@ -34,12 +34,12 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
hhy=new TH2D("hhy","hhy",nbeta,etamin,etamax,nbeta,etamin,etamax); hhy=new TH2D("hhy","hhy",nbeta,etamin,etamax,nbeta,etamin,etamax);
#endif #endif
#ifndef MYROOT1 #ifndef MYROOT1
delete [] heta; /* delete [] heta; */
delete [] hhx; /* delete [] hhx; */
delete [] hhy; /* delete [] hhy; */
heta=new int[nbeta*nbeta]; /* heta=new int[nbeta*nbeta]; */
hhx=new float[nbeta*nbeta]; /* hhx=new float[nbeta*nbeta]; */
hhy=new float[nbeta*nbeta]; /* hhy=new float[nbeta*nbeta]; */
#endif #endif
@ -48,7 +48,7 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
eta2InterpolationBase(eta2InterpolationBase *orig): etaInterpolationBase(orig){ }; eta2InterpolationBase(eta2InterpolationBase *orig): etaInterpolationBase(orig){ };
virtual eta2InterpolationBase* Clone()=0;/* { /* virtual eta2InterpolationBase* Clone()=0; {
return new eta2InterpolationBase(this); return new eta2InterpolationBase(this);
@ -56,8 +56,6 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
*/ */
////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////
//////////// /*It return position hit for the event in input */ ////////////// //////////// /*It return position hit for the event in input */ //////////////
virtual void getInterpolatedPosition(int x, int y, int *data, double &int_x, double &int_y) virtual void getInterpolatedPosition(int x, int y, int *data, double &int_x, double &int_y)

20
slsDetectorCalibration/interpolations/eta3InterpolationBase.h
View File

@ -15,9 +15,9 @@ class eta3InterpolationBase : public virtual etaInterpolationBase {
public: public:
eta3InterpolationBase(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(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 << "e3ib " << nb << " " << emin << " " << emax << endl; // cout << "e3ib " << nb << " " << emin << " " << emax << endl;
if (nbeta<=0) { /* if (nbeta<=0) { */
nbeta=nSubPixels*10; /* nbeta=nSubPixels*10; */
} /* } */
if (etamin>=etamax) { if (etamin>=etamax) {
etamin=-1; etamin=-1;
etamax=1; etamax=1;
@ -33,13 +33,13 @@ class eta3InterpolationBase : public virtual etaInterpolationBase {
hhy=new TH2D("hhy","hhy",nbeta,etamin,etamax,nbeta,etamin,etamax); hhy=new TH2D("hhy","hhy",nbeta,etamin,etamax,nbeta,etamin,etamax);
#endif #endif
#ifndef MYROOT1 #ifndef MYROOT1
delete [] heta; /* delete [] heta; */
delete [] hhx; /* delete [] hhx; */
delete [] hhy; /* delete [] hhy; */
heta=new int[nbeta*nbeta]; /* heta=new int[nbeta*nbeta]; */
hhx=new float[nbeta*nbeta]; /* hhx=new float[nbeta*nbeta]; */
hhy=new float[nbeta*nbeta]; /* hhy=new float[nbeta*nbeta]; */
#endif #endif
// cout << nbeta << " " << etamin << " " << etamax << endl; // cout << nbeta << " " << etamin << " " << etamax << endl;
@ -47,7 +47,7 @@ class eta3InterpolationBase : public virtual etaInterpolationBase {
eta3InterpolationBase(eta3InterpolationBase *orig): etaInterpolationBase(orig){ }; eta3InterpolationBase(eta3InterpolationBase *orig): etaInterpolationBase(orig){ };
virtual eta3InterpolationBase* Clone()=0; /* virtual eta3InterpolationBase* Clone()=0; */

36
slsDetectorCalibration/interpolations/etaInterpolationBase.h
View File

@ -19,11 +19,11 @@ class etaInterpolationBase : public slsInterpolation {
// cout << "eb " << nb << " " << emin << " " << emax << endl; // cout << "eb " << nb << " " << emin << " " << emax << endl;
// cout << nb << " " << etamin << " " << etamax << endl; // cout << nb << " " << etamin << " " << etamax << endl;
if (nbeta<=0) { if (nbeta<=0) {
cout << "aaa:" <<endl; //cout << "aaa:" <<endl;
nbeta=nSubPixels*10; nbeta=nSubPixels*10;
} }
if (etamin>=etamax) { if (etamin>=etamax) {
cout << "aaa:" <<endl; // cout << "aaa:" <<endl;
etamin=-1; etamin=-1;
etamax=2; etamax=2;
} }
@ -68,12 +68,31 @@ class etaInterpolationBase : public slsInterpolation {
}; };
virtual etaInterpolationBase* Clone()=0;/*{
/* virtual etaInterpolationBase* Clone()=0; {
return new etaInterpolationBase(this); return new etaInterpolationBase(this);
};*/
virtual void resetFlatField() {
#ifndef MYROOT1
for (int ibx=0; ibx<nbeta*nbeta; ibx++) {
heta[ibx]=0;
hhx[ibx]=0;
hhy[ibx]=0;
}
#endif
#ifdef MYROOT1
heta->Reset();
hhx->Reset();
hhy->Reset();
#endif
}; };
*/
#ifdef MYROOT1 #ifdef MYROOT1
TH2D *setEta(TH2D *h, int nb=-1, double emin=1, double emax=0) 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); return setEta(h, nb, emin, emax);
}; };
int *getFlatField(){return setEta(NULL);}; int *getFlatField(){return setEta(NULL);};
int *getFlatField(int &nb, double &emin, double &emax){ int *getFlatField(int &nb, double &emin, double &emax){
@ -183,6 +205,8 @@ class etaInterpolationBase : public slsInterpolation {
}; };
#endif #endif

4
slsDetectorCalibration/interpolations/etaInterpolationPosXY.h
View File

@ -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){ 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; // cout << "e2pxy " << nb << " " << emin << " " << emax << endl;
}; };
eta2InterpolationPosXY(eta2InterpolationPosXY *orig): etaInterpolationBase(orig), etaInterpolationPosXY(orig) {}; eta2InterpolationPosXY(eta2InterpolationPosXY *orig): etaInterpolationBase(orig), etaInterpolationPosXY(orig) {};
virtual eta2InterpolationPosXY* Clone() { return new eta2InterpolationPosXY(this);}; virtual eta2InterpolationPosXY* Clone() { return new eta2InterpolationPosXY(this);};
}; };
class eta3InterpolationPosXY : public virtual eta3InterpolationBase, public virtual etaInterpolationPosXY { class eta3InterpolationPosXY : public virtual eta3InterpolationBase, public virtual etaInterpolationPosXY {
public: 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){ 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){

20
slsDetectorCalibration/interpolations/slsInterpolation.h
View File

@ -56,18 +56,30 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
#endif #endif
#ifndef MYROOT1 #ifndef MYROOT1
hint=new int[nSubPixels*nPixelsX*nSubPixels*nPixelsY]; hint=new int[nSubPixels*nPixelsX*nSubPixels*nPixelsY];
memcpy(hint, orig->hint,nSubPixels*nPixelsX*nSubPixels*nPixelsY*sizeof(int)); memcpy(hint, orig->hint,nSubPixels*nPixelsX*nSubPixels*nPixelsY*sizeof(int));
#endif #endif
}; };
virtual int setId(int i) {id=i; return id;}; 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 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) { int getImageSize(int &nnx, int &nny, int &ns) {
nnx=nSubPixels*nPixelsX; nnx=nSubPixels*nPixelsX;
nny=nSubPixels*nPixelsY; 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();}; virtual int *getFlatField(int &nb, double &emin, double &emax){nb=0; emin=0; emax=0; return getFlatField();};
#endif #endif
virtual void resetFlatField()=0;
//virtual void Streamer(TBuffer &b); //virtual void Streamer(TBuffer &b);
static int calcQuad(int *cl, double &sum, double &totquad, double sDum[2][2]){ static int calcQuad(int *cl, double &sum, double &totquad, double sDum[2][2]){

369
slsDetectorCalibration/moenchExecutables/moenchZmqProcess.cpp
View File

@ -17,9 +17,12 @@
#include<iostream> #include<iostream>
//#include "analogDetector.h" //#include "analogDetector.h"
//#include "multiThreadedAnalogDetector.h"
//#include "singlePhotonDetector.h" //#include "singlePhotonDetector.h"
#include "interpolatingDetector.h" //#include "interpolatingDetector.h"
#include "multiThreadedAnalogDetector.h" //#include "multiThreadedCountingDetector.h"
#include "multiThreadedInterpolatingDetector.h"
#include "etaInterpolationPosXY.h"
#include "ansi.h" #include "ansi.h"
#include <iostream> #include <iostream>
using namespace std; using namespace std;
@ -36,27 +39,34 @@ int main(int argc, char *argv[]) {
*/ */
FILE *of=NULL; FILE *of=NULL;
int fifosize=1000; int fifosize=1000;
int nthreads=20; int etabins=1000;//nsubpix*2*100;
double etamin=-1, etamax=2;
// help // help
if (argc < 3 ) { if (argc < 3 ) {
cprintf(RED, "Help: ./trial [receive socket ip] [receive starting port number] [send_socket ip] [send starting port number]\n"); 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; return EXIT_FAILURE;
} }
// receive parameters // receive parameters
bool send = false; bool send = false;
char* socketip=argv[1]; char* socketip=argv[1];
uint32_t portnum = atoi(argv[2]); 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 // send parameters if any
char* socketip2 = 0; char* socketip2 = 0;
uint32_t portnum2 = 0; uint32_t portnum2 = 0;
if (argc > 3) {
send = true;
int ok;
if (argc > 4) {
socketip2 = argv[3]; socketip2 = argv[3];
portnum2 = atoi(argv[4]); portnum2 = atoi(argv[4]);
if (portnum2>0)
send = true;
} }
cout << "\nrx socket ip : " << socketip << cout << "\nrx socket ip : " << socketip <<
"\nrx port num : " << portnum ; "\nrx port num : " << portnum ;
@ -64,26 +74,56 @@ int main(int argc, char *argv[]) {
cout << "\ntx socket ip : " << socketip2 << cout << "\ntx socket ip : " << socketip2 <<
"\ntx port num : " << portnum2; "\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(); //slsDetectorData *det=new moench03T1ZmqDataNew();
moench03T1ZmqDataNew *det=new moench03T1ZmqDataNew(); moench03T1ZmqDataNew *det=new moench03T1ZmqDataNew();
cout << endl << " det" <<endl; cout << endl << " det" <<endl;
int npx, npy; int npx, npy;
det->getDetectorSize(npx, npy); 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<uint16_t> *filter=new analogDetector<uint16_t>(det,1,NULL,1000); //analogDetector<uint16_t> *filter=new analogDetector<uint16_t>(det,1,NULL,1000);
singlePhotonDetector *filter=new singlePhotonDetector(det,3, 5, 1, 0, 1000, 10); //singlePhotonDetector *filter=new singlePhotonDetector(det,3, 5, 1, 0, 1000, 10);
// interpolatingDetector *filter=new interpolatingDetector(det,NULL, 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" <<endl; cout << " filter" <<endl;
char* buff; char* buff;
multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize); // multiThreadedCountingDetector *mt=new multiThreadedCountingDetector(filter,nthreads,fifosize);
cout << " multi" <<endl; // multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize);
multiThreadedInterpolatingDetector *mt=new multiThreadedInterpolatingDetector(filter,nthreads,fifosize);
mt->setFrameMode(eFrame); mt->setFrameMode(eFrame);
mt->StartThreads(); mt->StartThreads();
cout << " start" <<endl;
mt->popFree(buff); mt->popFree(buff);
cout << " pop" <<endl;
ZmqSocket* zmqsocket=NULL; ZmqSocket* zmqsocket=NULL;
@ -120,7 +160,7 @@ int main(int argc, char *argv[]) {
// send socket // send socket
ZmqSocket* zmqsocket2 = 0; ZmqSocket* zmqsocket2 = 0;
cout << "zmq2 " << endl; // cout << "zmq2 " << endl;
if (send) { if (send) {
#ifdef NEWZMQ #ifdef NEWZMQ
// receive socket // receive socket
@ -133,25 +173,29 @@ int main(int argc, char *argv[]) {
#ifdef NEWZMQ #ifdef NEWZMQ
} catch (...) { } catch (...) {
cprintf(RED, "Error: Could not create Zmq socket server on port %d and ip %s\n", portnum2, socketip2); cprintf(RED, "Error: Could not create Zmq socket server on port %d and ip %s\n", portnum2, socketip2);
delete zmqsocket2; // delete zmqsocket2;
delete zmqsocket; // zmqsocket2=NULL;
return EXIT_FAILURE; // delete zmqsocket;
// return EXIT_FAILURE;
send = false;
} }
#endif #endif
#ifndef NEWZMQ #ifndef NEWZMQ
if (zmqsocket2->IsError()) { if (zmqsocket2->IsError()) {
cprintf(RED, "Error: Could not create Zmq socket server on port %d and ip %s\n", portnum2, socketip2); cprintf(RED, "AAA Error: Could not create Zmq socket server on port %d and ip %s\n", portnum2, socketip2);
delete zmqsocket2; // delete zmqsocket2;
delete zmqsocket; //delete zmqsocket;
return EXIT_FAILURE; // return EXIT_FAILURE;
send = false;
} }
#endif #endif
if (zmqsocket2->Connect()) { 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()); zmqsocket2->GetZmqServerAddress());
delete zmqsocket2; // delete zmqsocket2;
return EXIT_FAILURE; send = false;
// return EXIT_FAILURE;
} else } else
printf("Zmq Client at %s\n", zmqsocket2->GetZmqServerAddress()); printf("Zmq Client at %s\n", zmqsocket2->GetZmqServerAddress());
@ -176,7 +220,6 @@ int main(int argc, char *argv[]) {
int *detimage; int *detimage;
int nnx, nny,nns; int nnx, nny,nns;
uint32_t imageSize = 0, nPixelsX = 0, nPixelsY = 0, dynamicRange = 0; uint32_t imageSize = 0, nPixelsX = 0, nPixelsY = 0, dynamicRange = 0;
filter->getImageSize(nnx, nny,nns);
// infinite loop // infinite loop
uint32_t packetNumber = 0; uint32_t packetNumber = 0;
uint64_t bunchId = 0; uint64_t bunchId = 0;
@ -187,7 +230,7 @@ int main(int argc, char *argv[]) {
uint16_t zCoord = 0; uint16_t zCoord = 0;
uint32_t debug = 0; uint32_t debug = 0;
uint32_t dr = 16; uint32_t dr = 16;
int16_t *dout=new int16_t [nnx*nny]; int16_t *dout;//=new int16_t [nnx*nny];
//uint32_t dr = 32; //uint32_t dr = 32;
//int32_t *dout=new int32_t [nnx*nny]; //int32_t *dout=new int32_t [nnx*nny];
uint32_t nSigma=5; uint32_t nSigma=5;
@ -214,6 +257,13 @@ int main(int argc, char *argv[]) {
frameMode fMode; frameMode fMode;
double *ped; double *ped;
filter->getImageSize(nnx, nny,nns);
while(1) { while(1) {
@ -237,48 +287,76 @@ int main(int argc, char *argv[]) {
while (mt->isBusy()) {;}//wait until all data are processed from the queues while (mt->isBusy()) {;}//wait until all data are processed from the queues
detimage=mt->getImage(nnx,nny,nns); if (of) {
if (fMode==ePedestal) { fclose(of);
nns=1; of=NULL;
// cout << "get pedestal " << endl;
ped=mt->getPedestal();
// cout << "got pedestal " << endl;
for (int ix=0; ix<nnx; ix++) {
for (int iy=0; iy<nny; iy++) {
dout[iy*nnx+ix]=ped[(iy)*nnx+ix];
// cout << ped[(iy)*nnx+ix] << " " ;
}
}
// cout << endl ;
// cout << "done " << endl;
} else {
for (int ix=0; ix<nnx; ix++) {
for (int iy=0; iy<nny; iy++) {
for (int isx=0; isx<nns; isx++) {
for (int isy=0; isy<nns; isy++) {
if (isx==0 && isy==0)
dout[iy*nnx+ix]=detimage[(iy+isy)*nnx*nns+ix+isx];
else
dout[iy*nnx+ix]+=detimage[(iy+isy)*nnx*nns+ix+isx];
}
}
if (dout[iy*nnx+ix]<0) dout[iy*nnx+ix]=0;
}
} }
if (fMode==ePedestal) {
sprintf(ofname,"ped_%s_%d.tiff",fname,fileindex);
mt->writePedestal(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; // cout << nns*nnx*nny*nns*dr/8 << " " << length << endl;
if (send) { if (send) {
strcpy(fname,filename.c_str());
#ifdef NEWZMQ if (fMode==ePedestal) {
//zmqsocket2->SendHeaderData (0, false, SLS_DETECTOR_JSON_HEADER_VERSION, dynamicRange, fileindex, nns=1;
// nnx, nny, nns*dynamicRange/8,acqIndex, frameIndex, fname, acqIndex, subFrameIndex, packetNumber,bunchId, timestamp, modId, xCoord, yCoord, zCoord,debug, roundRNumber, detType, version, flippedData, additionalJsonHeader); 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; ix<nnx*nny; ix++) {
// zmqsocket2->SendHeaderData (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; ix<nb*nb; ix++) {
dout[ix]=ff[ix];
}
} else {
detimage=mt->getImage(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; ix<nnx*nns; ix++) {
// for (int iy=0; iy<nny*nns; iy++) {
// for (int isx=0; isx<nns; isx++) {
// for (int isy=0; isy<nns; isy++) {
// if (isx==0 && isy==0)
// dout[iy*nnx+ix]=detimage[(iy+isy)*nnx*nns+ix+isx];
// else
// dout[iy*nnx+ix]+=detimage[(iy+isy)*nnx*nns+ix+isx];
// }
// }
dout[ix]=detimage[ix];
if (dout[ix]<0) dout[ix]=0;
// }
}
}
#ifdef NEWZMQ
zmqsocket2->SendHeaderData (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); zmqsocket2->SendHeaderData (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 #endif
zmqsocket2->SendData((char*)dout,length);//nns*dr/8); zmqsocket2->SendData((char*)dout,length);//nns*dr/8);
cprintf(GREEN, "Sent Data\n"); cprintf(GREEN, "Sent Data\n");
}
zmqsocket2->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION);
cprintf(RED, "Sent Dummy\n");
// stream dummy to socket2 to signal end of acquisition delete [] dout;
if (send) {
zmqsocket2->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION);
cprintf(RED, "Sent Dummy\n");
} }
mt->clearImage(); mt->clearImage();
if (of) {
fclose(of);
of=NULL;
}
newFrame=1; newFrame=1;
continue; //continue to not get out continue; //continue to not get out
} }
#ifdef NEWZMQ #ifdef NEWZMQ
if (newFrame) { if (newFrame) {
// acqIndex, frameIndex, subframeIndex, filename, fileindex // acqIndex, frameIndex, subframeIndex, filename, fileindex
size = doc["size"].GetUint(); size = doc["size"].GetUint();
multisize = size;// * zmqsocket->size(); multisize = size;// * zmqsocket->size();
@ -334,36 +412,38 @@ int main(int argc, char *argv[]) {
detType=doc["detType"].GetUint(); detType=doc["detType"].GetUint();
version=doc["version"].GetUint(); version=doc["version"].GetUint();
dataSize=size;
strcpy(fname,filename.c_str());
cprintf(BLUE, "Header Info:\n" // cprintf(BLUE, "Header Info:\n"
"size: %u\n" // "size: %u\n"
"multisize: %u\n" // "multisize: %u\n"
"dynamicRange: %u\n" // "dynamicRange: %u\n"
"nPixelsX: %u\n" // "nPixelsX: %u\n"
"nPixelsY: %u\n" // "nPixelsY: %u\n"
"currentFileName: %s\n" // "currentFileName: %s\n"
"currentAcquisitionIndex: %lu\n" // "currentAcquisitionIndex: %lu\n"
"currentFrameIndex: %lu\n" // "currentFrameIndex: %lu\n"
"currentFileIndex: %lu\n" // "currentFileIndex: %lu\n"
"currentSubFrameIndex: %u\n" // "currentSubFrameIndex: %u\n"
"xCoordX: %u\n" // "xCoordX: %u\n"
"yCoordY: %u\n" // "yCoordY: %u\n"
"zCoordZ: %u\n" // "zCoordZ: %u\n"
"flippedDataX: %u\n" // "flippedDataX: %u\n"
"packetNumber: %u\n" // "packetNumber: %u\n"
"bunchId: %u\n" // "bunchId: %u\n"
"timestamp: %u\n" // "timestamp: %u\n"
"modId: %u\n" // "modId: %u\n"
"debug: %u\n" // "debug: %u\n"
"roundRNumber: %u\n" // "roundRNumber: %u\n"
"detType: %u\n" // "detType: %u\n"
"version: %u\n", // "version: %u\n",
size, multisize, dynamicRange, nPixelsX, nPixelsY, // size, multisize, dynamicRange, nPixelsX, nPixelsY,
filename.c_str(), acqIndex, // filename.c_str(), acqIndex,
frameIndex, fileindex, subFrameIndex, // frameIndex, fileindex, subFrameIndex,
xCoord, yCoord,zCoord, // xCoord, yCoord,zCoord,
flippedDataX, packetNumber, bunchId, timestamp, modId, debug, roundRNumber, detType, version); // flippedDataX, packetNumber, bunchId, timestamp, modId, debug, roundRNumber, detType, version);
/* Analog detector commands */ /* Analog detector commands */
isPedestal=0; isPedestal=0;
@ -386,7 +466,7 @@ int main(int argc, char *argv[]) {
fMode=eFlat; fMode=eFlat;
isFlat=1; isFlat=1;
} else if (frameMode_s == "newFlatfield") { } else if (frameMode_s == "newFlatfield") {
//mt->resetFlatfield(); mt->resetFlatField();
isFlat=1; isFlat=1;
//cprintf(MAGENTA, "Resetting flatfield\n"); //cprintf(MAGENTA, "Resetting flatfield\n");
fMode=eFlat; fMode=eFlat;
@ -397,16 +477,15 @@ int main(int argc, char *argv[]) {
cprintf(MAGENTA, "%s\n" , frameMode_s.c_str()); cprintf(MAGENTA, "%s\n" , frameMode_s.c_str());
mt->setFrameMode(fMode); mt->setFrameMode(fMode);
threshold=0; // threshold=0;
cprintf(MAGENTA, "Threshold: "); cprintf(MAGENTA, "Threshold: ");
if (doc.HasMember("threshold")) { if (doc.HasMember("threshold")) {
if (doc["threshold"].IsInt()) if (doc["threshold"].IsInt()) {
threshold=doc["threshold"].GetInt(); threshold=doc["threshold"].GetInt();
} mt->setThreshold(threshold);
mt->setThreshold(threshold); }
}
cprintf(MAGENTA, "%d\n", threshold); cprintf(MAGENTA, "%d\n", threshold);
xmin=0; xmin=0;
xmax=npx; 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); mt->setROI(xmin, xmax, ymin, ymax);
// if (doc.HasMember("dynamicRange")) { if (doc.HasMember("dynamicRange")) {
// dr=doc["dynamicRange"].GetUint(); dr=doc["dynamicRange"].GetUint();
// } dr=16;
}
dMode=eAnalog; dMode=eAnalog;
detectorMode_s="analog"; detectorMode_s="analog";
@ -463,28 +543,44 @@ int main(int argc, char *argv[]) {
/* Single Photon Detector commands */ /* Single Photon Detector commands */
nSigma=5;
// if (doc.HasMember("nSigma")) { if (doc.HasMember("nSigma")) {
// nSigma=doc["nSigma"].GetUint(); 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 */ /* interpolating detector commands */
// if (doc.HasMember("nSubPixels")) { if (doc.HasMember("nSubPixels")) {
// nsubPixels=doc["nSubPixels"].GetUint(); if (doc["nSubPixels"].IsUint())
// } nSubPixels=doc["nSubPixels"].GetUint();
mt->setNSubPixels(nSubPixels);
}
//if (doc.HasMember("interpolation")) {
// intMode_s=doc["intepolation"].GetString();
//}
newFrame=0; newFrame=0;
zmqsocket->CloseHeaderMessage(); zmqsocket->CloseHeaderMessage();
@ -514,6 +610,7 @@ int main(int argc, char *argv[]) {
iframe++; iframe++;
} // exiting infinite loop } // exiting infinite loop

209
slsDetectorCalibration/multiThreadedAnalogDetector.h
View File

@ -1,6 +1,8 @@
#ifndef MULTITHREADED_ANALOG_DETECTOR_H #ifndef MULTITHREADED_ANALOG_DETECTOR_H
#define MULTITHREADED_ANALOG_DETECTOR_H #define MULTITHREADED_ANALOG_DETECTOR_H
#define MAXTHREADS 1000
#include <vector> #include <vector>
#include <string> #include <string>
#include <sstream> #include <sstream>
@ -14,11 +16,8 @@
#include <cstdlib> #include <cstdlib>
#include <pthread.h> #include <pthread.h>
//#include "analogDetector.h" #include "analogDetector.h"
#include "singlePhotonDetector.h"
//#include "interpolatingDetector.h"
#include "circularFifo.h" #include "circularFifo.h"
#include "slsInterpolation.h"
#include <unistd.h> #include <unistd.h>
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
@ -57,6 +56,9 @@ public:
return fMode; return fMode;
}; };
virtual double setThreshold(double th) {return det->setThreshold(th);}; 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 void setROI(int xmin, int xmax, int ymin, int ymax) {det->setROI(xmin,xmax,ymin,ymax);};
virtual int setDetectorMode(int dm) { virtual int setDetectorMode(int dm) {
if (dm>=0) { if (dm>=0) {
@ -111,43 +113,7 @@ public:
virtual void clearImage(){det->clearImage();}; 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);}; 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();}; FILE *getFilePointer(){return det->getFilePointer();};
void setMutex(pthread_mutex_t *fmutex){det->setMutex(fmutex);};
private: protected:
analogDetector<uint16_t> *det; analogDetector<uint16_t> *det;
int fMode; int fMode;
int dMode; int dMode;
@ -219,7 +184,7 @@ private:
class multiThreadedAnalogDetector class multiThreadedAnalogDetector
{ {
public: public:
multiThreadedAnalogDetector(analogDetector<uint16_t> *d, int n, int fs=1000) : nThreads(n), ithread(0) { multiThreadedAnalogDetector(analogDetector<uint16_t> *d, int n, int fs=1000) : stop(0), nThreads(n), ithread(0) {
dd[0]=d; dd[0]=d;
if (nThreads==1) if (nThreads==1)
dd[0]->setId(100); dd[0]->setId(100);
@ -239,6 +204,7 @@ public:
image=new int[nn]; image=new int[nn];
ff=NULL; ff=NULL;
ped=NULL; ped=NULL;
cout << "Ithread is " << ithread << endl;
} }
~multiThreadedAnalogDetector() { ~multiThreadedAnalogDetector() {
@ -251,13 +217,15 @@ public:
} }
int setFrameMode(int fm) { int ret; for (int i=0; i<nThreads; i++) ret=dets[i]->setFrameMode(fm); return ret;}; virtual int setFrameMode(int fm) { int ret; for (int i=0; i<nThreads; i++) { ret=dets[i]->setFrameMode(fm);} return ret;};
double setThreshold(int fm) { double ret; for (int i=0; i<nThreads; i++) ret=dets[i]->setThreshold(fm); return ret;}; virtual double setThreshold(int fm) { double ret; for (int i=0; i<nThreads; i++) ret=dets[i]->setThreshold(fm); return ret;};
int setDetectorMode(int dm) { int ret; for (int i=0; i<nThreads; i++) ret=dets[i]->setDetectorMode(dm); return ret;}; virtual int setDetectorMode(int dm) { int ret; for (int i=0; i<nThreads; i++) ret=dets[i]->setDetectorMode(dm); return ret;};
void setROI(int xmin, int xmax, int ymin, int ymax) { for (int i=0; i<nThreads; i++) dets[i]->setROI(xmin, xmax,ymin,ymax);}; virtual void setROI(int xmin, int xmax, int ymin, int ymax) { for (int i=0; i<nThreads; i++) dets[i]->setROI(xmin, xmax,ymin,ymax);};
void newDataSet(){for (int i=0; i<nThreads; i++) dets[i]->newDataSet();};
int *getImage(int &nnx, int &nny, int &ns) { virtual void newDataSet(){for (int i=0; i<nThreads; i++) dets[i]->newDataSet();};
virtual int *getImage(int &nnx, int &nny, int &ns) {
int *img; int *img;
// int nnx, nny, ns; // int nnx, nny, ns;
int nn=dets[0]->getImageSize(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; ii<nThreads; ii++) { for (int ii=0; ii<nThreads; ii++) {
dets[ii]->clearImage(); dets[ii]->clearImage();
@ -315,10 +283,11 @@ public:
virtual void StartThreads() { virtual void StartThreads() {
for (int i=0; i<nThreads; i++) for (int i=0; i<nThreads; i++) {
dets[i]->StartThread(); dets[i]->StartThread();
} }
}
@ -329,7 +298,7 @@ public:
} }
int isBusy() { virtual int isBusy() {
int ret=0, ret1; int ret=0, ret1;
for (int i=0; i<nThreads; i++) { for (int i=0; i<nThreads; i++) {
ret1=dets[i]->isBusy(); ret1=dets[i]->isBusy();
@ -340,125 +309,21 @@ public:
} }
bool pushData(char* &ptr) { virtual bool pushData(char* &ptr) {
dets[ithread]->pushData(ptr); dets[ithread]->pushData(ptr);
} }
bool popFree(char* &ptr) { virtual bool popFree(char* &ptr) {
cout << ithread << endl;
dets[ithread]->popFree(ptr); dets[ithread]->popFree(ptr);
} }
int nextThread() { virtual int nextThread() {
ithread++; ithread++;
if (ithread==nThreads) ithread=0; if (ithread==nThreads) ithread=0;
return ithread; return ithread;
} }
virtual void prepareInterpolation(int &ok){
getFlatField(); //sum up all etas
setFlatField(); //set etas to all detectors
for (int i=0; i<nThreads; i++) {
dets[i]->prepareInterpolation(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; i<nThreads; i++) {
// cout << "mtgff* ff has " << nb << " bins " << endl;
inte=(slsInterpolation*)dets[i]->getInterpolation();
f0=inte->getFlatField();
if (f0) {
// cout << "ff " << i << endl;
for (int ib=0; ib<nb*nb; ib++) {
if (i==0)
ff[ib]=f0[ib];
else
ff[ib]+=f0[ib];
/* if (i==0 && f0[ib]>0) */
/* 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; i<nThreads; i++) {
dets[i]->setFlatField(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; ix<nb*nb; ix++) {
gm[ix]=ff[ix];
}
WriteToTiff(gm,imgname ,nb, nb);
delete [] gm;
} else cout << "Could not allocate float image " << endl;
}
} else
cout << "Detector without flat field! " << endl;
return NULL;
};
void *readFlatField(const char * imgname, int nb=-1, double emin=1, double emax=0){
int* inte=(int*)dets[0]->getFlatField(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; ix<nb*nb; ix++) {
ff[ix]=gm[ix];
}
delete [] gm;
return setFlatField(ff,nb,emin,emax);
} else
cout << "Not interpolating detector! " << endl;
return NULL;
};
virtual double *getPedestal(){ virtual double *getPedestal(){
int nx, ny; int nx, ny;
@ -500,7 +365,7 @@ public:
void *writePedestal(const char * imgname){ virtual void *writePedestal(const char * imgname){
int nx, ny; int nx, ny;
dets[0]->getDetectorSize(nx,ny); dets[0]->getDetectorSize(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; int nx, ny;
dets[0]->getDetectorSize(nx,ny); dets[0]->getDetectorSize(nx,ny);
@ -549,7 +414,7 @@ public:
\param f file pointer \param f file pointer
\returns current file pointer \returns current file pointer
*/ */
FILE *setFilePointer(FILE *f){ virtual FILE *setFilePointer(FILE *f){
for (int i=0; i<nThreads; i++) { for (int i=0; i<nThreads; i++) {
dets[i]->setFilePointer(f); dets[i]->setFilePointer(f);
//dets[i]->setMutex(&fmutex); //dets[i]->setMutex(&fmutex);
@ -557,18 +422,18 @@ public:
return dets[0]->getFilePointer(); return dets[0]->getFilePointer();
}; };
/** gets file pointer where to write the clusters to /** gets file pointer where to write the clusters to
\returns current file pointer \returns current file pointer
*/ */
FILE *getFilePointer(){return dets[0]->getFilePointer();}; virtual FILE *getFilePointer(){return dets[0]->getFilePointer();};
private: protected:
bool stop; bool stop;
const int nThreads; const int nThreads;
threadedAnalogDetector *dets[20]; threadedAnalogDetector *dets[MAXTHREADS];
analogDetector<uint16_t> *dd[20]; analogDetector<uint16_t> *dd[MAXTHREADS];
int ithread; int ithread;
int *image; int *image;
int *ff; int *ff;

55
slsDetectorCalibration/multiThreadedCountingDetector.h Normal file
View File

@ -0,0 +1,55 @@
#ifndef MULTITHREADED_COUNTING_DETECTOR_H
#define MULTITHREADED_COUNTING_DETECTOR_H
#include "singlePhotonDetector.h"
#include "multiThreadedAnalogDetector.h"
//#include <mutex>
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; i<nThreads; i++) ret=((threadedCountingDetector*)dets[i])->setNSigma(n); return ret;};
virtual void setEnergyRange(double emi, double ema) {for (int i=0; i<nThreads; i++) ((threadedCountingDetector*)dets[i])->setEnergyRange(emi,ema);};
};
#endif

136
slsDetectorCalibration/multiThreadedInterpolatingDetector.h Normal file
View File

@ -0,0 +1,136 @@
#ifndef MULTITHREADED_INTERPOLATING_DETECTOR_H
#define MULTITHREADED_INTERPOLATING_DETECTOR_H
#include "interpolatingDetector.h"
#include "multiThreadedCountingDetector.h"
//#include <mutex>
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; i<nThreads; i++) { */
((threadedInterpolatingDetector*)dets[0])->prepareInterpolation(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

483
slsDetectorCalibration/singlePhotonDetector.h
View File

@ -63,6 +63,8 @@ public analogDetector<uint16_t> {
fm=new pthread_mutex_t ;
eventMask=new eventType*[ny]; eventMask=new eventType*[ny];
for (int i=0; i<ny; i++) { for (int i=0; i<ny; i++) {
eventMask[i]=new eventType[nx]; eventMask[i]=new eventType[nx];
@ -114,6 +116,7 @@ public analogDetector<uint16_t> {
// cluster=clusters; // cluster=clusters;
setClusterSize(clusterSize); setClusterSize(clusterSize);
fm=orig->fm;
quad=UNDEFINED_QUADRANT; quad=UNDEFINED_QUADRANT;
tot=0; tot=0;
@ -138,7 +141,7 @@ public analogDetector<uint16_t> {
\param n number of sigma to be set (0 or negative gets) \param n number of sigma to be set (0 or negative gets)
\returns actual number of sigma parameter \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 /** sets/gets cluster size
\param n cluster size to be set, (0 or negative gets). If even is incremented by 1. \param n cluster size to be set, (0 or negative gets). If even is incremented by 1.
@ -217,230 +220,224 @@ public analogDetector<uint16_t> {
} }
for (int ix=xmin; ix<xmax; ix++) { for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) { for (int iy=ymin; iy<ymax; iy++) {
if (det->isGood(ix,iy)) {
val=subtractPedestal(data,ix,iy, cm); val=subtractPedestal(data,ix,iy, cm);
nn=analogDetector<uint16_t>::getNPhotons(data,ix,iy);//val/thr;// nn=analogDetector<uint16_t>::getNPhotons(data,ix,iy);//val/thr;//
if (nn>0) { if (nn>0) {
/* if (nph[ix+nx*iy]>0) { */ nph[ix+nx*iy]+=nn;
/* cout << nph[ix+nx*iy] << " => "; */ rest[iy][ix]=(val-nn*thr);//?+0.5*thr
/* cout << ix << " " << iy << " " << val << " "<< nn << " " << nph[ix+nx*iy]+nn << endl; */ nphFrame+=nn;
/* } */ nphTot+=nn;
nph[ix+nx*iy]+=nn; } else
/* if (nph[ix+nx*iy]>1) { */ rest[iy][ix]=val;
/* 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;
} }
} }
// }
for (int ix=xmin; ix<xmax; ix++) { for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) { for (int iy=ymin; iy<ymax; iy++) {
// for (int ix=clusterSize/2; ix<clusterSize/2-1; ix++) { if (det->isGood(ix,iy)) {
// for (int iy=clusterSizeY/2; iy<ny-clusterSizeY/2; iy++) { eventMask[iy][ix]=PEDESTAL;
eventMask[iy][ix]=PEDESTAL; max=0;
max=0; tl=0;
tl=0; tr=0;
tr=0; bl=0;
bl=0; br=0;
br=0; tot=0;
tot=0; quadTot=0;
quadTot=0;
if (rest[iy][ix]>0.25*thr) {
if (rest[iy][ix]>0.25*thr) { eventMask[iy][ix]=NEIGHBOUR;
eventMask[iy][ix]=NEIGHBOUR; for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) {
for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) { for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) {
for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) { if ((iy+ir)>=0 && (iy+ir)<ny && (ix+ic)>=0 && (ix+ic)<nx) {
if ((iy+ir)>=0 && (iy+ir)<ny && (ix+ic)>=0 && (ix+ic)<nx) { //clusters->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];//clusters->get_data(ic,ir); v=rest[iy+ir][ix+ic];//clusters->get_data(ic,ir);
tot+=v; tot+=v;
if (ir<=0 && ic<=0) if (ir<=0 && ic<=0)
bl+=v; bl+=v;
if (ir<=0 && ic>=0) if (ir<=0 && ic>=0)
br+=v; br+=v;
if (ir>=0 && ic<=0) if (ir>=0 && ic<=0)
tl+=v; tl+=v;
if (ir>=0 && ic>=0) if (ir>=0 && ic>=0)
tr+=v; tr+=v;
if (v>max) { if (v>max) {
max=v; max=v;
}
// if (ir==0 && ic==0) {
//}
} }
// if (ir==0 && ic==0) {
//}
} }
} }
}
if (rest[iy][ix]>=max) {
if (rest[iy][ix]>=max) { if (bl>=br && bl>=tl && bl>=tr) {
if (bl>=br && bl>=tl && bl>=tr) { quad=BOTTOM_LEFT;
quad=BOTTOM_LEFT; quadTot=bl;
quadTot=bl; } else if (br>=bl && br>=tl && br>=tr) {
} else if (br>=bl && br>=tl && br>=tr) { quad=BOTTOM_RIGHT;
quad=BOTTOM_RIGHT; quadTot=br;
quadTot=br; } else if (tl>=br && tl>=bl && tl>=tr) {
} else if (tl>=br && tl>=bl && tl>=tr) { quad=TOP_LEFT;
quad=TOP_LEFT; quadTot=tl;
quadTot=tl; } else if (tr>=bl && tr>=tl && tr>=br) {
} else if (tr>=bl && tr>=tl && tr>=br) { quad=TOP_RIGHT;
quad=TOP_RIGHT;
quadTot=tr; quadTot=tr;
} }
rms=getPedestalRMS(ix,iy); if (nSigma==0) {
tthr=nSigma*rms; tthr=thr;
tthr1=thr;
tthr2=thr;
} else {
rms=getPedestalRMS(ix,iy);
tthr=nSigma*rms;
tthr1=nSigma*sqrt(clusterSize*clusterSizeY)*rms; tthr1=nSigma*sqrt(clusterSize*clusterSizeY)*rms;
tthr2=nSigma*sqrt((clusterSize+1)/2.*((clusterSizeY+1)/2.))*rms; tthr2=nSigma*sqrt((clusterSize+1)/2.*((clusterSizeY+1)/2.))*rms;
if (thr>tthr) tthr=thr-tthr; if (thr>2*tthr) tthr=thr-tthr;
if (thr>tthr1) tthr1=tthr-tthr1; if (thr>2*tthr1) tthr1=tthr-tthr1;
if (thr>tthr2) tthr2=tthr-tthr2; if (thr>2*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++;
} 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); } else return getClusters(data, nph);
} }
return NULL; return NULL;
}; };
/** finds event type for pixel and fills cluster structure. The algorithm loops only if the evenMask for this pixel is still undefined. /* /\** 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 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. /* 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 data pointer to the data */
/param ix pixel x coordinate /* /param ix pixel x coordinate */
/param iy pixel y coordinate /* /param iy pixel y coordinate */
/param cm enable(1)/disable(0) common mode subtraction (if defined). /* /param cm enable(1)/disable(0) common mode subtraction (if defined). */
/returns event type for the given pixel /* /returns event type for the given pixel */
*/ /* *\/ */
eventType getEventType(char *data, int ix, int iy, int cm=0) { /* eventType getEventType(char *data, int ix, int iy, int cm=0) { */
// eventType ret=PEDESTAL; /* // eventType ret=PEDESTAL; */
double max=0, tl=0, tr=0, bl=0,br=0, v; /* double max=0, tl=0, tr=0, bl=0,br=0, v; */
// cout << iframe << endl; /* // cout << iframe << endl; */
int cy=(clusterSizeY+1)/2; /* int cy=(clusterSizeY+1)/2; */
int cs=(clusterSize+1)/2; /* int cs=(clusterSize+1)/2; */
double val; /* double val; */
tot=0; /* tot=0; */
quadTot=0; /* quadTot=0; */
quad=UNDEFINED_QUADRANT; /* quad=UNDEFINED_QUADRANT; */
if (iframe<nDark) { /* if (iframe<nDark) { */
addToPedestal(data, ix,iy); /* addToPedestal(data, ix,iy); */
return UNDEFINED_EVENT; /* return UNDEFINED_EVENT; */
} /* } */
// if (eventMask[iy][ix]==UNDEFINED) { /* // if (eventMask[iy][ix]==UNDEFINED) { */
eventMask[iy][ix]=PEDESTAL; /* eventMask[iy][ix]=PEDESTAL; */
clusters->x=ix; /* clusters->x=ix; */
clusters->y=iy; /* clusters->y=iy; */
clusters->rms=getPedestalRMS(ix,iy); /* clusters->rms=getPedestalRMS(ix,iy); */
clusters->ped=getPedestal(ix,iy, cm); /* clusters->ped=getPedestal(ix,iy, cm); */
for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) { /* for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) { */
for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) { /* for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) { */
if ((iy+ir)>=0 && (iy+ir)<ny && (ix+ic)>=0 && (ix+ic)<nx) { /* if ((iy+ir)>=0 && (iy+ir)<ny && (ix+ic)>=0 && (ix+ic)<nx) { */
v=subtractPedestal(data, ix+ic, iy+ir); /* v=subtractPedestal(data, ix+ic, iy+ir); */
clusters->set_data(v, ic, ir); /* clusters->set_data(v, ic, ir); */
// v=clusters->get_data(ic,ir); /* // v=clusters->get_data(ic,ir); */
tot+=v; /* tot+=v; */
if (ir<=0 && ic<=0) /* if (ir<=0 && ic<=0) */
bl+=v; /* bl+=v; */
if (ir<=0 && ic>=0) /* if (ir<=0 && ic>=0) */
br+=v; /* br+=v; */
if (ir>=0 && ic<=0) /* if (ir>=0 && ic<=0) */
tl+=v; /* tl+=v; */
if (ir>=0 && ic>=0) /* if (ir>=0 && ic>=0) */
tr+=v; /* tr+=v; */
if (v>max) { /* if (v>max) { */
max=v; /* max=v; */
} /* } */
if (ir==0 && ic==0) { /* if (ir==0 && ic==0) { */
if (v<-nSigma*clusters->rms) /* if (v<-nSigma*clusters->rms) */
eventMask[iy][ix]=NEGATIVE_PEDESTAL; /* eventMask[iy][ix]=NEGATIVE_PEDESTAL; */
} /* } */
} /* } */
} /* } */
} /* } */
if (bl>=br && bl>=tl && bl>=tr) { /* if (bl>=br && bl>=tl && bl>=tr) { */
quad=BOTTOM_LEFT; /* quad=BOTTOM_LEFT; */
quadTot=bl; /* quadTot=bl; */
} else if (br>=bl && br>=tl && br>=tr) { /* } else if (br>=bl && br>=tl && br>=tr) { */
quad=BOTTOM_RIGHT; /* quad=BOTTOM_RIGHT; */
quadTot=br; /* quadTot=br; */
} else if (tl>=br && tl>=bl && tl>=tr) { /* } else if (tl>=br && tl>=bl && tl>=tr) { */
quad=TOP_LEFT; /* quad=TOP_LEFT; */
quadTot=tl; /* quadTot=tl; */
} else if (tr>=bl && tr>=tl && tr>=br) { /* } else if (tr>=bl && tr>=tl && tr>=br) { */
quad=TOP_RIGHT; /* quad=TOP_RIGHT; */
quadTot=tr; /* quadTot=tr; */
} /* } */
if (max>nSigma*clusters->rms || tot>sqrt(clusterSizeY*clusterSize)*nSigma*clusters->rms || quadTot>cy*cs*nSigma*clusters->rms) { /* 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) { /* if (clusters->get_data(0,0)>=max) { */
eventMask[iy][ix]=PHOTON_MAX; /* eventMask[iy][ix]=PHOTON_MAX; */
} else { /* } else { */
eventMask[iy][ix]=PHOTON; /* eventMask[iy][ix]=PHOTON; */
} /* } */
} else if (eventMask[iy][ix]==PEDESTAL) { /* } else if (eventMask[iy][ix]==PEDESTAL) { */
if (cm==0) { /* if (cm==0) { */
if (det) /* if (det) */
val=dataSign*det->getValue(data, ix, iy); /* val=dataSign*det->getValue(data, ix, iy); */
else /* else */
val=((double**)data)[iy][ix]; /* val=((double**)data)[iy][ix]; */
addToPedestal(val,ix,iy); /* 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; ix<xmax; ix++) { for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) { for (int iy=ymin; iy<ymax; iy++) {
if (det->isGood(ix,iy)) {
max=0; max=0;
tl=0; tl=0;
tr=0; tr=0;
bl=0; bl=0;
br=0; br=0;
tot=0; tot=0;
quadTot=0; quadTot=0;
quad=UNDEFINED_QUADRANT; quad=UNDEFINED_QUADRANT;
eventMask[iy][ix]=PEDESTAL;
(clusters+nph)->rms=getPedestalRMS(ix,iy); eventMask[iy][ix]=PEDESTAL;
// cluster=clusters+nph;
(clusters+nph)->rms=getPedestalRMS(ix,iy);
// cluster=clusters+nph;
for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) {
for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) { for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) {
for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) {
if ((iy+ir)>=iy && (iy+ir)<ny && (ix+ic)>=ix && (ix+ic)<nx) {
if ((iy+ir)>=iy && (iy+ir)<ny && (ix+ic)>=ix && (ix+ic)<nx) { val[iy+ir][ix+ic]=subtractPedestal(data,ix+ic,iy+ir, cm);
val[iy+ir][ix+ic]=subtractPedestal(data,ix+ic,iy+ir, cm); }
}
v=&(val[iy+ir][ix+ic]);
v=&(val[iy+ir][ix+ic]); tot+=*v;
tot+=*v; if (ir<=0 && ic<=0)
if (ir<=0 && ic<=0) bl+=*v;
bl+=*v; if (ir<=0 && ic>=0)
if (ir<=0 && ic>=0) br+=*v;
br+=*v; if (ir>=0 && ic<=0)
if (ir>=0 && ic<=0) tl+=*v;
tl+=*v; if (ir>=0 && ic>=0)
if (ir>=0 && ic>=0) tr+=*v;
tr+=*v; if (*v>max) {
if (*v>max) { max=*v;
max=*v; }
}
if (ir==0 && ic==0) {
if (ir==0 && ic==0) { if (*v<-nSigma*(clusters+nph)->rms)
if (*v<-nSigma*(clusters+nph)->rms) eventMask[iy][ix]=NEGATIVE_PEDESTAL;
eventMask[iy][ix]=NEGATIVE_PEDESTAL; }
}
} }
} }
if (bl>=br && bl>=tl && bl>=tr) { if (bl>=br && bl>=tl && bl>=tr) {
(clusters+nph)->quad=BOTTOM_LEFT; (clusters+nph)->quad=BOTTOM_LEFT;
(clusters+nph)->quadTot=bl; (clusters+nph)->quadTot=bl;
} else if (br>=bl && br>=tl && br>=tr) { } else if (br>=bl && br>=tl && br>=tr) {
(clusters+nph)->quad=BOTTOM_RIGHT; (clusters+nph)->quad=BOTTOM_RIGHT;
(clusters+nph)->quadTot=br; (clusters+nph)->quadTot=br;
} else if (tl>=br && tl>=bl && tl>=tr) { } else if (tl>=br && tl>=bl && tl>=tr) {
(clusters+nph)->quad=TOP_LEFT; (clusters+nph)->quad=TOP_LEFT;
(clusters+nph)->quadTot=tl; (clusters+nph)->quadTot=tl;
} else if (tr>=bl && tr>=tl && tr>=br) { } else if (tr>=bl && tr>=tl && tr>=br) {
(clusters+nph)->quad=TOP_RIGHT; (clusters+nph)->quad=TOP_RIGHT;
(clusters+nph)->quadTot=tr; (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 (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 (val[iy][ix]>=max) {
eventMask[iy][ix]=PHOTON_MAX; eventMask[iy][ix]=PHOTON_MAX;
(clusters+nph)->tot=tot; (clusters+nph)->tot=tot;
(clusters+nph)->x=ix; (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); (clusters+nph)->set_data(val[iy+ir][ix+ic],ic,ir);
} }
} }
// cout << (clusters+nph)->iframe << " " << ix << " " << nph << " " << tot << " " << (clusters+nph)->quadTot << endl;
good=1; good=1;
if (eMin>0 && tot<eMin) good=0; if (eMin>0 && tot<eMin) good=0;
if (eMax>0 && tot>eMax) good=0; if (eMax>0 && tot>eMax) good=0;
@ -566,13 +562,14 @@ int *getClusters(char *data, int *ph=NULL) {
nph++; nph++;
image[iy*nx+ix]++; image[iy*nx+ix]++;
} }
} else { } else {
eventMask[iy][ix]=PHOTON; eventMask[iy][ix]=PHOTON;
} }
} else if (eventMask[iy][ix]==PEDESTAL) { } else if (eventMask[iy][ix]==PEDESTAL) {
addToPedestal(data,ix,iy,cm); addToPedestal(data,ix,iy,cm);
}
} }
} }
} }
@ -704,6 +701,7 @@ void writeClusters(FILE *f){for (int i=0; i<nphFrame; i++) (clusters+i)->write(f
void setEnergyRange(double emi, double ema){eMin=emi; eMax=ema;}; void setEnergyRange(double emi, double ema){eMin=emi; eMax=ema;};
void getEnergyRange(double &emi, double &ema){emi=eMin; ema=eMax;}; void getEnergyRange(double &emi, double &ema){emi=eMin; ema=eMax;};
void setMutex(pthread_mutex_t *m){fm=m;};
protected: protected:
@ -721,6 +719,7 @@ void writeClusters(FILE *f){for (int i=0; i<nphFrame; i++) (clusters+i)->write(f
int nphTot; int nphTot;
int nphFrame; int nphFrame;
pthread_mutex_t *fm;
}; };