From b3e425b965ca7a30cf3f13824003e8204270a49e Mon Sep 17 00:00:00 2001 From: koennecke Date: Thu, 13 Aug 2009 07:29:25 +0000 Subject: [PATCH] - Added an autoamtic window calculation routine for RITA --- amorset.c | 2 +- autowin.c | 178 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ autowin.h | 14 +++++ lmd200.c | 47 +++++++++----- make_gen | 2 +- psi.c | 3 + tdchm.c | 1 + 7 files changed, 229 insertions(+), 18 deletions(-) create mode 100644 autowin.c create mode 100644 autowin.h diff --git a/amorset.c b/amorset.c index 32a9478..1ae3fa3 100644 --- a/amorset.c +++ b/amorset.c @@ -7,7 +7,7 @@ Mark Koennecke, October 2005 - Commented supprot for slit 5 away as this is gone now + Commented support for slit 5 away as this is gone now Mark Koennecke, March 2009 --------------------------------------------------------------------*/ diff --git a/autowin.c b/autowin.c new file mode 100644 index 0000000..f60ced8 --- /dev/null +++ b/autowin.c @@ -0,0 +1,178 @@ +/** + * This is a routine which calculates automatic windows for RITA's + * monochromatic window mode. In this mode data from 9 analyzer blades + * is displayed on the detector. RITA needs the windows in which to + * sum the contributions of each detector blade. The algorithm to find + * them goes like this: + * 1) the HM is summed in y, resulting in a 1D histogram + * 2) In the 1D histogram peaks are located. + * 3) For each peak from 2, the COG is calculated to get a good peak + * center + * 4) The window is the peak center +- 9 bixels. + * + * copyright: see file COPYRIGHT + * + * Mark Koennecke, July 2009 + */ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +extern float nintf(float x); +/*------------------------------------------------------------------------*/ +static HistInt *sumY(pHistMem hm, SConnection *pCon) +{ + int dim[10], i; + char command[30]; + HistInt *data = NULL; + + GetHistDim(hm,dim,&i); + snprintf(command,30,"sum:0:0:%d", dim[1]); + data = hm->pDriv->SubSample(hm->pDriv, pCon, 0, command); + /* + for(i = 0; i < dim[0]; i++){ + printf(" %d %d\n", i, data[i] ); + } + */ + return data; +} +/*-----------------------------------------------------------------------*/ +static int findPeaks(HistInt *data, int length){ + int result, i, j; + float sum, avg; + + result = LLDcreate(sizeof(int)); + + for(i = 3; i < length - 3; i++){ + sum = 0.; + for(j = i -3; j < i +3; j++){ + if(j != 4){ + sum += data[j]; + } + } + avg = sum/6.; + if(data[i] > 8.*sqrt(avg)){ + if(data[i+1] < data[i] && data[i] > data[i-1]){ + LLDnodeAppendFrom(result,&i); + } + } + } + + return result; +} +/*------------------------------------------------------------------------*/ +static float calcCOG(HistInt *data, int length, int startPos) +{ + float threshold, sum = .0, xsum = .0, cog; + int i; + + threshold = data[startPos] * .33; + + /* sum to left */ + i = startPos; + while(i > 0 && data[i] > threshold){ + sum += data[i]; + xsum += data[i]*i; + i--; + } + + /* sum to right */ + i = startPos+1; + while(i < length && data[i] > threshold){ + sum += data[i]; + xsum += data[i]*i; + i++; + } + + cog = xsum/sum; + return cog; +} +/*-----------------------------------------------------------------------*/ +static void printPeaks(SConnection *pCon, int peakList, + HistInt *data, int length) +{ + int status, startPos, oldStartPos = 0; + float cog; + pDynString buf; + + SCStartBuffering(pCon); + status = LLDnodePtr2First(peakList); + while(status > 0){ + LLDnodeDataTo(peakList,&startPos); + cog = calcCOG(data, length, startPos); + startPos = (int)nintf(cog); + if((startPos - oldStartPos) > 2){ + SCPrintf(pCon,eValue,"%d,%d,%d", startPos - 5, startPos + 5, + startPos); + oldStartPos = startPos; + } + status = LLDnodePtr2Next(peakList); + } + buf = SCEndBuffering(pCon); + SCWrite(pCon,GetCharArray(buf), eValue); +} +/*------------------------------------------------------------------------*/ +static int CalcCmd(pSICSOBJ self, SConnection *pCon, pHdb commandNode, + pHdb par[], int nPar) +{ + + HistInt *data = NULL; + pHistMem hm = NULL; + int peakList, dim[10]; + char hmName[30]; + + if(nPar < 1){ + strcpy(hmName,"hm"); + } else { + strncpy(hmName,par[0]->value.v.text,30); + } + + hm = FindCommandData(pServ->pSics, + hmName, "HistMem"); + if(hm == NULL || hm->pDriv->SubSample == NULL){ + SCPrintf(pCon,eError,"ERROR: hm %s not found", hmName); + return 0; + } + + data = sumY(hm,pCon); + if(data == NULL){ + SCWrite(pCon,"ERROR: subsampling failed", eError); + return 0; + } + + GetHistDim(hm,dim,&peakList); + + peakList = findPeaks(data, dim[0]); + + printPeaks(pCon,peakList, data, dim[0]); + free(data); + + return 1; +} +/*----------------------------------------------------------------*/ +int MakeRitaWin(SConnection *pCon, SicsInterp *pSics, void *pData, + int argc, char *argv[]) +{ + pSICSOBJ pNew = NULL; + pHdb cmd, node; + + pNew = SetupSICSOBJ(pCon,pSics,pData,argc,argv); + if(pNew == NULL){ + return 0; + } + + cmd = AddSICSHdbPar(pNew->objectNode, "calc", usSpy, + MakeSICSFunc(CalcCmd)); + SetHdbProperty(cmd,"type","command"); + SetHdbProperty(cmd,"priv","spy"); + node = MakeSICSHdbPar("hm",usSpy,MakeHdbText("banana")); + AddHipadabaChild(cmd,node,NULL); + + return 1; +} diff --git a/autowin.h b/autowin.h new file mode 100644 index 0000000..d4cd611 --- /dev/null +++ b/autowin.h @@ -0,0 +1,14 @@ +/** + * This is the calculation of automatic windows for RITA's + * monochromatic imaging mode. For more see comment in autowin.c. + * + * copyright: see file COPYRIGHT + * + * Mark Koennecke, July 2009 + */ +#ifndef AUTOWIN_H_ +#define AUTOWIN_H_ +int MakeRitaWin(SConnection *pCon, SicsInterp *pSics, void *pData, + int argc, char *argv[]); + +#endif /*AUTOWIN_H_*/ diff --git a/lmd200.c b/lmd200.c index a7356fe..2c9279b 100644 --- a/lmd200.c +++ b/lmd200.c @@ -87,8 +87,8 @@ static void storeData(pSICSOBJ self, int start, char *lineBuffer) pPtr = trim(pPtr); pPtr = stptok(pPtr, number, 80, " "); while (pPtr != NULL) { - if (strstr(number, "noinp") != NULL) { - node->value.v.floatArray[start] = -9999.99; + if (strstr(number, "noinp") != NULL || strstr(number,"----") != NULL) { + node->value.v.floatArray[start] = -99.99; start++; } else { pKomma = strchr(number, ','); @@ -134,30 +134,38 @@ static int countTokens(char *lineBuffer) /*--------------------------------------------------------------------------*/ static void interpretLine(pSICSOBJ self, pLMD200 priv) { - pHdb node = NULL; - + pHdb node = NULL, stamp = NULL; + hdbValue timeVal; + node = GetHipadabaNode(self->objectNode, "data"); assert(node != NULL); - /* - * TODO: check the message headers against what is now coming out - * of the device. - */ + /* printf("Interpreting line: %s\n", priv->lineBuffer); */ if (strstr(priv->lineBuffer, "Alarm") != NULL) { doAlarm(self, priv->lineBuffer); return; } else if (strstr(priv->lineBuffer, "Pre") != NULL) { setAlarm(self, priv->lineBuffer); - } else if (strstr(priv->lineBuffer, "...0") == NULL) { + } else if (strstr(priv->lineBuffer, "001...008") != NULL) { storeData(self, 0, priv->lineBuffer); NotifyHipadabaPar(node, NULL); return; - } else if(strstr(priv->lineBuffer,"..000") == NULL) { + } else if(strstr(priv->lineBuffer,"009...016") != NULL) { storeData(self, 8, priv->lineBuffer); NotifyHipadabaPar(node, NULL); + } else if(strstr(priv->lineBuffer,"017...019") != NULL) { + storeData(self, 16, priv->lineBuffer); + NotifyHipadabaPar(node, NULL); + } else { + if(strlen(priv->lineBuffer) > 5){ + timeVal = MakeHdbText(priv->lineBuffer); + stamp = GetHipadabaNode(self->objectNode, "timestamp"); + if(stamp != NULL){ + UpdateHipadabaPar(stamp, timeVal, NULL); + } + } } } - /*----------------------------------------------------------------------------*/ static int LMD200Callback(int handle, void *userData) { @@ -168,13 +176,19 @@ static int LMD200Callback(int handle, void *userData) pPtr = ANETreadPtr(priv->asChannel,&length); pEnd = strchr(pPtr,(int)'\r'); - if(pEnd != NULL){ - length = pEnd - pPtr; + while(pEnd != NULL){ *pEnd = '\0'; - strncpy(priv->lineBuffer,pPtr,BUFLEN); - ANETreadConsume(priv->asChannel,length+1); + strcpy(priv->lineBuffer,pPtr); interpretLine(self,priv); + ANETreadConsume(priv->asChannel,(pEnd - pPtr)+1); + pPtr = ANETreadPtr(priv->asChannel,&length); + if(length > 0){ + pEnd = strchr(pPtr,(int)'\r'); + } else { + pEnd = NULL; + } } + return 1; } /*---------------------------------------------------------------------------*/ @@ -206,7 +220,7 @@ int MakeLMD200(SConnection * pCon, SicsInterp * pSics, void *pData, priv->port = atoi(argv[3]); priv->asChannel = ANETconnect(priv->host, priv->port); if (priv->asChannel < 0) { - SCWrite(pCon, "ERROR: failed to connect to LMD200", eError); + SCWrite(pCon, "ERROR: failed to connect to LMD200", eLogError); killLMD200(priv); return 0; } @@ -218,6 +232,7 @@ int MakeLMD200(SConnection * pCon, SicsInterp * pSics, void *pData, textValue = MakeHdbText("I do not feel very alarmed"); dataValue = makeHdbValue(HIPFLOATAR, 32); AddSICSHdbPar(self->objectNode, "alarm", usInternal, textValue); + AddSICSHdbPar(self->objectNode, "timestamp", usInternal, textValue); AddSICSHdbPar(self->objectNode, "data", usInternal, dataValue); self->pPrivate = priv; self->KillPrivate = killLMD200; diff --git a/make_gen b/make_gen index 1331b96..db48f44 100644 --- a/make_gen +++ b/make_gen @@ -23,7 +23,7 @@ OBJ=psi.o buffer.o ruli.o dmc.o nxsans.o nextrics.o sps.o pimotor.o \ dgrambroadcast.o sinq.o tabledrive.o tcpdocho.o julcho.o \ ritastorage.o poldizug.o audinelib.o delcam.o el737hpdrivsps.o \ rebin.o sanslirebin.o lmd200.o slsvme.o julprot.o sinqhttpprot.o \ - pmacprot.o pfeifferprot.o termprot.o phytron.o + pmacprot.o pfeifferprot.o termprot.o phytron.o autowin.o .SECONDARY.: sanslirebin.c diff --git a/psi.c b/psi.c index ff49ded..aaf1bf1 100644 --- a/psi.c +++ b/psi.c @@ -56,6 +56,7 @@ #include "amorset.h" #include "sinqhttp.h" #include "poldizug.h" +#include "autowin.h" /* from tcpdornier.c */ @@ -157,6 +158,7 @@ static void AddPsiCommands(SicsInterp * pInter) AddCommand(pInter, "MakePoldiReiss", MakePoldiReiss, NULL, NULL); AddCommand(pInter, "MakeSansliRebin", MakeSansliRebin, NULL, NULL); AddCommand(pInter, "MakeLMD200", MakeLMD200, NULL, NULL); + AddCommand(pInter, "MakeRitaWin", MakeRitaWin, NULL, NULL); /* AddCommand(pInter,"MakeDifrac",MakeDifrac,NULL,NULL); */ @@ -196,6 +198,7 @@ static void RemovePsiCommands(SicsInterp * pSics) RemoveCommand(pSics, "MakeSinq"); RemoveCommand(pSics, "MakeTableDrive"); RemoveCommand(pSics, "MakeAmorSet"); + RemoveCommand(pSics, "MakeRitaWin"); } /*---------------------------------------------------------------------*/ diff --git a/tdchm.c b/tdchm.c index f541d18..ac241c8 100644 --- a/tdchm.c +++ b/tdchm.c @@ -598,6 +598,7 @@ pHistDriver MakeTDCHM(pStringDict pOption) pNew->GetMonitor = TDCGetMonitor; pNew->GetTime = TDCGetTime; pNew->Preset = TDCPreset; + pNew->SubSample = DefaultSubSample; pNew->FreePrivate = TDCFree; pNew->Pause = TDCPause; pNew->Continue = TDCContinue;