standalone version of angularConversion and postProcessing added

git-svn-id: file:///afs/psi.ch/project/sls_det_software/svn/slsDetectorSoftware@205 951219d9-93cf-4727-9268-0efd64621fa3
This commit is contained in:
bergamaschi 2012-07-17 15:18:31 +00:00
parent aa480795c6
commit 3795dcedfe
5 changed files with 2018 additions and 0 deletions

View File

@ -0,0 +1,68 @@
CFLAGS= -DC_ONLY
#FLAGS= -DVERBOSE -DVERYVERBOSE
INCLUDES= -IcommonFiles -IslsDetector -IMySocketTCP -IusersFunctions -ImultiSlsDetector -IslsDetectorUtils -IslsDetectorCommand -IslsDetectorAnalysis
#EPICSFLAGS=-D EPICS -I/usr/local/epics/base/include/ -I /usr/local/epics/base/include/os/Linux/ -L /usr/local/epics/base/lib/SL5-x86/ -Wl,-R/usr/local/epics/base/lib/SL5-x86 -lca -lCom
CC=gcc
CXX=g++
SRC_CLNT= slsDetectorAnalysis/fileIO.cpp MySocketTCP/MySocketTCP.cpp usersFunctions/usersFunctions.cpp slsDetector/slsDetectorUtils.cpp slsDetector/slsDetectorCommand.cpp slsDetectorAnalysis/AngularConversion_Standalone.cpp slsDetectorAnalysis/energyConversion.cpp slsDetector/slsDetectorActions.cpp slsDetectorAnalysis/postProcessing_Standalone.cpp slsDetector/slsDetector.cpp multiSlsDetector/multiSlsDetector.cpp
OBJS = $(SRC_CLNT:.cpp=.o)
HEADERS = $(SRC_CLNT:.cpp=.h) commonFiles/sls_detector_defs.h slsDetectorAnalysis/detectorData.h slsDetector/slsDetectorBase.h slsDetector/slsDetectorUsers.h multiSlsDetector/multiSlsDetectorCommand.h
SRC_MYTHEN_SVC = mythenDetectorServer/server.c mythenDetectorServer/server_funcs.c mythenDetectorServer/communication_funcs.c mythenDetectorServer/firmware_funcs.c mythenDetectorServer/mcb_funcs.c mythenDetectorServer/trimming_funcs.c
all: package $(SRC_CLNT)
echo "compiling all"
doc: $(SRC_H) $(SRC_CLNT)
doxygen doxy.config
mythenServer: $(SRC_MYTHEN_SVC)
$(CC) $(SRC_MYTHEN_SVC) $(CFLAGS) $(FLAGS) $(INCLUDES) -ImythenDetectorServer -DVIRTUAL -lm -D MCB_FUNCS -DC_ONLY -DVERBOSE
mv a.out mythenServer
picassoServer: $(SRC_MYTHEN_SVC)
$(CC) $(SRC_MYTHEN_SVC) $(CFLAGS) $(FLAGS) $(INCLUDES) -ImythenDetectorServer -D VIRTUAL -lm -DMCB_FUNCS -DPICASSOD -DC_ONLY
mv a.out picassoServer
%.o : %.cpp %.h
$(CXX) -Wall -o $@ -c $< $(INCLUDES) $(FLAGS) $(EPICSFLAGS)
package: $(OBJS)
$(CXX) -shared -Wl,-soname,libSlsDetector.so -o libSlsDetector.so $(OBJS) -lc $(INCLUDES) $(FLAGS) $(EPICSFLAGS)
ar rcs libSlsDetector.a $(OBJS)
clean:
rm -rf libSlsDetector.a libSlsDetector.so core docs/* $(OBJS)
#-------------------------------------------------------------------------------
lib: package
# added install target, HBl
install_lib: lib
$(shell test -d $(DESTDIR) || mkdir -p $(DESTDIR))
cp -P libSlsDetector.so $(DESTDIR)
install_inc:
$(shell test -d $(DESTDIR) || mkdir -p $(DESTDIR))
cp -P $(HEADERS) $(DESTDIR)
install_doc: doc
$(shell test -d $(DESTDIR) || mkdir -p $(DESTDIR))
cp -Pr docs/* $(DESTDIR)

View File

@ -0,0 +1,529 @@
#include "angularConversion.h"
#include <iostream>
#include <fstream>
#include <sstream>
#include <math.h>
#include "usersFunctions.h"
using namespace std;
angularConversion::angularConversion( int *np, float *pos, float *bs, float *fo, float *go): currentPosition(0),
currentPositionIndex(0),
numberOfPositions(np),
detPositions(pos),
binSize(bs),
fineOffset(fo),
globalOffset(go)
{
//angleFunctionPointer=0;
registerAngleFunctionCallback(&defaultAngleFunction);
}
angularConversion::~angularConversion(){
}
// int angularConversion::setAngularConversionPointer(angleConversionConstant *p, int *nm, int nch, int idet) {
// if (p) {
// angOff[idet]=p;
// nMods[idet]=nm;
// nCh[idet]=nch;
// } else {
// angOff[idet]=NULL;
// nMods[idet]=NULL;
// }
// return OK;
// }
float* angularConversion::convertAngles(float pos) {
int imod=0;
float *ang=new float[totalNumberOfChannels];
float enc=pos;
angleConversionConstant *p=NULL;
int ch0=0;
int chlast=chansPerMod[imod];
int nchmod=chansPerMod[imod];
p=angConvs[imod];
if (moveFlag[imod]==0)
enc=0;
else
enc=pos;
for (int ip=0; ip<totalNumberOfChannels; ip++) {
#ifdef VERBOSE
// cout << "ip " << ip << " ch0 " << ch0 << " chlast " << chlast << " imod " << imod << endl;
#endif
if (ip>=chlast) {
imod++;
p=angConvs[imod];
if (moveFlag[imod]==0)
enc=0;
else
enc=pos;
ch0=chlast;
nchmod=chansPerMod[imod];
if (nchmod>0)
chlast+=nchmod;
}
if (p)
ang[ip]=angle(ip-ch0, \
enc, \
(*fineOffset)+(*globalOffset), \
p->r_conversion, \
p->center, \
p->offset, \
p->tilt, \
*angDirection );
}
return ang;
}
//static!
int angularConversion::readAngularConversion(string fname, int nmod, angleConversionConstant *angOff) {
ifstream infile;
string ss;
#ifdef VERBOSE
std::cout<< "Opening file "<< fname << std::endl;
#endif
infile.open(fname.c_str(), ios_base::in);
if (infile.is_open()) {
readAngularConversion(infile, nmod, angOff);
infile.close();
} else {
std::cout<< "Could not open calibration file "<< fname << std::endl;
return -1;
}
return 0;
}
//static
int angularConversion::readAngularConversion( ifstream& infile, int nmod, angleConversionConstant *angOff) {
string str;
int mod;
float center, ecenter;
float r_conv, er_conv;
float off, eoff;
string ss;
int interrupt=0;
int nm=0;
//" module %i center %E +- %E conversion %E +- %E offset %f +- %f \n"
while (infile.good() and interrupt==0) {
getline(infile,str);
#ifdef VERBOSE
cout << "** mod " << nm << " " ;
std::cout<< str << std::endl;
#endif
istringstream ssstr(str);
ssstr >> ss >> mod;
ssstr >> ss >> center;
ssstr >> ss >> ecenter;
ssstr >> ss >> r_conv;
ssstr >> ss >> er_conv;
ssstr >> ss >> off;
ssstr >> ss >> eoff;
if (nm<nmod && nm>=0 ) {
angOff[nm].center=center;
angOff[nm].r_conversion=r_conv;
angOff[nm].offset=off;
angOff[nm].ecenter=ecenter;
angOff[nm].er_conversion=er_conv;
angOff[nm].eoffset=eoff;
} else
break;
//cout << nm<<" " << angOff[nm].offset << endl;
nm++;
if (nm>=nmod)
break;
}
return nm;
}
//static
int angularConversion:: writeAngularConversion(string fname, int nmod, angleConversionConstant *angOff) {
ofstream outfile;
outfile.open (fname.c_str(),ios_base::out);
if (outfile.is_open())
{
writeAngularConversion(outfile, nmod, angOff);
outfile.close();
} else {
std::cout<< "Could not open file " << fname << "for writing"<< std::endl;
return -1;
}
//" module %i center %E +- %E conversion %E +- %E offset %f +- %f \n"
return 0;
}
//static
int angularConversion:: writeAngularConversion(ofstream& outfile, int nmod, angleConversionConstant *angOff) {
for (int imod=0; imod<nmod; imod++) {
outfile << " module " << imod << " center "<< angOff[imod].center<<" +- "<< angOff[imod].ecenter<<" conversion "<< angOff[imod].r_conversion << " +- "<< angOff[imod].er_conversion << " offset "<< angOff[imod].offset << " +- "<< angOff[imod].eoffset << std::endl;
}
return 0;
}
//static
int angularConversion::resetMerging(float *mp, float *mv, float *me, int *mm, int nb) {
#ifdef VERBOSE
cout << "creating merging arrays "<< nb << endl;
#endif
for (int ibin=0; ibin<nb; ibin++) {
mp[ibin]=0;
mv[ibin]=0;
me[ibin]=0;
mm[ibin]=0;
}
return OK;
}
//static
int angularConversion::finalizeMerging(float *mp, float *mv, float *me, int *mm,int nb) {
int np=0;
for (int ibin=0; ibin<nb; ibin++) {
if (mm[ibin]>0) {
#ifdef VERBOSE
cout << "finalize " << ibin << " "<< mm[ibin] << " " << mp[ibin]<< mv[ibin] << me[ibin] << endl;
#endif
mp[np]=mp[ibin]/mm[ibin];
mv[np]=mv[ibin]/mm[ibin];
me[np]=me[ibin]/mm[ibin];
me[np]=sqrt(me[ibin]);
mm[np]=mm[ibin];
np++;
}
}
return np;
}
//static
int angularConversion::addToMerging(float *p1, float *v1, float *e1, float *mp, float *mv,float *me, int *mm, int nchans, float binsize,int nbins, int *badChanMask ) {
float binmi=-180.;
int ibin=0;
if (p1==NULL)
return 0;
if (v1==NULL)
return FAIL;
if (mp==NULL) //can be changed if we want to use a fixed bin algorithm!
return FAIL;
if (mv==NULL)
return FAIL;
if (me==NULL)
return FAIL;
if (mm==NULL)
return FAIL;
if (nchans==0)
return FAIL;
if (binsize<=0)
return FAIL;
if (nbins<=0)
return FAIL;
for (int ip=0; ip<nchans; ip++) {
if (badChanMask) {
if (badChanMask[ip]) {
#ifdef VERBOSE
cout << "channel " << ip << " is bad " << endl;
#endif
continue;
}
}
ibin=(int)((p1[ip]-binmi)/binsize);
if (ibin<nbins && ibin>=0) {
mp[ibin]+=p1[ip];
mv[ibin]+=v1[ip];
if (e1)
me[ibin]+=(e1[ip]*e1[ip]);
else
me[ibin]+=v1[ip];
mm[ibin]++;
#ifdef VERBOSE
cout << "add " << ibin << " "<< mm[ibin] << " " << mp[ibin]<< mv[ibin] << me[ibin] << endl;
#endif
} else
return FAIL;
}
return OK;
}
int angularConversion::deleteMerging() {
if (mergingBins)
delete [] mergingBins;
if (mergingCounts)
delete [] mergingCounts;
if (mergingErrors)
delete [] mergingErrors;
}
int angularConversion::resetMerging() {
getAngularConversionParameter(BIN_SIZE);
mergingBins=new float[nBins];
mergingCounts=new float[nBins];
mergingErrors=new float[nBins];
mergingMultiplicity=new int[nBins];
return resetMerging(mergingBins, mergingCounts, mergingErrors, mergingMultiplicity);
}
int angularConversion::resetMerging(float *mp, float *mv, float *me, int *mm) {
getAngularConversionParameter(BIN_SIZE);
if (nBins)
return resetMerging(mp, mv, me, mm,nBins);
else
return FAIL;
}
int angularConversion::finalizeMerging() {
int np=finalizeMerging(mergingBins, mergingCounts, mergingErrors, mergingMultiplicity);
if (mergingMultiplicity)
delete [] mergingMultiplicity;
return np;
}
int angularConversion::finalizeMerging(float *mp, float *mv, float *me, int *mm) {
if (nBins)
return finalizeMerging(mp, mv, me, mm, nBins);
else
return FAIL;
}
int angularConversion::addToMerging(float *p1, float *v1, float *e1, int *badChanMask ) {
return addToMerging(p1,v1,e1,mergingBins,mergingCounts, mergingErrors, mergingMultiplicity, badChanMask);
}
int angularConversion::addToMerging(float *p1, float *v1, float *e1, float *mp, float *mv,float *me, int *mm, int *badChanMask ) {
int del=0;
if (getAngularConversionParameter(BIN_SIZE)==0){
cout << "no bin size " << endl;
return FAIL;
}
if (nBins==0) {
cout << "no bins " << endl;
return FAIL;
}
if (p1==NULL) {
del=1;
p1=convertAngles();
}
int ret=addToMerging(p1, v1, e1, mp, mv,me, mm,totalNumberOfChannels, *binSize,nBins, badChanMask );
if (del) {
delete [] p1;
p1=NULL;
}
return ret;
}
/**
sets the value of s angular conversion parameter
\param c can be ANGULAR_DIRECTION, GLOBAL_OFFSET, FINE_OFFSET, BIN_SIZE
\param v the value to be set
\returns the actual value
*/
float angularConversion::setAngularConversionParameter(angleConversionParameter c, float v){
switch (c) {
case ANGULAR_DIRECTION:
if (v<0)
*angDirection=-1;
else
*angDirection=1;
return *angDirection;
case GLOBAL_OFFSET:
*globalOffset=v;
return *globalOffset;
case FINE_OFFSET:
*fineOffset=v;
return *fineOffset;
case BIN_SIZE:
if (v>0) {
*binSize=v;
nBins=360./(*binSize);
}
return *binSize;
case MOVE_FLAG:
if (moveFlag) {
if (v>0)
*moveFlag=1;
else if (v==0)
*moveFlag=0;
return *moveFlag;
}
return -1;
default:
return 0;
}
}
/**
returns the value of an angular conversion parameter
\param c can be ANGULAR_DIRECTION, GLOBAL_OFFSET, FINE_OFFSET, BIN_SIZE
\returns the actual value
*/
float angularConversion::getAngularConversionParameter(angleConversionParameter c) {
switch (c) {
case ANGULAR_DIRECTION:
return *angDirection;
case GLOBAL_OFFSET:
return *globalOffset;
case FINE_OFFSET:
return *fineOffset;
case BIN_SIZE:
if (*binSize>0)
nBins=360./(*binSize);
else
nBins=0;
return *binSize;
case MOVE_FLAG:
if (moveFlag)
return *moveFlag;
else
return -1;
default:
return 0;
}
}
int angularConversion::setAngularConversionFile(string fname) {
if (fname=="") {
setAngularCorrectionMask(0);
#ifdef VERBOSE
std::cout << "Unsetting angular conversion" << std::endl;
#endif
} else {
if (fname=="default") {
fname=string(angConvFile);
}
#ifdef VERBOSE
std::cout << "Setting angular conversion to " << fname << std:: endl;
#endif
if (readAngularConversionFile(fname)>=0) {
setAngularCorrectionMask(1);
strcpy(angConvFile,fname.c_str());
}
}
return setAngularCorrectionMask();
}
/*
set positions for the acquisition
\param nPos number of positions
\param pos array with the encoder positions
\returns number of positions
*/
int angularConversion::setPositions(int nPos, float *pos){
if (nPos>=0)
*numberOfPositions=nPos;
for (int ip=0; ip<nPos; ip++)
detPositions[ip]=pos[ip];
return *numberOfPositions;
}
/*
get positions for the acquisition
\param pos array which will contain the encoder positions
\returns number of positions
*/
int angularConversion::getPositions(float *pos){
if (pos) {
for (int ip=0; ip<(*numberOfPositions); ip++)
pos[ip]=detPositions[ip];
}
return *numberOfPositions;
}

