583 lines
14 KiB
C
583 lines
14 KiB
C
/*------------------------------------------------------------------------
|
|
F I T C E N T E R
|
|
|
|
A simple peak finding and center of gravity determination module for SICS.
|
|
|
|
copyright: see copyright.h
|
|
|
|
Mark Koennecke, October 1997
|
|
|
|
Added center of edge finding
|
|
|
|
Mark Koennecke, Octover 2011
|
|
|
|
-----------------------------------------------------------------------------*/
|
|
#include <limits.h>
|
|
#include <math.h>
|
|
#include "sics.h"
|
|
#include "counter.h"
|
|
#include "motor.h"
|
|
#include "scan.h"
|
|
#include "fitcenter.h"
|
|
#define THRESHOLD .1
|
|
/*--------------------------------------------------------------------------*/
|
|
typedef struct __FitCenter {
|
|
pObjectDescriptor pDes;
|
|
pScanData pScan;
|
|
long *lCounts;
|
|
float *fAxis;
|
|
char pName[132];
|
|
float fCenter;
|
|
float fStddev;
|
|
long lPeak;
|
|
float FWHM;
|
|
int iNP;
|
|
} FitCenter;
|
|
/*-------------------------------------------------------------------------*/
|
|
pFit CreateFitCenter(pScanData pScan)
|
|
{
|
|
pFit pNew = NULL;
|
|
|
|
pNew = (pFit) malloc(sizeof(FitCenter));
|
|
if (!pNew) {
|
|
return NULL;
|
|
}
|
|
memset(pNew, 0, sizeof(FitCenter));
|
|
|
|
pNew->pDes = CreateDescriptor("FitCenter");
|
|
if (!pNew->pDes) {
|
|
free(pNew);
|
|
return NULL;
|
|
}
|
|
pNew->pScan = pScan;
|
|
return pNew;
|
|
}
|
|
|
|
/*-------------------------------------------------------------------------*/
|
|
void DeleteFitCenter(void *pData)
|
|
{
|
|
pFit self = NULL;
|
|
|
|
self = (pFit) pData;
|
|
assert(self);
|
|
|
|
if (self->pDes) {
|
|
DeleteDescriptor(self->pDes);
|
|
}
|
|
if (self->lCounts) {
|
|
free(self->lCounts);
|
|
}
|
|
if (self->fAxis) {
|
|
free(self->fAxis);
|
|
}
|
|
free(self);
|
|
}
|
|
|
|
/*--------------------------------------------------------------------------*/
|
|
static int Init(pFit self)
|
|
{
|
|
pMotor pMot = NULL;
|
|
float fZero, fSign;
|
|
int i;
|
|
|
|
/* collect data first */
|
|
self->iNP = GetScanNP(self->pScan);
|
|
if (self->lCounts) {
|
|
free(self->lCounts);
|
|
self->lCounts = NULL;
|
|
}
|
|
if (self->fAxis) {
|
|
free(self->fAxis);
|
|
self->fAxis = NULL;
|
|
}
|
|
self->lCounts = (long *) malloc(self->iNP * sizeof(long));
|
|
self->fAxis = (float *) malloc(self->iNP * sizeof(float));
|
|
if ((!self->lCounts) || (!self->fAxis)) {
|
|
return 0;
|
|
}
|
|
|
|
GetScanCounts(self->pScan, self->lCounts, self->iNP);
|
|
GetScanVar(self->pScan, 0, self->fAxis, self->iNP);
|
|
GetScanVarName(self->pScan, 0, self->pName, 131);
|
|
|
|
return 1;
|
|
}
|
|
|
|
/*-------------------------------------------------------------------------*/
|
|
static int FindMax(pFit self)
|
|
{
|
|
long lMax = LONG_MIN;
|
|
int i, iMax;
|
|
|
|
for (i = 0, iMax = 0; i < self->iNP; i++) {
|
|
if (self->lCounts[i] > lMax) {
|
|
lMax = self->lCounts[i];
|
|
iMax = i;
|
|
}
|
|
}
|
|
|
|
return iMax;
|
|
}
|
|
|
|
/*-------------------------------------------------------------------------*/
|
|
static int CalculateFitIntern(pFit self)
|
|
{
|
|
int i, iLeft, iRight, iRet, iPeak, iLast;
|
|
long lHalf;
|
|
float fFrac, fDelta, fStep, fLeft, fRight, fNenner, fTest;
|
|
float fLFrac, fRFrac, fSum, fCI;
|
|
|
|
iLeft = 0;
|
|
|
|
/* find the maximum counts */
|
|
iRet = FindMax(self);
|
|
self->lPeak = self->lCounts[iRet];
|
|
/* a default fit is the peak maximum. This
|
|
helps optimise to do the right thing if no proper calculation
|
|
could be performed
|
|
*/
|
|
self->fCenter = self->fAxis[iRet];
|
|
iPeak = iRet;
|
|
|
|
if (self->lPeak < 3) {
|
|
return -3;
|
|
}
|
|
|
|
/* find FWHM, first from Right then from left */
|
|
iRet = 1;
|
|
lHalf = self->lPeak / 2;
|
|
iLeft = iPeak;
|
|
while ((self->lCounts[iLeft] > lHalf) && (iLeft > 0)) {
|
|
iLeft--;
|
|
}
|
|
if (iLeft == 0) {
|
|
iRet = -1;
|
|
iLeft++;
|
|
}
|
|
fNenner =
|
|
(float) self->lCounts[iLeft + 1] - (float) self->lCounts[iLeft];
|
|
if (fNenner < 0.) {
|
|
fTest = -fNenner;
|
|
} else {
|
|
fTest = fNenner;
|
|
}
|
|
if (fTest < .00001) {
|
|
return -4;
|
|
}
|
|
fFrac = ((float) self->lCounts[iLeft + 1] - (float) lHalf) / fNenner;
|
|
fLeft = self->fAxis[iLeft] +
|
|
fFrac * (self->fAxis[iLeft + 1] - self->fAxis[iLeft]);
|
|
fLFrac = fFrac;
|
|
|
|
iRight = iPeak;
|
|
while ((self->lCounts[iRight] > lHalf) && (iRight < self->iNP)) {
|
|
iRight++;
|
|
}
|
|
if (iRight >= self->iNP - 1) {
|
|
iRet = -2;
|
|
iRight = self->iNP - 1;
|
|
}
|
|
fNenner =
|
|
(float) self->lCounts[iRight - 1] - (float) self->lCounts[iRight];
|
|
if (fNenner < 0.) {
|
|
fTest = -fNenner;
|
|
} else {
|
|
fTest = fNenner;
|
|
}
|
|
if (fTest < .00001) {
|
|
return -4;
|
|
}
|
|
fFrac = ((float) self->lCounts[iRight - 1] - (float) lHalf) / fNenner;
|
|
fRight = self->fAxis[iRight] - fFrac * (self->fAxis[1] - self->fAxis[0]);
|
|
fRFrac = fFrac;
|
|
|
|
if (iRight != iLeft) {
|
|
self->FWHM = fRight - fLeft;
|
|
} else {
|
|
self->FWHM = 0.07;
|
|
}
|
|
/*
|
|
self->fCenter = fLeft + (fRight - fLeft)/2.;
|
|
*/
|
|
fSum = lHalf * (fLFrac);
|
|
fNenner = lHalf;
|
|
/* center of gravity calculation between iLeft and iRight */
|
|
for (i = iLeft + 1; i < iRight; i++) {
|
|
fSum += self->lCounts[i] * (fLFrac + i - iLeft);
|
|
fNenner += self->lCounts[i];
|
|
}
|
|
fSum += (iRight - 1 + fRFrac + fLFrac - iLeft) * lHalf;
|
|
fNenner += lHalf;
|
|
fCI = fSum / fNenner;
|
|
/* iLeft is zero, zero point shift by lFrac */
|
|
fCI -= fLFrac;
|
|
fTest = fNenner;
|
|
fNenner = modff(fCI, &fSum);
|
|
self->fCenter = self->fAxis[(int) (fSum) + iLeft] +
|
|
fNenner * (self->fAxis[iLeft + 1] - self->fAxis[iLeft]);
|
|
|
|
/* error calculation */
|
|
fTest = 2 * (sqrtf(fTest)) / fTest;
|
|
/* because we are really summing from iLeft */
|
|
fTest = fCI * fTest;
|
|
if (fTest < 0)
|
|
fTest = -fTest;
|
|
self->fStddev = fTest * (self->fAxis[iLeft + 1] - self->fAxis[iLeft]);
|
|
|
|
return iRet;
|
|
}
|
|
/*-------------------------------------------------------------------------*/
|
|
static int CalculateEdgeIntern(pFit self)
|
|
{
|
|
int i, iRet, iPeak, nSlope = 0;;
|
|
long lCount;
|
|
float fNenner,fSum, fCI;
|
|
|
|
|
|
/* find the maximum counts */
|
|
iRet = FindMax(self);
|
|
self->lPeak = self->lCounts[iRet];
|
|
/* a default fit is the peak maximum. This
|
|
helps optimise to do the right thing if no proper calculation
|
|
could be performed
|
|
*/
|
|
if (self->lPeak < 3) {
|
|
return -3;
|
|
}
|
|
|
|
/*
|
|
* calculate the COG between .9*max < counts < .1*max
|
|
*/
|
|
fSum = 0.;
|
|
fNenner = 0.;
|
|
for (i = 0; i < self->iNP; i++) {
|
|
lCount = self->lCounts[i];
|
|
if(lCount > self->lPeak *.20 && lCount < self->lPeak * .80){
|
|
fSum += lCount * self->fAxis[i];
|
|
fNenner += lCount;
|
|
nSlope++;
|
|
}
|
|
}
|
|
if(fNenner > .0001){
|
|
fCI = fSum / fNenner;
|
|
} else {
|
|
return -3;
|
|
}
|
|
self->fCenter = fCI;
|
|
self->fStddev = 1.0;
|
|
|
|
if(nSlope < 3) {
|
|
iRet = -7;
|
|
}
|
|
return iRet;
|
|
}
|
|
|
|
/*--------------------------------------------------------------------------*/
|
|
int CalculateFit(pFit self)
|
|
{
|
|
int iRet;
|
|
|
|
|
|
/* initialise */
|
|
iRet = Init(self);
|
|
if (!iRet) {
|
|
return iRet;
|
|
}
|
|
return CalculateFitIntern(self);
|
|
}
|
|
|
|
/*------------------------------------------------------------------------*/
|
|
int CalculateFitFromData(pFit self, float *fAxis, long *lCounts,
|
|
int iLength)
|
|
{
|
|
int iRet;
|
|
float *fOld;
|
|
long *lOld;
|
|
int iOld;
|
|
|
|
/* save old data. The normal handling is to have a local copy of
|
|
data with fit. This is handled differently in this function, where
|
|
it is assumed, that the data is owned by the caller. In order not
|
|
to upset fit's scheme for data handling data pointers are swapped.
|
|
This is more efficient then copying the data around without real
|
|
need.
|
|
*/
|
|
iOld = self->iNP;
|
|
fOld = self->fAxis;
|
|
lOld = self->lCounts;
|
|
|
|
/* put in new */
|
|
self->iNP = iLength;
|
|
self->fAxis = fAxis;
|
|
self->lCounts = lCounts;
|
|
|
|
iRet = CalculateFitIntern(self);
|
|
|
|
/* reset data pointers */
|
|
self->iNP = iOld;
|
|
self->fAxis = fOld;
|
|
self->lCounts = lOld;
|
|
|
|
return iRet;
|
|
}
|
|
|
|
/*------------------------------------------------------------------------*/
|
|
void GetFitResults(pFit self, float *fCenter, float *fStdDev,
|
|
float *FWHM, float *fMax)
|
|
{
|
|
assert(self);
|
|
*fCenter = self->fCenter;
|
|
*FWHM = self->FWHM;
|
|
*fMax = (float) self->lPeak;
|
|
*fStdDev = self->fStddev;
|
|
}
|
|
|
|
/*--------------------------------------------------------------------------*/
|
|
int DriveCenter(pFit self, SConnection * pCon, SicsInterp * pSics)
|
|
{
|
|
int iRet;
|
|
|
|
assert(self);
|
|
assert(pCon);
|
|
assert(pSics);
|
|
|
|
iRet =
|
|
StartMotor(pServ->pExecutor, pSics, pCon, self->pName,
|
|
RUNDRIVE, self->fCenter);
|
|
if (!iRet) {
|
|
return 0;
|
|
}
|
|
|
|
iRet = Wait4Success(pServ->pExecutor);
|
|
switch (iRet) {
|
|
case DEVINT:
|
|
if (SCGetInterrupt(pCon) == eAbortOperation) {
|
|
SCSetInterrupt(pCon, eContinue);
|
|
SCWrite(pCon, "Driving Aborted", eError);
|
|
}
|
|
return 0;
|
|
break;
|
|
case DEVDONE:
|
|
SCWrite(pCon, "Driving to center done", eValue);
|
|
break;
|
|
default:
|
|
SCWrite(pCon,
|
|
"WARNING: driving to center finished with problems", eWarning);
|
|
break;
|
|
|
|
}
|
|
return 1;
|
|
}
|
|
|
|
/*--------------------------------------------------------------------------*/
|
|
int FitFactory(SConnection * pCon, SicsInterp * pSics, void *pData,
|
|
int argc, char *argv[])
|
|
{
|
|
pFit self = NULL;
|
|
pScanData pScan = NULL;
|
|
pDummy pDum = NULL;
|
|
CommandList *pCom = NULL;
|
|
char pBueffel[132];
|
|
int iRet, iRet1, iRet2;
|
|
|
|
if (argc < 2) {
|
|
SCWrite(pCon,
|
|
"ERROR: MakeFit expected the name of a scan command as a parameter",
|
|
eError);
|
|
return 0;
|
|
}
|
|
|
|
pCom = FindCommand(pSics, argv[1]);
|
|
if (!pCom) {
|
|
snprintf(pBueffel,sizeof(pBueffel), "ERROR: scan command %s not found", argv[1]);
|
|
SCWrite(pCon, pBueffel, eError);
|
|
return 0;
|
|
}
|
|
pDum = (pDummy) pCom->pData;
|
|
if (!pDum) {
|
|
snprintf(pBueffel,sizeof(pBueffel)-1, "ERROR: scan command %s not found", argv[1]);
|
|
SCWrite(pCon, pBueffel, eError);
|
|
return 0;
|
|
}
|
|
if (strcmp(pDum->pDescriptor->name, "ScanObject") != 0) {
|
|
snprintf(pBueffel,sizeof(pBueffel)-1, "ERROR: %s is no scan command", argv[1]);
|
|
SCWrite(pCon, pBueffel, eError);
|
|
return 0;
|
|
}
|
|
|
|
pScan = (pScanData) pDum;
|
|
self = CreateFitCenter(pScan);
|
|
if (!self) {
|
|
SCWrite(pCon, "ERROR: Failure to allocate Peak & Center command",
|
|
eError);
|
|
return 0;
|
|
}
|
|
|
|
iRet = AddCommand(pSics, "peak", FitWrapper, DeleteFitCenter, self);
|
|
iRet1 = AddCommand(pSics, "center", CenterWrapper, NULL, self);
|
|
iRet2 = AddCommand(pSics, "edge", EdgeWrapper, NULL, self);
|
|
if (!iRet || !iRet1) {
|
|
snprintf(pBueffel,sizeof(pBueffel)-1,
|
|
"ERROR: duplicate commands peak and center not created");
|
|
SCWrite(pCon, pBueffel, eError);
|
|
DeleteFitCenter((void *) self);
|
|
return 0;
|
|
}
|
|
return 1;
|
|
}
|
|
/*--------------------------------------------------------------------------*/
|
|
int FitWrapper(SConnection * pCon, SicsInterp * pSics, void *pData,
|
|
int argc, char *argv[])
|
|
{
|
|
pFit self = NULL;
|
|
int iRet;
|
|
char pBueffel[256];
|
|
pDynString buf = NULL;
|
|
|
|
self = (pFit) pData;
|
|
assert(self);
|
|
|
|
iRet = CalculateFit(self);
|
|
switch (iRet) {
|
|
case 0:
|
|
SCWrite(pCon, "ERROR: failure to fit your data!", eError);
|
|
return 0;
|
|
break;
|
|
case -1:
|
|
SCWrite(pCon, "WARNING: could not find left hand half width",
|
|
eWarning);
|
|
SCWrite(pCon,
|
|
"Fit Results most certainly dodgy, SICS suggests measuring a full peak",
|
|
eWarning);
|
|
break;
|
|
case -2:
|
|
SCWrite(pCon, "WARNING: could not find right hand half width", eError);
|
|
SCWrite(pCon,
|
|
"Fit Results most certainly dodgy, SICS suggests measuring a full peak",
|
|
eWarning);
|
|
break;
|
|
case -3:
|
|
SCWrite(pCon, "ERROR: No counts found in Fit!", eError);
|
|
return 0;
|
|
break;
|
|
case -4:
|
|
SCWrite(pCon, "ERROR: Insufficient counts in peak", eError);
|
|
return 0;
|
|
break;
|
|
}
|
|
|
|
/*
|
|
This is a little feature to get the peak without rubbish for
|
|
the TAS routines
|
|
*/
|
|
if (argc > 1) {
|
|
strtolower(argv[1]);
|
|
if (strcmp(argv[1], "value") == 0) {
|
|
snprintf(pBueffel,sizeof(pBueffel)-1, "%f", self->fCenter);
|
|
SCWrite(pCon, pBueffel, eValue);
|
|
return 1;
|
|
}
|
|
if (strcmp(argv[1], "data") == 0) {
|
|
snprintf(pBueffel, 255, "%f,%f,%ld",
|
|
self->fCenter, self->FWHM, self->lPeak);
|
|
SCWrite(pCon, pBueffel, eValue);
|
|
return 1;
|
|
}
|
|
/* print peak position and its uncertainty */
|
|
if (strcmp(argv[1], "valerr") == 0) {
|
|
snprintf(pBueffel, sizeof(pBueffel)-1, "%f,%f", self->fCenter, self->fStddev);
|
|
SCWrite(pCon, pBueffel, eValue);
|
|
return 1;
|
|
}
|
|
}
|
|
|
|
/* print results */
|
|
buf = CreateDynString(256,256);
|
|
if(buf != NULL){
|
|
snprintf(pBueffel,sizeof(pBueffel)-1, "Estimated Peak Center: %f, StdDev: %f \n",
|
|
self->fCenter, self->fStddev);
|
|
DynStringConcat(buf,pBueffel);
|
|
DynStringConcatChar(buf,'\n');
|
|
snprintf(pBueffel,sizeof(pBueffel)-1, "Maximum peak height: %ld\n", self->lPeak);
|
|
DynStringConcat(buf,pBueffel);
|
|
DynStringConcatChar(buf,'\n');
|
|
snprintf(pBueffel,sizeof(pBueffel)-1, "Approximate FWHM: %f\n", self->FWHM);
|
|
DynStringConcat(buf,pBueffel);
|
|
DynStringConcatChar(buf,'\n');
|
|
SCWrite(pCon, GetCharArray(buf), eValue);
|
|
DeleteDynString(buf);
|
|
} else {
|
|
SCWrite(pCon,"ERROR: out of memory printing fitcenter results",eError);
|
|
return 0;
|
|
}
|
|
|
|
return 1;
|
|
}
|
|
/*---------------------------------------------------------------------------*/
|
|
int EdgeWrapper(SConnection * pCon, SicsInterp * pSics, void *pData,
|
|
int argc, char *argv[])
|
|
{
|
|
pFit self = NULL;
|
|
int iRet;
|
|
char pBueffel[256];
|
|
pDynString buf = NULL;
|
|
|
|
self = (pFit) pData;
|
|
assert(self);
|
|
|
|
Init(self);
|
|
iRet = CalculateEdgeIntern(self);
|
|
switch (iRet) {
|
|
case 0:
|
|
SCWrite(pCon, "ERROR: failure to fit your data!", eError);
|
|
return 0;
|
|
break;
|
|
case -3:
|
|
SCWrite(pCon, "ERROR: No counts found in Fit!", eError);
|
|
return 0;
|
|
break;
|
|
case -4:
|
|
SCWrite(pCon, "ERROR: Insufficient counts in peak", eError);
|
|
return 0;
|
|
break;
|
|
case -7:
|
|
SCWrite(pCon,
|
|
"WARNING: not enough points in slope, results may be inreliable, remeasure with smaller step width",
|
|
eWarning);
|
|
break;
|
|
}
|
|
|
|
/*
|
|
This is a little feature to get the peak without rubbish for
|
|
the TAS routines
|
|
*/
|
|
if (argc > 1) {
|
|
strtolower(argv[1]);
|
|
if (strcmp(argv[1], "value") == 0) {
|
|
snprintf(pBueffel,sizeof(pBueffel)-1, "%f", self->fCenter);
|
|
SCWrite(pCon, pBueffel, eValue);
|
|
return 1;
|
|
}
|
|
}
|
|
|
|
/* print results */
|
|
snprintf(pBueffel,sizeof(pBueffel)-1, "Estimated Edge Center: %f\n",
|
|
self->fCenter);
|
|
SCWrite(pCon, pBueffel, eValue);
|
|
|
|
return 1;
|
|
}
|
|
/*---------------------------------------------------------------------------*/
|
|
int CenterWrapper(SConnection * pCon, SicsInterp * pSics, void *pData,
|
|
int argc, char *argv[])
|
|
{
|
|
pFit self = NULL;
|
|
int iRet;
|
|
|
|
self = (pFit) pData;
|
|
assert(self);
|
|
|
|
CalculateFit(self);
|
|
return DriveCenter(self, pCon, pSics);
|
|
}
|