postprocessing based on external functions - can work with f90 interface

git-svn-id: file:///afs/psi.ch/project/sls_det_software/svn/slsDetectorSoftware@285 951219d9-93cf-4727-9268-0efd64621fa3
This commit is contained in:
bergamaschi
2012-10-08 09:34:14 +00:00
parent e5c5b76236
commit 05181fa618
27 changed files with 1291 additions and 864 deletions

View File

@ -4,6 +4,10 @@
class angleConversionConstant {
public:
angleConversionConstant(){};
angleConversionConstant(double c, double r, double o, double t){center=c; r_conversion=r; offset=o; tilt=t;};
//typedef struct {
double center; /**< center of the module (channel at which the radius is perpendicular to the module surface) */
double ecenter; /**< error in the center determination */
@ -20,7 +24,7 @@ class angleConversionConstant {
double getOffset(){return offset;};
double getTilt(){return tilt;};
int setAngConvConstant(angleConversionConstant *acc) {\
void setAngConvConstant(angleConversionConstant *acc) {\
center=acc->center; \
ecenter=acc->ecenter; \
r_conversion=acc->r_conversion; \

View File

@ -190,6 +190,18 @@ double angularConversion::setAngularConversionParameter(angleConversionParameter
return *moveFlag;
}
return -1;
case SAMPLE_X:
if (sampleDisplacement) {
sampleDisplacement[X]=v;
return sampleDisplacement[X];
}
return 0;
case SAMPLE_Y:
if (sampleDisplacement) {
sampleDisplacement[Y]=v;
return sampleDisplacement[Y];
}
return 0;
default:
return 0;
}

View File

@ -345,6 +345,8 @@ class angularConversion : public virtual slsDetectorBase, public angularConversi
double *sampleDisplacement;
/**
current position of the detector
*/
@ -367,6 +369,8 @@ class angularConversion : public virtual slsDetectorBase, public angularConversi
*/
int getCurrentPositionIndex() {return currentPositionIndex;};
void incrementPositionIndex() {currentPositionIndex++;};
void resetPositionIndex() {currentPositionIndex=0;};
@ -385,7 +389,6 @@ class angularConversion : public virtual slsDetectorBase, public angularConversi
/** merging bins */
double *mergingBins;

View File

@ -30,7 +30,7 @@ double* angularConversionStatic::convertAngles(double pos, int nch, int *chansPe
angleConversionConstant *p=NULL;
int ch0=0;
int chlast=chansPerMod[0];
int chlast=chansPerMod[0]-1;
int nchmod=chansPerMod[0];
p=angOff[imod];
if (mF[imod]==0)
@ -42,7 +42,7 @@ double* angularConversionStatic::convertAngles(double pos, int nch, int *chansPe
#ifdef VERBOSE
// cout << "ip " << ip << " ch0 " << ch0 << " chlast " << chlast << " imod " << imod << endl;
#endif
if (ip>=chlast) {
if (ip>chlast) {
imod++;
p=angOff[imod];
if (mF[imod]==0)
@ -50,10 +50,10 @@ double* angularConversionStatic::convertAngles(double pos, int nch, int *chansPe
else
enc=pos;
ch0=chlast;
ch0=chlast+1;
nchmod=chansPerMod[imod];
if (nchmod>0)
chlast+=nchmod;
chlast=ch0+nchmod-1;
}
if (p)
@ -69,6 +69,72 @@ double* angularConversionStatic::convertAngles(double pos, int nch, int *chansPe
return ang;
}
double angularConversionStatic::convertAngle(double pos, int ich, int *chansPerMod, angleConversionConstant **angOff, int *mF, double fo, double go, int angdir) {
int imod=0;
double ang;
double enc=0, trans=0;
angleConversionConstant *p=NULL;
int ch0=0;
int chlast=chansPerMod[0]-1;
int nchmod=chansPerMod[0];
while (ich>chlast) {
imod++;
ch0=chlast+1;
nchmod=chansPerMod[imod];
chlast=ch0+nchmod-1;
}
p=angOff[imod];
switch (mF[imod]) {
case 0:
enc=0;
trans=0;
break;
case 1:
enc=pos;
trans=0;
break;
case -1:
enc=-pos;
trans=0;
break;
case 2:
enc=0;
trans=pos;
break;
case -2:
enc=0;
trans=-pos;
break;
default:
enc=0;
trans=0;
}
if (p)
ang=angle(ich-ch0, \
enc, \
fo+go, \
p->r_conversion, \
p->center, \
p->offset, \
trans, \
angdir );
return ang;
}
//static!
@ -265,6 +331,53 @@ int angularConversionStatic::addToMerging(double *p1, double *v1, double *e1, d
}
return slsDetectorDefs::OK;
}
int angularConversionStatic::addPointToMerging(double p1, double v1, double e1, double *mp, double *mv,double *me, int *mm, double binsize,int nbins) {
double binmi=-180.;
int ibin=0;
if (mp==NULL) //can be changed if we want to use a fixed bin algorithm!
return slsDetectorDefs::FAIL;
if (mv==NULL)
return slsDetectorDefs::FAIL;
if (me==NULL)
return slsDetectorDefs::FAIL;
if (mm==NULL)
return slsDetectorDefs::FAIL;
if (binsize<=0)
return slsDetectorDefs::FAIL;
if (nbins<=0)
return slsDetectorDefs::FAIL;
ibin=(int)((p1-binmi)/binsize);
if (ibin<nbins && ibin>=0) {
mp[ibin]+=p1;
mv[ibin]+=v1;
if (e1)
me[ibin]+=(e1*e1);
else
me[ibin]+=v1;
mm[ibin]++;
#ifdef VERBOSE
cout << "add " << ibin << " "<< mm[ibin] << " " << mp[ibin]<< mv[ibin] << me[ibin] << endl;
#endif
} else
return slsDetectorDefs::FAIL;
return slsDetectorDefs::OK;
}

View File

@ -4,7 +4,9 @@
#include <fstream>
#include <sstream>
#include <math.h>
#include "angleConversionConstant.h"
#include "sls_detector_defs.h"
#include "angleFunction.h"
using namespace std;
@ -180,7 +182,7 @@ int angularConversionStatic::resetMerging(double *mp, double *mv, double *me, in
me[ibin]=0;
mm[ibin]=0;
}
return OK;
return slsDetectorDefs::OK;
}
@ -214,25 +216,25 @@ int angularConversionStatic::addToMerging(double *p1, double *v1, double *e1, d
if (p1==NULL)
return 0;
if (v1==NULL)
return FAIL;
return slsDetectorDefs::FAIL;
if (mp==NULL) //can be changed if we want to use a fixed bin algorithm!
return FAIL;
return slsDetectorDefs::FAIL;
if (mv==NULL)
return FAIL;
return slsDetectorDefs::FAIL;
if (me==NULL)
return FAIL;
return slsDetectorDefs::FAIL;
if (mm==NULL)
return FAIL;
return slsDetectorDefs::FAIL;
if (nchans==0)
return FAIL;
return slsDetectorDefs::FAIL;
if (binsize<=0)
return FAIL;
return slsDetectorDefs::FAIL;
if (nbins<=0)
return FAIL;
return slsDetectorDefs::FAIL;
for (int ip=0; ip<nchans; ip++) {
if (badChanMask) {
@ -259,10 +261,10 @@ int angularConversionStatic::addToMerging(double *p1, double *v1, double *e1, d
cout << "add " << ibin << " "<< mm[ibin] << " " << mp[ibin]<< mv[ibin] << me[ibin] << endl;
#endif
} else
return FAIL;
return slsDetectorDefs::FAIL;
}
return OK;
return slsDetectorDefs::OK;
}

View File

@ -116,7 +116,26 @@ class angularConversionStatic
\returns OK or FAIL
*/
static int addToMerging(double *p1, double *v1, double *e1, double *mp, double *mv,double *me, int *mm, int nchans, double binsize,int nb, int *badChanMask );
static int addToMerging(double *p1, double *v1, double *e1, double *mp, double *mv,double *me, int *mm, int nchans, double binsize,int nb, int *badChanMask=NULL);
/**
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 addPointToMerging(double p1, double v1, double e1, double *mp, double *mv,double *me, int *mm, double binsize, int nb);
/**
@ -139,6 +158,8 @@ class angularConversionStatic
double* convertAngles(double pos, int nch, int *chansPerMod, angleConversionConstant **angOff, int *mF, double fo, double go, int angdir);
double convertAngle(double pos, int ich, int *chansPerMod, angleConversionConstant **angOff, int *mF, double fo, double go, int angdir);
protected:

View File

@ -7,17 +7,17 @@
#endif
#include "sls_detector_defs.h"
#include <string>
#include <fstream>
#include "angleConversionConstant.h"
//#include "angleConversionConstant.h"
//double angle(int ichan, double encoder, double totalOffset, double conv_r, double center, double offset, double tilt, int direction)
class angleConversionConstant;
using namespace std;
@ -35,7 +35,8 @@ class containing the methods to set/unset the angular conversion and merge the d
The angular conversion itself is defined by the angle() function defined in usersFunctions.cpp
*/
class angularConversionStatic : public virtual slsDetectorDefs
class angularConversionStatic
// : public virtual slsDetectorDefs
{
public:

View File

@ -45,6 +45,13 @@ class badChannelCorrections{
} \
} \
return nbad; };
static int setBadChannelCorrection(ifstream &infile, int &nbad, int *badlist, int moff){ \
int retval=readBadChannelCorrectionFile(infile,nbad,badlist); \
for (int ich=0; ich<nbad; ich++) badlist[ich]=badlist[ich]+moff; \
return retval; \
};
protected:

View File

@ -1,5 +1,9 @@
#include "postProcessing.h"
#include "postProcessingFuncs.h"
#include "angleConversionConstant.h"
#ifdef VERBOSE
#include "usersFunctions.h"
#endif
@ -12,63 +16,21 @@ postProcessing::postProcessing(){
//cout << "reg callback "<< endl;
dataReady = 0;
pCallbackArg = 0;
#ifdef VERBOSE
registerDataCallback(&defaultDataReadyFunc, NULL);
#endif
//cout << "done "<< endl;
rawDataReady = 0;
pRawDataArg = 0;
ppFun=new postProcessingFuncs();
}
int postProcessing::flatFieldCorrect(double datain, double errin, double &dataout, double &errout, double ffcoefficient, double fferr){
double 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(double datain, double errin, double &dataout, double &errout, double tau, double t){
// double data;
double 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;
};
postProcessing::~postProcessing(){delete ppFun;};
@ -82,47 +44,62 @@ void postProcessing::processFrame(int *myData, int delflag) {
string fname;
// double *fdata=NULL;
incrementProgress();
#ifdef VERBOSE
cout << "start processing"<< endl;
#endif
incrementProgress();
#ifdef VERBOSE
cout << "prog incremented"<< endl;
#endif
/** decode data */
fdata=decodeData(myData, fdata);
#ifdef VERBOSE
cout << "decode"<< endl;
#endif
fname=createFileName();
//Checking for write flag
if(*correctionMask&(1<<WRITE_FILE))
{
fname=createFileName();
//Checking for write flag
if(*correctionMask&(1<<WRITE_FILE)) {
#ifdef VERBOSE
cout << "writing raw data " << endl;
cout << "writing raw data " << endl;
#endif
//uses static function?!?!?!?
writeDataFile (fname+string(".raw"),fdata, NULL, NULL, 'i');
//uses static function?!?!?!?
writeDataFile (fname+string(".raw"),fdata, NULL, NULL, 'i');
#ifdef VERBOSE
cout << "done " << endl;
cout << "done " << endl;
#endif
}
}
if (rawDataReady)
rawDataReady(fdata,pRawDataArg);
if (*correctionMask & ~(1<<WRITE_FILE))
doProcessing(fdata,delflag, fname);
doProcessing(fdata,delflag, fname);
delete [] myData;
myData=NULL;
fdata=NULL;
delete [] myData;
delete [] fdata;
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);
if(*correctionMask&(1<<WRITE_FILE))
IncrementFileIndex();
popDataQueue(); //remove the data from the queue
queuesize=dataQueueSize();
}
@ -132,210 +109,53 @@ if(*correctionMask&(1<<WRITE_FILE))
void postProcessing::doProcessing(double *lfdata, int delflag, string fname) {
double *ang=new double[arraySize];
double *val=new double[arraySize];
double *err=new double[arraySize];
int np;
detectorData *thisData;
// /** write raw data file */
// if (*correctionMask==0 && delflag==1) {
// // delete [] fdata;
// ;
// } else {
int npos=getNumberOfPositions();
double *rcdata=NULL, *rcerr=NULL;
double *ffcdata=NULL, *ffcerr=NULL;
double *ang=NULL;
// int imod;
int np;
//string fname;
detectorData *thisData;
string ext=".dat";
// fname=createFileName();
/** rate correction */
if (*correctionMask&(1<<RATE_CORRECTION)) {
string ext=".dat";
double t=(*expTime)*1E-9;
if (GetCurrentPositionIndex()<=1) {
#ifdef VERBOSE
cout << "rate correcting " << endl;
cout << "init dataset" << endl;
#endif
rcdata=new double[getTotalNumberOfChannels()];
rcerr=new double[getTotalNumberOfChannels()];
rateCorrect(lfdata,NULL,rcdata,rcerr);
delete [] lfdata;
#ifdef VERBOSE
cout << "done " << endl;
#endif
} else {
rcdata=lfdata;
ppFun->initDataset();
}
lfdata=NULL;
#ifdef VERBOSE
cout << "add frame" << endl;
#endif
ppFun->addFrame(lfdata, &currentPosition, &currentI0, &t, (fname+ext).c_str(), NULL);
if ((GetCurrentPositionIndex()>=npos && positionFinished() && dataQueueSize()) || npos==0) {
#ifdef VERBOSE
cout << "finalize dataset" << endl;
#endif
ppFun->finalizeDataset(ang, val, err, &np);
IncrementPositionIndex();
/** flat field correction */
if (*correctionMask&(1<<FLAT_FIELD_CORRECTION)) {
#ifdef VERBOSE
cout << "ff correcting " << endl;
#endif
ffcdata=new double[getTotalNumberOfChannels()];
ffcerr=new double[getTotalNumberOfChannels()];
flatFieldCorrect(rcdata,rcerr,ffcdata,ffcerr);
delete [] rcdata;
rcdata=NULL;
if (rcerr)
delete [] rcerr;
rcerr=NULL;
#ifdef VERBOSE
cout << "done " << endl;
#endif
} else {
ffcdata=rcdata;
ffcerr=rcerr;
rcdata=NULL;
rcerr=NULL;
}
// writes angualr converted files
if (*correctionMask!=0) {
if (*correctionMask&(1<< ANGULAR_CONVERSION)) {
#ifdef VERBOSE
cout << "ang conv " << endl;
#endif
ang=convertAngles();
#ifdef VERBOSE
cout << "done - write position data file " << endl;
#endif
}
writeDataFile (fname+ext, ffcdata, ffcerr,ang);
#ifdef VERBOSE
cout << "done " << endl;
#endif
}
// if (*correctionMask&(1<< ANGULAR_CONVERSION) && getNumberOfPositions()>0) {
if (*correctionMask&(1<< ANGULAR_CONVERSION)) {
#ifdef VERBOSE
cout << "**************Current position index is " << getCurrentPositionIndex() << endl;
#endif
// if (*numberOfPositions>0) {
if (getCurrentPositionIndex()<=1) {
#ifdef VERBOSE
cout << "reset merging " << endl;
#endif
resetMerging();
}
#ifdef VERBOSE
cout << "add to merging "<< getCurrentPositionIndex() << endl;
#endif
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=finalizeMerging();
/** file writing */
incrementPositionIndex();
// cout << "unlock 1" << endl;
pthread_mutex_lock(&mp);
fname=createFileName();
pthread_mutex_unlock(&mp);
fname=createFileName();
#ifdef VERBOSE
cout << "writing merged data file" << endl;
#endif
writeDataFile (fname+ext,np,getMergedCounts(), getMergedErrors(), getMergedPositions(),'f');
#ifdef VERBOSE
cout << " done" << endl;
#endif
// if (delflag) {
// deleteMerging();
// } else {
thisData=new detectorData(getMergedCounts(),getMergedErrors(),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(*correctionMask&(1<<WRITE_FILE))
writeDataFile (fname+ext,np,ang, val, err,'f');
thisData=new detectorData(ang,val,err,getCurrentProgress(),(fname+ext).c_str(),np);
if (dataReady) {
dataReady(thisData, pCallbackArg);
delete thisData;
}
// pthread_mutex_lock(&mg);
// finalDataQueue.push(thisData);
// pthread_mutex_unlock(&mg);
// }
delete thisData;
}
//}
incrementFileIndex();
#ifdef VERBOSE
cout << "fdata is " << fdata << endl;
#endif
}
@ -422,70 +242,89 @@ void* postProcessing::processData(int delflag) {
#endif
setTotalProgress();
pthread_mutex_lock(&mp);
queuesize=dataQueue.size();
pthread_mutex_unlock(&mp);
queuesize=dataQueueSize();
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) {
while((queuesize=dataQueueSize())>0) {
/** Pop data queue */
myData=dataQueue.front(); // get the data from the queue
pthread_mutex_unlock(&mp);
#ifdef VERBOSE
cout << "data found"<< endl;
#endif
myData=dataQueueFront(); // get the data from the queue
#ifdef VERBOSE
cout << "got them"<< endl;
#endif
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);
if (checkJoinThread()) {
if (dataQueueSize()==0) {
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::checkJoinThread() {
int retval;
pthread_mutex_lock(&mp);
retval=jointhread;
pthread_mutex_unlock(&mp);
return retval;
}
void postProcessing::setJoinThread( int v) {
pthread_mutex_lock(&mp);
jointhread=v;
pthread_mutex_unlock(&mp);
}
int* postProcessing::dataQueueFront() {
int *retval=NULL;
pthread_mutex_lock(&mp);
if( !dataQueue.empty() ) {
retval=dataQueue.front();
}
pthread_mutex_unlock(&mp);
return retval;
}
int postProcessing::dataQueueSize() {
int retval;
pthread_mutex_lock(&mp);
retval=dataQueue.size();
pthread_mutex_unlock(&mp);
return retval;
}
int* postProcessing::popDataQueue() {
int *retval=NULL;
pthread_mutex_lock(&mp);
if( !dataQueue.empty() ) {
retval=dataQueue.front();
dataQueue.pop();
}
pthread_mutex_unlock(&mp);
return retval;
}
@ -502,11 +341,13 @@ detectorData* postProcessing::popFinalDataQueue() {
void postProcessing::resetDataQueue() {
int *retval=NULL;
pthread_mutex_lock(&mp);
while( !dataQueue.empty() ) {
retval=dataQueue.front();
dataQueue.pop();
delete [] retval;
}
pthread_mutex_unlock(&mp);
}
@ -523,12 +364,125 @@ void postProcessing::resetFinalDataQueue() {
void postProcessing::startThread(int delflag) {
/////////////////////////////////// Initialize dataset
//resetDataQueue();
setTotalProgress();
int nmod=getNMods();
int *chPM=new int[nmod];
int *mM=new int[nmod];
int totch=0;
#ifdef VERBOSE
cout << "init dataset stuff" << endl;
#endif
for (int im=0; im<nmod; im++) {
#ifdef VERBOSE
cout << "MODULE " << im << endl;
#endif
chPM[im]=getChansPerMod(im);
#ifdef VERBOSE
cout << "chanspermod" << endl;
#endif
mM[im]=getMoveFlag(im);
#ifdef VERBOSE
cout << "moveflag" << endl;
#endif
totch+=chPM[im];
}
#ifdef VERBOSE
cout << "total channels is " << totch << endl;
#endif
double *ffcoeff=NULL, *fferr=NULL;
if (*correctionMask&(1<<FLAT_FIELD_CORRECTION)) {
#ifdef VERBOSE
cout << "get ff " << endl;
#endif
ffcoeff=new double[totch];
fferr=new double[totch];
getFlatFieldCorrection(ffcoeff,fferr);
}
double tdead;
if (*correctionMask&(1<<RATE_CORRECTION)) {
#ifdef VERBOSE
cout << "get tau " << endl;
#endif
tdead=getRateCorrectionTau();
} else
tdead=0;
int angdir=getAngularDirection();
double to=0;
double bs=0;
double sx=0, sy=0;
double *angRad=NULL;
double *angOff=NULL;
double *angCenter=NULL;
angleConversionConstant *p=NULL;
if (*correctionMask&(1<< ANGULAR_CONVERSION)) {
#ifdef VERBOSE
cout << "get angconv " << endl;
#endif
bs=getBinSize();
to=getGlobalOffset()+getFineOffset();
angRad=new double[nmod];
angOff=new double[nmod];
angCenter=new double[nmod];
for (int im=0; im<nmod; im++) {
p=getAngularConversionPointer(im);
angRad[im]=p->r_conversion;
angOff[im]=p->offset;
angCenter[im]=p->center;
}
sx=getAngularConversionParameter(SAMPLE_X);
sy=getAngularConversionParameter(SAMPLE_Y);
}
#ifdef VERBOSE
cout << "init dataset" << endl;
#endif
ppFun->initDataset(&nmod,chPM,mM,badChannelMask, ffcoeff, fferr, &tdead, &angdir, angRad, angOff, angCenter, &to, &bs, &sx, &sy);
#ifdef VERBOSE
cout << "done" << endl;
#endif
if (*correctionMask&(1<< ANGULAR_CONVERSION)) {
arraySize=getNumberOfAngularBins();
if (arraySize<=0)
arraySize=totch;
} else {
arraySize=totch;
}
queuesize=dataQueueSize();
resetFinalDataQueue();
resetDataQueue();
/////////////////////////////////// Start thread ////////////////////////////////////////////////////////
#ifdef VERBOSE
cout << "start thread stuff" << endl;
#endif
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;
@ -553,7 +507,8 @@ void postProcessing::startThread(int delflag) {
ret = pthread_create(&dataProcessingThread, &tattr,startProcessDataNoDelete, (void*)this);
pthread_attr_destroy(&tattr);
// scheduling parameters of target thread
// scheduling parameters of target thread
ret = pthread_setschedparam(dataProcessingThread, policy, &param);
}

View File

@ -18,6 +18,9 @@
#include <queue>
#include <math.h>
class postProcessingFuncs;
using namespace std;
@ -39,7 +42,7 @@ class postProcessing : public angularConversion, public fileIO, public badChanne
public:
postProcessing();
virtual ~postProcessing(){};
virtual ~postProcessing();
@ -75,7 +78,6 @@ class postProcessing : public angularConversion, public fileIO, public badChanne
*/
virtual int setBadChannelCorrection(string fname="")=0;
static int setBadChannelCorrection(ifstream &infile, int &nbad, int *badlist, int moff){int retval=readBadChannelCorrectionFile(infile,nbad,badlist); for (int ich=0; ich<nbad; ich++) badlist[ich]=badlist[ich]+moff; return retval;};
/**
set bad channels correction
@ -85,7 +87,7 @@ class postProcessing : public angularConversion, public fileIO, public badChanne
\returns 0 if bad channel disabled, >0 otherwise
*/
virtual int setBadChannelCorrection(string fname, int &nbad, int *badlist, int off=0)=0;
using badChannelCorrections::setBadChannelCorrection;
/**
set bad channels correction
@ -96,30 +98,7 @@ class postProcessing : public angularConversion, public fileIO, public badChanne
*/
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(double datain, double errin, double &dataout, double &errout, double ffcoefficient, double 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(double datain, double errin, double &dataout, double &errout, double tau, double t);
int enableWriteToFile(int i=-1) {if (i>0) ((*correctionMask)|=(1<<WRITE_FILE)); else if (i==0) ((*correctionMask)&=~(1<< WRITE_FILE)); return (((*correctionMask)&(1<< WRITE_FILE ))>>WRITE_FILE) ;};
@ -201,6 +180,10 @@ s
*/
int* popDataQueue();
int* dataQueueFront();
int dataQueueSize();
/**
pops the data from thepostprocessed data queue
\returns pointer to the popped data or NULL if the queue is empty.
@ -209,6 +192,10 @@ s
detectorData* popFinalDataQueue();
int checkJoinThread();
void setJoinThread(int v);
/**
resets the raw data queue
\sa dataQueue
@ -234,9 +221,14 @@ s
virtual int flatFieldCorrect(double*, double*, double*, double*)=0;
virtual int getNMods()=0;
int GetCurrentPositionIndex(){ pthread_mutex_lock(&mp); int retval=getCurrentPositionIndex(); pthread_mutex_unlock(&mp); return retval;};
void IncrementPositionIndex(){ pthread_mutex_lock(&mp); incrementPositionIndex(); pthread_mutex_unlock(&mp);};
void IncrementFileIndex(){ pthread_mutex_lock(&mp); incrementFileIndex(); pthread_mutex_unlock(&mp);};
void ResetPositionIndex(){pthread_mutex_lock(&mp); resetPositionIndex(); pthread_mutex_unlock(&mp);};
void registerDataCallback(int( *userCallback)(detectorData*, void*), void *pArg) {dataReady = userCallback; pCallbackArg = pArg;};
@ -246,8 +238,14 @@ s
virtual double getRateCorrectionTau()=0;
int positionFinished(int v=-1){pthread_mutex_lock(&mp); if (v>=0) posfinished=v; int retval=posfinished; pthread_mutex_unlock(&mp); return retval;};
double getCurrentPosition() {double p; pthread_mutex_lock(&mp); p=currentPosition; pthread_mutex_unlock(&mp); return p;}
int setCurrentPosition(double v) { pthread_mutex_lock(&mp); currentPosition=v; pthread_mutex_unlock(&mp); }
protected:
@ -257,7 +255,8 @@ s
char *flatFieldDir;
char *flatFieldFile;
int64_t *expTime;
/** mutex to synchronize main and data processing threads */
pthread_mutex_t mp;
@ -309,26 +308,22 @@ s
*/
double currentI0;
int arraySize;
private:
double *fdata;
// private:
/* virtual void incrementProgress()=0; */
/* virtual double getCurrentProgress()=0; */
/* virtual void incrementFileIndex()=0; */
/* virtual int setTotalProgress()=0; */
/* virtual double* decodeData(int *datain, double *fdata=NULL)=0; */
/* virtual int getTotalNumberOfChannels()=0; */
int (*dataReady)(detectorData*,void*);
void *pCallbackArg;
int (*rawDataReady)(double*,void*);
int (*rawDataReady)(double*,void*);
void *pRawDataArg;
postProcessingFuncs *ppFun;
};

View File

@ -0,0 +1,438 @@
#include "postProcessingFuncs.h"
#include "angleConversionConstant.h"
postProcessingFuncs::postProcessingFuncs(int *nModules,int *chPerMod,int modMask[],int badCh[], double ffcoeff[], double fferr[], double* t, int *dir, double angRadius[], double angOffset[], double angCentre[], double* to, double* bs, double *sX, double *sY):
nMods(0), chansPerMod(NULL), moduleMask(NULL), badChannelMask(NULL), ffCoeff(NULL), ffErr(NULL), tDead(0), angDir(1), angConv(NULL), totalOffset(0), binSize(0), sampleX(0), sampleY(0), totalChans(0), nBins(0), mp(NULL), mv(NULL), me(NULL), mm(NULL)
{
initDataset(nModules, chPerMod,modMask,badCh, ffcoeff, fferr, t, dir, angRadius, angOffset, angCentre, to, bs, sX, sY);
}
int postProcessingFuncs::initDataset() {
if (nBins) {
mp=new double[nBins];
mv=new double[nBins];
me=new double[nBins];
mm=new int[nBins];
resetMerging(mp,mv,me,mm, nBins);
} else {
mv=new double[totalChans];
me=new double[totalChans];
}
totalI0=0;
return 0;
}
int postProcessingFuncs::finalizeDataset(double *ang, double *val, double *err, int *np) {
if (nBins)
*np=finalizeMerging(mp,mv,me,mm,nBins);
else
*np=totalChans;
if (totalI0<=0)
totalI0=1.;
for (int ip=0; ip<(*np); ip++) {
if (mp)
ang[ip]=mp[ip];
if (mv)
val[ip]=mv[ip]*totalI0;
if (me)
err[ip]=me[ip]*totalI0;
}
if (mp)
delete [] mp;
if (mv)
delete [] mv;
if (me)
delete [] me;
if (mm)
delete [] mm;
return 0;
}
int postProcessingFuncs::addFrame(double *data, double *pos, double *I0, double *expTime, const char *filename, int *var) {
double p1, vin, ein, vout, eout;
double e=0.;
double i0=*I0;
if (i0>0)
totalI0+=i0;
for (int ich=0; ich<totalChans; ich++) {
vin=data[ich];
ein=0;
vout=data[ich];
if (vout>=0)
eout=sqrt(vout);
else
eout=0;
//ratecorrect
if (tDead) {
rateCorrect(vin, ein, vout, eout, tDead, *expTime);
vin=vout;
ein=eout;
}
//ffcorrect
if (ffCoeff) {
if (ffErr)
e=ffErr[ich];
else
e=0;
flatFieldCorrect(vin, ein, vout, eout, ffCoeff[ich], e);
}
//i0correct
if (i0>0) {
vout/=i0;
eout/=i0;
}
if (badChannelMask)
if (badChannelMask[ich])
continue;
if (nBins) {
//angconv
//check module mask?!?!?!?
p1=convertAngle(*pos,ich,chansPerMod,angConv,moduleMask,totalOffset,0,angDir);
addPointToMerging(p1,vout,eout,mp,mv,me,mm, binSize, nBins);
} else {
mv[ich]+=vout;
me[ich]+=eout*eout;
}
}
}
int postProcessingFuncs::initDataset(int *nModules,int *chPerMod,int modMask[],int badCh[], double ffcoeff[], double fferr[], double* t, int *dir, double angRadius[], double angOffset[], double angCenter[], double* to, double* bs, double *sX, double *sY) {
#ifdef VERBOSE
cout << "delete pointers " << endl;
#endif
deletePointers();
#ifdef VERBOSE
cout << "nmod " << endl;
#endif
if (nModules)
nMods=*nModules;
else
nMods=0;
#ifdef VERBOSE
cout << "tdead " << endl;
#endif
if (t)
tDead=*t;
else
t=0;
#ifdef VERBOSE
cout << "toffset " << endl;
#endif
if (to)
totalOffset=*to;
else
to=0;
#ifdef VERBOSE
cout << "binsize " << endl;
#endif
if (bs)
binSize=*bs;
else
binSize=0;
#ifdef VERBOSE
cout << "samplex " << endl;
#endif
if (sX)
sampleX=*sX;
else
sampleX=0;
#ifdef VERBOSE
cout << "sampley " << endl;
#endif
if (sY)
sampleY=*sY;
else
sampleY=0;
#ifdef VERBOSE
cout << "angdir " << endl;
#endif
if (dir)
angDir=*dir;
else
angDir=1;
totalChans=0;
chansPerMod=new int [nMods];
moduleMask=new int [nMods];
nBins=0;
if (angRadius && angOffset && angCenter && (binSize>0)) {
angConv=new angleConversionConstant*[nMods];
nBins=360./binSize+1;
}
for (int im=0; im<nMods; im++) {
#ifdef VERBOSE
cout << "MODULE "<< im << endl;
#endif
chansPerMod[im]=chPerMod[im];
moduleMask[im]=modMask[im];
totalChans+=chansPerMod[im];
if (angConv) {
angConv[im]=new angleConversionConstant(angCenter[im], angRadius[im], angOffset[im], 0);
}
}
#ifdef VERBOSE
cout << "badchans " << endl;
#endif
if (badCh)
badChannelMask= new int[totalChans];
#ifdef VERBOSE
cout << "ffcoeff " << endl;
#endif
if (ffcoeff)
ffCoeff=new double[totalChans];
#ifdef VERBOSE
cout << "fferr " << endl;
#endif
if (fferr)
ffErr=new double[totalChans];
for (int ich=0; ich<totalChans; ich++) {
if (badChannelMask)
badChannelMask[ich]=badCh[ich];
if (ffCoeff)
ffCoeff[ich]=ffcoeff[ich];
if (ffErr)
ffErr[ich]=fferr[ich];
}
return 0;
}
void postProcessingFuncs::deletePointers() {
delete [] chansPerMod;
delete [] moduleMask;
if (badChannelMask)
delete [] badChannelMask;
if (ffCoeff)
delete [] ffCoeff;
if (ffErr)
delete [] ffErr;
if (angConv) {
for (int im=0; im<nMods; im++) {
if (angConv[im])
delete angConv[im];
}
delete [] angConv;
}
}
postProcessingFuncs::~postProcessingFuncs(){
deletePointers();
}
int postProcessingFuncs::flatFieldCorrect(double datain, double errin, double &dataout, double &errout, double ffcoefficient, double fferr){
double 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 postProcessingFuncs::rateCorrect(double datain, double errin, double &dataout, double &errout, double tau, double t){
// double data;
double 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 postProcessingFuncs::calculateFlatField(int* nModules, int *chPerMod, int *moduleMask, int *badChannelMask, double *ffData, double *ffCoeff, double *ffErr) {
int nmed=0, im=0;
double *xmed;
if (chPerMod==NULL)
return -1;
// if (moduleMask==NULL)
// return -1;
if (ffData==NULL)
return -1;
if (ffErr==NULL)
return -1;
int totch=0;
for (int im=0; im<*nModules; im++) {
totch+=chPerMod[im];
}
xmed=new double[totch];
for (int ich=0; ich<totch; ich++) {
if (badChannelMask)
if (badChannelMask[ich])
continue;
// calculate median if ffData is positive and channel is good
if (ffData[ich]>0) {
im=0;
while ((im<nmed) && (xmed[im]<ffData[ich]))
im++;
for (int i=nmed; i>im; i--)
xmed[i]=xmed[i-1];
xmed[im]=ffData[ich];
nmed++;
}
}
if (nmed>1 && xmed[nmed/2]>0) {
#ifdef VERBOSE
std::cout<< "Flat field median is " << xmed[nmed/2] << " calculated using "<< nmed << " points" << std::endl;
#endif
for (int ich=0; ich<totch; ich++) {
if (badChannelMask)
if (badChannelMask[ich]) {
ffCoeff[ich]=0.;
ffErr[ich]=1.;
}
if (ffData[ich]>0) {
ffCoeff[ich]=xmed[nmed/2]/ffData[ich];
ffErr[ich]=ffCoeff[ich]*sqrt(ffData[ich])/ffData[ich];
} else {
ffCoeff[ich]=0.;
ffErr[ich]=1.;
}
}
}
delete [] xmed;
return 0;
}

View File

@ -0,0 +1,80 @@
#ifndef POSTPROCESSINGFUNCS_H
#define POSTPROCESSINGFUNC_H
#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cstring>
#include <string>
#include <sstream>
#include <queue>
#include <math.h>
#include "angularConversionStatic.h"
class angleConversionConstant;
using namespace std;
class postProcessingFuncs : public virtual angularConversionStatic
{
public:
postProcessingFuncs(int *nModules=NULL,int *chPerMod=NULL,int *modMask=NULL,int *badChMask=NULL, double *ffcoeff=NULL, double *fferr=NULL, double* t=NULL, int *dir=NULL, double *angRadius=NULL, double *angOffset=NULL, double *angCentre=NULL, double* to=NULL, double* bs=NULL, double *sX=NULL, double *sY=NULL);
~postProcessingFuncs();
int initDataset(int *nModules,int *chPerMod,int modMask[],int badCh[], double ffcoeff[], double fferr[], double* tDead, int *dir, double angRadius[], double angOffset[], double angCentre[], double* to, double* bs, double *sX, double *sY);
int initDataset();
int finalizeDataset(double ang[], double val[], double err[], int *np);
int addFrame(double data[], double *pos, double *IO, double *expTime, const char *filename, int *var=0);
static int calculateFlatField(int* nModules, int *chPerMod, int moduleMask[], int badChannelMask[], double ffData[], double ffCoeff[], double ffErr[]);
static int flatFieldCorrect(double datain, double errin, double &dataout, double &errout, double ffcoefficient, double fferr);
static int rateCorrect(double datain, double errin, double &dataout, double &errout, double tau, double t);
private:
void deletePointers();
int nMods;
int *chansPerMod;
int *moduleMask;
int *badChannelMask;
double *ffCoeff;
double *ffErr;
double tDead;
int angDir;
angleConversionConstant **angConv;
double totalOffset;
double binSize;
double sampleX;
double sampleY;
int totalChans;
int nBins;
double totalI0;
double *mp, *mv,*me;
int *mm;
};
#endif