View File

@ -0,0 +1,467 @@
#ifndef ANGULARCONVERSION_H
#define ANGULARCONVERSION_H
//#include "slsDetectorBase.h"
#include "sls_detector_defs.h"
#include <string>
#include <fstream>
//float angle(int ichan, float encoder, float totalOffset, float conv_r, float center, float offset, float tilt, int direction)
using namespace std;
/**
@short Angular conversion constants needed for a detector module
*/
typedef struct {
float center; /**< center of the module (channel at which the radius is perpendicular to the module surface) */
float ecenter; /**< error in the center determination */
float r_conversion; /**< detector pixel size (or strip pitch) divided by the diffractometer radius */
float er_conversion; /**< error in the r_conversion determination */
float offset; /**< the module offset i.e. the position of channel 0 with respect to the diffractometer 0 */
float eoffset; /**< error in the offset determination */
float tilt; /**< ossible tilt in the orthogonal direction (unused)*/
float etilt; /**< error in the tilt determination */
} angleConversionConstant;
/**
@short methods to set/unset the angular conversion and merge the data
class containing the methods to set/unset the angular conversion and merge the data
The angular conversion itself is defined by the angle() function defined in usersFunctions.cpp
*/
class angularConversion : public virtual slsDetectorDefs { //: public virtual slsDetectorBase {
public:
/** default constructor */
angularConversion(int*, float*, float*, float*, float*);
/** virtual destructor */
virtual ~angularConversion();
//virtual int readAngularConversion(string fname)=0;
/**
reads an angular conversion file
\param fname file to be read
\param nmod number of modules (maximum) to be read
\param angOff pointer to array of angleConversionConstants
\returns OK or FAIL
*/
static int readAngularConversion(string fname, int nmod, angleConversionConstant *angOff);
/**
reads an angular conversion file
\param ifstream input file stream to be read
\param nmod number of modules (maximum) to be read
\param angOff pointer to array of angleConversionConstants
\returns OK or FAIL
*/
static int readAngularConversion(ifstream& ifs, int nmod, angleConversionConstant *angOff);
/**
writes an angular conversion file
\param fname file to be written
\param nmod number of modules to be written
\param angOff pointer to array of angleConversionConstants
\returns OK or FAIL
*/
static int writeAngularConversion(string fname, int nmod, angleConversionConstant *angOff);
/**
writes an angular conversion file
\param ofstream output file stream
\param nmod number of modules to be written
\param angOff pointer to array of angleConversionConstants
\returns OK or FAIL
*/
static int writeAngularConversion(ofstream& ofs, int nmod, angleConversionConstant *angOff);
/**
sets the arrays of the merged data to 0. NB The array should be created with size nbins >= 360./getBinSize();
\param mp already merged postions
\param mv already merged data
\param me already merged errors (squared sum)
\param mm multiplicity of merged arrays
\param nbins number of bins
\returns OK or FAIL
*/
static int resetMerging(float *mp, float *mv,float *me, int *mm, int nbins);
/**
sets the arrays of the merged data to 0. NB The array should be created with size >= 360./getBinSize();
\param mp already merged postions
\param mv already merged data
\param me already merged errors (squared sum)
\param mm multiplicity of merged arrays
\returns OK or FAIL
*/
int resetMerging(float *mp, float *mv,float *me, int *mm);
/**
creates the arrays for merging the data and sets them to 0.
*/
int resetMerging();
/**
merge dataset
\param p1 angular positions of dataset
\param v1 data
\param e1 errors
\param mp already merged postions
\param mv already merged data
\param me already merged errors (squared sum)
\param mm multiplicity of merged arrays
\param nchans number of channels
\param binsize size of angular bin
\param nb number of angular bins
\param badChanMask badchannelmask (if NULL does not correct for bad channels)
\returns OK or FAIL
*/
static int addToMerging(float *p1, float *v1, float *e1, float *mp, float *mv,float *me, int *mm, int nchans, float binsize,int nb, int *badChanMask );
/**
merge dataset
\param p1 angular positions of dataset
\param v1 data
\param e1 errors
\param mp already merged postions
\param mv already merged data
\param me already merged errors (squared sum)
\param mm multiplicity of merged arrays
\param badChanMask badchannelmask (if NULL does not correct for bad channels)
\returns OK or FAIL
*/
int addToMerging(float *p1, float *v1, float *e1, float *mp, float *mv,float *me, int *mm, int *badChanMask);
/**
merge dataset
\param p1 angular positions of dataset
\param v1 data
\param e1 errors
\param badChanMask badchannelmask (if NULL does not correct for bad channels)
\returns OK or FAIL
*/
int addToMerging(float *p1, float *v1, float *e1,int *badChanMask);
/**
calculates the "final" positions, data value and errors for the merged data
\param mp already merged postions
\param mv already merged data
\param me already merged errors (squared sum)
\param mm multiplicity of merged arrays
\param nb number of bins
\returns FAIL or the number of non empty bins (i.e. points belonging to the pattern)
*/
static int finalizeMerging(float *mp, float *mv,float *me, int *mm, int nb);
/**
calculates the "final" positions, data value and errors for the merged data
\param mp already merged postions
\param mv already merged data
\param me already merged errors (squared sum)
\param mm multiplicity of merged arrays
\returns FAIL or the number of non empty bins (i.e. points belonging to the pattern)
*/
int finalizeMerging(float *mp, float *mv,float *me, int *mm);
/**
calculates the "final" positions, data value and errors for the merged data
\returns FAIL or the number of non empty bins (i.e. points belonging to the pattern)
*/
int finalizeMerging();
/**
set detector global offset
\param f global offset to be set
\returns actual global offset
*/
float setGlobalOffset(float f){return setAngularConversionParameter(GLOBAL_OFFSET,f);};
/**
set detector fine offset
\param f global fine to be set
\returns actual fine offset
*/
float setFineOffset(float f){return setAngularConversionParameter(FINE_OFFSET,f);};
/**
get detector fine offset
\returns actual fine offset
*/
float getFineOffset(){return getAngularConversionParameter(FINE_OFFSET);};
/**
get detector global offset
\returns actual global offset
*/
float getGlobalOffset(){return getAngularConversionParameter(GLOBAL_OFFSET);};
/**
set detector bin size
\param bs bin size to be set
\returns actual bin size
*/
float setBinSize(float bs){if (bs>0) nBins=360/bs; return setAngularConversionParameter(BIN_SIZE,bs);};
/**
get detector bin size
\returns detector bin size used for merging (approx angular resolution)
*/
float getBinSize() {return getAngularConversionParameter(BIN_SIZE);};
/**
get angular direction
\returns actual angular direction (1 is channel number increasing with angle, -1 decreasing)
*/
int getAngularDirection(){return (int)getAngularConversionParameter(ANGULAR_DIRECTION);};
/**
set angular direction
\param d angular direction to be set (1 is channel number increasing with angle, -1 decreasing)
\returns actual angular direction (1 is channel number increasing with angle, -1 decreasing)
*/
int setAngularDirection(int d){return (int)setAngularConversionParameter(ANGULAR_DIRECTION, (float)d);};
/**
\returns number of angular bins in the merging (360./binsize)
*/
int getNumberOfAngularBins(){return nBins;};
/**
set angular conversion parameter
\param c parameter type (globaloffset, fineoffset, binsize, angular direction, move flag)
\param v value to be set
\returns actual value
*/
float setAngularConversionParameter(angleConversionParameter c, float v);
/**
get angular conversion parameter
\param c parameter type (globaloffset, fineoffset, binsize, angular direction, move flag)
\returns actual value
*/
float getAngularConversionParameter(angleConversionParameter c);
/**
set positions for the acquisition
\param nPos number of positions
\param pos array with the encoder positions
\returns number of positions
*/
virtual int setPositions(int nPos, float *pos);
/**
get positions for the acquisition
\param pos array which will contain the encoder positions
\returns number of positions
*/
virtual int getPositions(float *pos=NULL);
/**
deletes the array of merged data
\returns OK
*/
int deleteMerging();
/**
\returns pointer to the array o merged positions
*/
float *getMergedPositions(){return mergingBins;};
/**
\returns pointer to the array of merged counts
*/
float *getMergedCounts(){return mergingCounts;};
/**
\returns pointer to the array of merged errors
*/
float *getMergedErrors(){return mergingErrors;};
/* /\** */
/* reads teh angular conversion file for the (multi)detector and writes it to shared memory */
/* *\/ */
/* virtual int readAngularConversionFile(string fname="")=0; */
/* /\** */
/* get angular conversion */
/* \param direction reference to diffractometer direction */
/* \param angconv array that will be filled with the angular conversion constants */
/* \returns 0 if angular conversion disabled, >0 otherwise */
/* *\/ */
/* virtual int getAngularConversion(int &direction, angleConversionConstant *angconv=NULL)=0; */
/* /\** */
/* pure virtual function */
/* \param file name to be written (nmod and array of angular conversion constants default to the ones ot the slsDetector */
/* *\/ */
/* virtual int writeAngularConversion(string fname)=0; */
/* /\** */
/* \returns number of modules of the (multi)detector */
/* *\/ */
/* virtual int getNMods()=0; */
/* /\** */
/* returns number of channels in the module */
/* \param imod module number */
/* \returns number of channels in the module */
/* *\/ */
/* virtual int getChansPerMod(int imod=0)=0; */
/* /\** */
/* get the angular conversion contant of one modules */
/* \param imod module number */
/* \returns pointer to the angular conversion constant */
/* *\/ */
/* virtual angleConversionConstant *getAngularConversionPointer(int imod=0)=0; */
/* /\** */
/* \param imod module number */
/* \returns move flag of the module (1 encoder is added to the angle, 0 not) */
/* *\/ */
/* virtual int getMoveFlag(int imod)=0; */
/* /\** */
/* enables/disable the angular conversion */
/* \param i 1 sets, 0 unsets,, -1 gets */
/* \returns actual angular conversion flag */
/* *\/ */
/* virtual int setAngularCorrectionMask(int i=-1)=0; */
/**
converts channel number to angle
\param pos encoder position
\returns array of angles corresponding to the channels
*/
float* convertAngles(float pos);
/**
converts channel number to angle for the current encoder position
\returns array of angles corresponding to the channels
*/
float *convertAngles(){return convertAngles(currentPosition);};
/**
returns number of positions
*/
int getNumberOfPositions() {return *numberOfPositions;};
int setTotalNumberOfChannels(int i){if (i>=0) totalNumberOfChannels=i; return totalNumberOfChannels;};
void incrementPositionIndex() {currentPositionIndex++;};
protected:
private:
/** merging bins */
float *mergingBins;
/** merging counts */
float *mergingCounts;
/** merging errors */
float *mergingErrors;
/** merging multiplicity */
int *mergingMultiplicity;
float (*angle)(float, float, float, float, float, float, float, int);
int totalNumberOfChannels;
int moveFlag[MAXMODS*MAXDET];
/** pointer to number of positions for the acquisition*/
int *numberOfPositions;
/** pointer to the detector positions for the acquisition*/
float *detPositions;
/** pointer to angular bin size*/
float *binSize;
int *correctionMask;
int chansPerMod;
int nMod;
angleConversionConstant angConvs[MAXMODS*MAXDET];
int directions[MAXMODS*MAXDET];
/** pointer to beamlien fine offset*/
float *fineOffset;
/** pointer to beamlien global offset*/
float *globalOffset;
/** pointer to beamlien angular direction*/
int *angDirection;
/** pointer to detector move flag (1 moves with encoder, 0 not)*/
// int *moveFlag;
/** number of bins for angular conversion (360./binsize)*/
int nBins;
/**
current position of the detector
*/
float currentPosition;
/**
current position index of the detector
*/
int currentPositionIndex;
/**
returns current position index
*/
int getCurrentPositionIndex() {return currentPositionIndex;};
void registerAngleFunctionCallback(float( *fun)(float, float, float, float, float, float, float, int)) {angle = fun;};
};
#endif

