- Added an autoamtic window calculation routine for RITA

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 <math.h>
+#include <sics.h>
+#include <sicsobj.h>
+#include <HistMem.h>
+#include <HistMem.i>
+#include <HistDriv.i>
+#include <lld.h>
+#include <dynstring.h>
+#include <sicshipadaba.h>
+
+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;
+}