View File

@ -0,0 +1,591 @@
#include "postProcessing.h"
#include "usersFunctions.h"
#include "angularConversion.h"
postProcessing::postProcessing(){
pthread_mutex_t mp1 = PTHREAD_MUTEX_INITIALIZER;
mp=mp1;
pthread_mutex_init(&mp, NULL);
mg=mp1;
pthread_mutex_init(&mg, NULL);
//cout << "reg callback "<< endl;
dataReady = 0;
pCallbackArg = 0;
registerDataCallback(&defaultDataReadyFunc, NULL);
//cout << "done "<< endl;
angConv=new angularConversion(numberOfPositions,detPositions,binSize, fineOffset, globalOffset);
}
int postProcessing::flatFieldCorrect(float datain, float errin, float &dataout, float &errout, float ffcoefficient, float fferr){
float e;
dataout=datain*ffcoefficient;
if (errin==0 && datain>=0)
e=sqrt(datain);
else
e=errin;
if (dataout>0)
errout=sqrt(e*ffcoefficient*e*ffcoefficient+datain*fferr*datain*fferr);
else
errout=1.0;
return 0;
};
int postProcessing::rateCorrect(float datain, float errin, float &dataout, float &errout, float tau, float t){
// float data;
float e;
dataout=(datain*exp(tau*datain/t));
if (errin==0 && datain>=0)
e=sqrt(datain);
else
e=errin;
if (dataout>0)
errout=e*dataout*sqrt((1/(datain*datain)+tau*tau/(t*t)));
else
errout=1.;
return 0;
};
int postProcessing::setBadChannelCorrection(ifstream &infile, int &nbad, int *badlist, int moff){
int interrupt=0;
int ich;
int chmin,chmax;
string str;
nbad=0;
while (infile.good() and interrupt==0) {
getline(infile,str);
#ifdef VERBOSE
std::cout << str << std::endl;
#endif
istringstream ssstr;
ssstr.str(str);
if (ssstr.bad() || ssstr.fail() || infile.eof()) {
interrupt=1;
break;
}
if (str.find('-')!=string::npos) {
ssstr >> chmin ;
ssstr.str(str.substr(str.find('-')+1,str.size()));
ssstr >> chmax;
#ifdef VERBOSE
std::cout << "channels between"<< chmin << " and " << chmax << std::endl;
#endif
for (ich=chmin; ich<=chmax; ich++) {
if (nbad<MAX_BADCHANS) {
badlist[nbad]=ich;
nbad++;
#ifdef VERBOSE
std::cout<< nbad << " Found bad channel "<< ich << std::endl;
#endif
} else
interrupt=1;
}
} else {
ssstr >> ich;
#ifdef VERBOSE
std::cout << "channel "<< ich << std::endl;
#endif
if (nbad<MAX_BADCHANS) {
badlist[nbad]=ich;
nbad++;
#ifdef VERBOSE
std::cout << nbad << " Found bad channel "<< ich << std::endl;
#endif
} else
interrupt=1;
}
}
for (int ich=0; ich<nbad; ich++) {
badlist[ich]=badlist[ich]+moff;
}
return nbad;
}
void postProcessing::processFrame(int *myData, int delflag) {
string fname;
// float *fdata=NULL;
incrementProgress();
/** decode data */
fdata=decodeData(myData, fdata);
fname=createFileName();
//Checking for write flag
if(*correctionMask&(1<<WRITE_FILE))
{
//uses static function?!?!?!?
writeDataFile (fname+string(".raw"),fdata, NULL, NULL, 'i');
}
doProcessing(fdata,delflag, fname);
delete [] myData;
myData=NULL;
fdata=NULL;
#ifdef VERBOSE
cout << "Pop data queue " << *fileIndex << endl;
#endif
pthread_mutex_lock(&mp);
dataQueue.pop(); //remove the data from the queue
queuesize=dataQueue.size();
pthread_mutex_unlock(&mp);
}
void postProcessing::doProcessing(float *lfdata, int delflag, string fname) {
// /** write raw data file */
// if (*correctionMask==0 && delflag==1) {
// // delete [] fdata;
// ;
// } else {
float *rcdata=NULL, *rcerr=NULL;
float *ffcdata=NULL, *ffcerr=NULL;
float *ang=NULL;
// int imod;
int np;
//string fname;
detectorData *thisData;
string ext=".dat";
// fname=createFileName();
/** rate correction */
if (*correctionMask&(1<<RATE_CORRECTION)) {
rcdata=new float[getTotalNumberOfChannels()];
rcerr=new float[getTotalNumberOfChannels()];
rateCorrect(lfdata,NULL,rcdata,rcerr);
delete [] lfdata;
} else {
rcdata=lfdata;
}
lfdata=NULL;
/** flat field correction */
if (*correctionMask&(1<<FLAT_FIELD_CORRECTION)) {
ffcdata=new float[getTotalNumberOfChannels()];
ffcerr=new float[getTotalNumberOfChannels()];
flatFieldCorrect(rcdata,rcerr,ffcdata,ffcerr);
delete [] rcdata;
rcdata=NULL;
if (rcerr)
delete [] rcerr;
rcerr=NULL;
} else {
ffcdata=rcdata;
ffcerr=rcerr;
rcdata=NULL;
rcerr=NULL;
}
// writes angualr converted files
if (*correctionMask!=0) {
if (*correctionMask&(1<< ANGULAR_CONVERSION))
ang=convertAngles();
writeDataFile (fname+ext, ffcdata, ffcerr,ang);
}
if (*correctionMask&(1<< ANGULAR_CONVERSION) && getNumberOfPositions()>0) {
#ifdef VERBOSE
cout << "**************Current position index is " << getCurrentPositionIndex() << endl;
#endif
// if (*numberOfPositions>0) {
if (getCurrentPositionIndex()<=1) {
#ifdef VERBOSE
cout << "reset merging " << endl;
#endif
angConv->resetMerging();
}
#ifdef VERBOSE
cout << "add to merging "<< getCurrentPositionIndex() << endl;
#endif
angConv->addToMerging(ang, ffcdata, ffcerr, badChannelMask );
#ifdef VERBOSE
cout << getCurrentPositionIndex() << " " << getNumberOfPositions() << endl;
#endif
// cout << "lock 1" << endl;
pthread_mutex_lock(&mp);
if ((getCurrentPositionIndex()>=getNumberOfPositions() && posfinished==1 && queuesize==1)) {
#ifdef VERBOSE
cout << "finalize merging " << getCurrentPositionIndex()<< endl;
#endif
np=angConv->finalizeMerging();
/** file writing */
angConv->incrementPositionIndex();
// cout << "unlock 1" << endl;
pthread_mutex_unlock(&mp);
fname=createFileName();
#ifdef VERBOSE
cout << "writing merged data file" << endl;
#endif
writeDataFile (fname+ext,np,angConv->getMergedCounts(), angConv->getMergedErrors(), angConv->getMergedPositions(),'f');
#ifdef VERBOSE
cout << " done" << endl;
#endif
// if (delflag) {
// deleteMerging();
// } else {
thisData=new detectorData(angConv->getMergedCounts(),angConv->getMergedErrors(),angConv->getMergedPositions(),getCurrentProgress(),(fname+ext).c_str(),np);
// // cout << "lock 2" << endl;
// pthread_mutex_lock(&mg);
// finalDataQueue.push(thisData);
// // cout << "unlock 2" << endl;
// pthread_mutex_unlock(&mg);
if (dataReady) {
dataReady(thisData, pCallbackArg);
delete thisData;
}
// }
// cout << "lock 3" << endl;
pthread_mutex_lock(&mp);
}
// cout << "unlock 3" << endl;
pthread_mutex_unlock(&mp);
if (ffcdata)
delete [] ffcdata;
ffcdata=NULL;
if (ffcerr)
delete [] ffcerr;
ffcerr=NULL;
if (ang)
delete [] ang;
ang=NULL;
} else {
// if (delflag) {
// if (ffcdata)
// delete [] ffcdata;
// if (ffcerr)
// delete [] ffcerr;
// if ( ang)
// delete [] ang;
// } else {
thisData=new detectorData(ffcdata,ffcerr,NULL,getCurrentProgress(),(fname+ext).c_str(),getTotalNumberOfChannels());
if (dataReady) {
dataReady(thisData, pCallbackArg);
delete thisData;
}
// pthread_mutex_lock(&mg);
// finalDataQueue.push(thisData);
// pthread_mutex_unlock(&mg);
// }
}
//}
incrementFileIndex();
#ifdef VERBOSE
cout << "fdata is " << fdata << endl;
#endif
}
int postProcessing::fillBadChannelMask() {
int nbad=0;
if (*correctionMask&(1<< DISCARD_BAD_CHANNELS)) {
nbad=getBadChannelCorrection();
#ifdef VERBOSE
cout << "number of bad channels is " << nbad << endl;
#endif
if (nbad>0) {
int *badChansList=new int[nbad];
getBadChannelCorrection(badChansList);
if (badChannelMask)
delete [] badChannelMask;
badChannelMask=new int[getTotalNumberOfChannels()];
#ifdef VERBOSE
cout << " pointer to bad channel mask is " << badChannelMask << endl;
#endif
for (int ichan=0; ichan<getTotalNumberOfChannels(); ichan++)
badChannelMask[ichan]=0;
#ifdef VERBOSE
cout << " badChanMask has be reset" << badChannelMask << endl;
#endif
for (int ichan=0; ichan<nbad; ichan++) {
if (badChansList[ichan]<getTotalNumberOfChannels() && badChansList[ichan]>=0 ) {
if (badChannelMask[badChansList[ichan]]==0)
nbad++;
badChannelMask[badChansList[ichan]]=1;
}
}
delete [] badChansList;
} else {
if (badChannelMask) {
#ifdef VERBOSE
cout << "deleting bad channel mask beacuse number of bad channels is 0" << endl;
#endif
delete [] badChannelMask;
badChannelMask=NULL;
}
}
} else {
#ifdef VERBOSE
cout << "bad channel correction is disabled " << nbad << endl;
#endif
if (badChannelMask) {
#ifdef VERBOSE
cout << "deleting bad channel mask beacuse no bad channel correction is selected" << endl;
#endif
delete [] badChannelMask;
badChannelMask=NULL;
}
}
#ifdef VERBOSE
cout << "number of bad channels is " << nbad << endl;
#endif
return nbad;
}
void* postProcessing::processData(int delflag) {
#ifdef VERBOSE
std::cout<< " processing data - threaded mode " << *threadedProcessing << endl;
#endif
angConv->setTotalNumberOfChannels(getTotalNumberOfChannels());
setTotalProgress();
pthread_mutex_lock(&mp);
queuesize=dataQueue.size();
pthread_mutex_unlock(&mp);
int *myData;
int dum=1;
fdata=NULL;
while(dum | *threadedProcessing) { // ????????????????????????
/* IF THERE ARE DATA PROCESS THEM*/
pthread_mutex_lock(&mp);
while((queuesize=dataQueue.size())>0) {
/** Pop data queue */
myData=dataQueue.front(); // get the data from the queue
pthread_mutex_unlock(&mp);
if (myData) {
processFrame(myData,delflag);
//usleep(1000);
}
pthread_mutex_lock(&mp);
}
pthread_mutex_unlock(&mp);
/* IF THERE ARE NO DATA look if acquisition is finished */
pthread_mutex_lock(&mp);
if (jointhread) {
if (dataQueue.size()==0) {
pthread_mutex_unlock(&mp);
break;
}
pthread_mutex_unlock(&mp);
} else {
pthread_mutex_unlock(&mp);
}
dum=0;
}
if (fdata) {
#ifdef VERBOSE
cout << "delete fdata" << endl;
#endif
delete [] fdata;
#ifdef VERBOSE
cout << "done " << endl;
#endif
}
return 0;
}
int* postProcessing::popDataQueue() {
int *retval=NULL;
if( !dataQueue.empty() ) {
retval=dataQueue.front();
dataQueue.pop();
}
return retval;
}
detectorData* postProcessing::popFinalDataQueue() {
detectorData *retval=NULL;
pthread_mutex_unlock(&mg);
if( !finalDataQueue.empty() ) {
retval=finalDataQueue.front();
finalDataQueue.pop();
}
pthread_mutex_unlock(&mg);
return retval;
}
void postProcessing::resetDataQueue() {
int *retval=NULL;
while( !dataQueue.empty() ) {
retval=dataQueue.front();
dataQueue.pop();
delete [] retval;
}
}
void postProcessing::resetFinalDataQueue() {
detectorData *retval=NULL;
pthread_mutex_lock(&mg);
while( !finalDataQueue.empty() ) {
retval=finalDataQueue.front();
finalDataQueue.pop();
delete retval;
}
pthread_mutex_unlock(&mg);
}
void postProcessing::startThread(int delflag) {
pthread_attr_t tattr;
int ret;
sched_param param, mparam;
int policy= SCHED_OTHER;
// set the priority; others are unchanged
//newprio = 30;
mparam.sched_priority =1;
param.sched_priority =1;
/* Initialize and set thread detached attribute */
pthread_attr_init(&tattr);
pthread_attr_setdetachstate(&tattr, PTHREAD_CREATE_JOINABLE);
// param.sched_priority = 5;
// scheduling parameters of main thread
ret = pthread_setschedparam(pthread_self(), policy, &mparam);
//#ifdef VERBOSE
// printf("current priority is %d\n",param.sched_priority);
//#endif
if (delflag)
ret = pthread_create(&dataProcessingThread, &tattr,startProcessData, (void*)this);
else
ret = pthread_create(&dataProcessingThread, &tattr,startProcessDataNoDelete, (void*)this);
pthread_attr_destroy(&tattr);
// scheduling parameters of target thread
ret = pthread_setschedparam(dataProcessingThread, policy, &param);
}

View File

@ -0,0 +1,363 @@
#ifndef POSTPROCESSING_H
#define POSTPROCESSING_H
#include "detectorData.h"
#include "fileIO.h"
#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cstring>
#include <string>
#include <sstream>
#include <queue>
#include <math.h>
class angularConversion;
using namespace std;
#define MAX_BADCHANS 2000
#define defaultTDead {170,90,750} /**< should be changed in order to have it separate for the different detector types */
/**
@short methods for data postprocessing
(including thread for writing data files and plotting in parallel with the acquisition)
*/
class postProcessing : public virtual fileIO {
//: public angularConversion, public fileIO
public:
postProcessing();
virtual ~postProcessing(){};
/**
get bad channels correction
\param bad pointer to array that if bad!=NULL will be filled with the bad channel list
\returns 0 if bad channel disabled or no bad channels, >0 otherwise
*/
virtual int getBadChannelCorrection(int *bad=NULL)=0;
/**
get flat field corrections
\param corr if !=NULL will be filled with the correction coefficients
\param ecorr if !=NULL will be filled with the correction coefficients errors
\returns 0 if ff correction disabled, >0 otherwise
*/
virtual int getFlatFieldCorrection(float *corr=NULL, float *ecorr=NULL)=0;
/**
set flat field corrections
\param corr if !=NULL the flat field corrections will be filled with corr (NULL usets ff corrections)
\param ecorr if !=NULL the flat field correction errors will be filled with ecorr (1 otherwise)
\returns 0 if ff correction disabled, >0 otherwise
*/
virtual int setFlatFieldCorrection(float *corr, float *ecorr=NULL)=0;
/**
set bad channels correction
\param fname file with bad channel list ("" disable)
\returns 0 if bad channel disabled, >0 otherwise
*/
virtual int setBadChannelCorrection(string fname="")=0;
static int setBadChannelCorrection(ifstream &infile, int &nbad, int *badlist, int moff=0);
/**
set bad channels correction
\param fname file with bad channel list ("" disable)
\param nbad reference to number of bad channels
\param badlist array of badchannels
\returns 0 if bad channel disabled, >0 otherwise
*/
virtual int setBadChannelCorrection(string fname, int &nbad, int *badlist, int off=0)=0;
/**
set bad channels correction
\param nch number of bad channels
\param chs array of channels
\param ff 0 if normal bad channels, 1 if ff bad channels
\returns 0 if bad channel disabled, >0 otherwise
*/
virtual int setBadChannelCorrection(int nch, int *chs, int ff=0)=0;
/**
flat field correct data
\param datain data
\param errin error on data (if<=0 will default to sqrt(datain)
\param dataout corrected data
\param errout error on corrected data
\param ffcoefficient flat field correction coefficient
\param fferr erro on ffcoefficient
\returns 0
*/
static int flatFieldCorrect(float datain, float errin, float &dataout, float &errout, float ffcoefficient, float fferr);
/**
rate correct data
\param datain data
\param errin error on data (if<=0 will default to sqrt(datain)
\param dataout corrected data
\param errout error on corrected data
\param tau dead time 9in ns)
\param t acquisition time (in ns)
\returns 0
*/
static int rateCorrect(float datain, float errin, float &dataout, float &errout, float tau, float t);
int enableWriteToFile(int i=-1) {if (i>0) ((*correctionMask)|=(1<<WRITE_FILE)); if(i==0) ((*correctionMask)&=~(1<< WRITE_FILE)); return ((*correctionMask)&(1<< WRITE_FILE));};
int setAngularCorrectionMask(int i=-1){if (i==0) (*correctionMask)&=~(1<< ANGULAR_CONVERSION); if (i>0) (*correctionMask)|=(1<< ANGULAR_CONVERSION); return ((*correctionMask)&(1<< ANGULAR_CONVERSION));};
int enableAngularConversion(int i=-1) {if (i>0) return setAngularConversionFile("default"); if (i==0) return setAngularConversionFile(""); return setAngularCorrectionMask();};
int enableBadChannelCorrection(int i=-1) {if (i>0) return setBadChannelCorrection("default"); if (i==0) return setBadChannelCorrection(""); return ((*correctionMask)&(1<< DISCARD_BAD_CHANNELS));};
/** returns the bad channel list file */
string getBadChannelCorrectionFile() {if ((*correctionMask)&(1<< DISCARD_BAD_CHANNELS)) return string(badChanFile); else return string("none");};
/**
get flat field corrections file directory
\returns flat field correction file directory
*/
string getFlatFieldCorrectionDir(){return string(flatFieldDir);};
/**
set flat field corrections file directory
\param flat field correction file directory
\returns flat field correction file directory
*/
string setFlatFieldCorrectionDir(string dir){strcpy(flatFieldDir,dir.c_str()); return string(flatFieldDir);};
/**
get flat field corrections file name
\returns flat field correction file name
*/
string getFlatFieldCorrectionFile(){ if ((*correctionMask)&(1<<FLAT_FIELD_CORRECTION)) return string(flatFieldFile); else return string("none");};
/**
set/get if the data processing and file writing should be done by a separate thread
s
\param b 0 sequencial data acquisition and file writing, 1 separate thread, -1 get
\returns thread flag
*/
int setThreadedProcessing(int b=-1) {if (b>=0) *threadedProcessing=b; return *threadedProcessing;};
/** processes the data
\param delflag 0 leaves the data in the final data queue
\returns nothing
*/
void *processData(int delflag);
/** processes the data
\param delflag 0 leaves the data in the final data queue
\returns nothing
*/
void processFrame(int* myData, int delflag);
/** processes the data
\param delflag 0 leaves the data in the final data queue
\returns nothing
*/
void doProcessing(float* myData, int delflag, string fname);
/**
pops the data from the data queue
\returns pointer to the popped data or NULL if the queue is empty.
\sa dataQueue
*/
int* popDataQueue();
/**
pops the data from thepostprocessed data queue
\returns pointer to the popped data or NULL if the queue is empty.
\sa finalDataQueue
*/
detectorData* popFinalDataQueue();
/**
resets the raw data queue
\sa dataQueue
*/
void resetDataQueue();
/**
resets the postprocessed data queue
\sa finalDataQueue
*/
void resetFinalDataQueue();
int fillBadChannelMask();
virtual int rateCorrect(float*, float*, float*, float*)=0;
virtual int flatFieldCorrect(float*, float*, float*, float*)=0;
void registerDataCallback(int( *userCallback)(detectorData*, void*), void *pArg) {dataReady = userCallback; pCallbackArg = pArg;};
/**
sets the angular conversion file
\param fname file to read
\returns angular conversion flag
*/
int setAngularConversionFile(string fname);
/**
returns the angular conversion file
*/
string getAngularConversionFile(){if (setAngularCorrectionMask()) return string(angConvFile); else return string("none");};
protected:
int *threadedProcessing;
int *correctionMask;
char *flatFieldDir;
char *flatFieldFile;
char *badChanFile;
int *nBadChans;
int *badChansList;
int *nBadFF;
int *badFFList;
/** pointer to angular conversion file name*/
char *angConvFile;
/** mutex to synchronize main and data processing threads */
pthread_mutex_t mp;
/** mutex to synchronizedata processing and plotting threads */
pthread_mutex_t mg;
/** sets when the acquisition is finished */
int jointhread;
/** sets when the position is finished */
int posfinished;
/**
data queue
*/
queue<int*> dataQueue;
/**
queue containing the postprocessed data
*/
queue<detectorData*> finalDataQueue;
/** data queue size */
int queuesize;
/**
start data processing thread
*/
void startThread(int delflag=1); //
/** the data processing thread */
pthread_t dataProcessingThread;
/** pointer to bad channel mask 0 is channel is good 1 if it is bad \sa fillBadChannelMask() */
int *badChannelMask;
/**
I0 measured
*/
float currentI0;
float *fdata;
int (*dataReady)(detectorData*,void*);
void *pCallbackArg;
private:
angularConversion *angConv;
};
static void* startProcessData(void *n){\
postProcessing *myDet=(postProcessing*)n;\
myDet->processData(1);\
pthread_exit(NULL);\
};
static void* startProcessDataNoDelete(void *n){\
postProcessing *myDet=(postProcessing*)n;\
myDet->processData(0);\
pthread_exit(NULL);\
};
#endif