sics/fomerge.c

/*---------------------------------------------------------------------------
                           F O M E R G E

 This module is for the FOCUS instrument. FOCUS has three detector 
 banks. These are sent to SICS as one big array by the histogram memory.
 However, depending on the situation these need to be accessed either:
 - as single banks.
 - or as a merged bank, where detector counts from all three banks are
   summed onto a common two theta scale. 

 This module is responsible for accessing single banks and for doing
 the merging operation.

 copyright: see copyright.h


 Mark Koennecke, March 2000
 
 extended to support nxscripted file writing: Mark Koennecke, May 2004
 
 extended to support GTSE, Mark Koennecke, May 2008
  
 modifed to support second generation HM's via the bridge, 
 Mark Koennecke, December 2012
 --------------------------------------------------------------------------*/
#include <stdlib.h>
#include <assert.h>
#include <stdio.h>
#include "sics.h"
#include "nxscript.h"
#include "HistMem.h"
#include "fortify.h"
#include "scan.h"
#include "fitcenter.h"
#include "sicsdata.h"

static pFit fitter = NULL;


#include "fomerge.h"


/*================== module static data items ==========================*/
static HistInt *masterData, *upperData, *mediumData, *lowerData,
    *mergedData = NULL;
static int timeBin, nUpper, nLower, nMedium, nMerged;
static float *upperTheta, *lowerTheta, *mediumTheta, *mergedTheta;
static int *mergeUp, *mergeMed, *mergeLow;
int iMerged;

int medium, upper, lower;
/*------------------------------------------------------------------------
   The tables with the theta values for the detector banks and the merge
   table living in mergeUp, mergeMed, mergeLow is initialized from a
   data file, the mergefile. mergeUp, mergeLow, mergeMed contain a 0 if
   the detector from this bank does not need to be summed for this merged
   two theta value or the number of the detector from  this bank to use. 
   -------------------------------------------------------------------------*/
int initializeFM(char *mergefile)
{
  FILE *fd = NULL;
  int i, j, iTest;
  char pBuffer[132];
  char *pPtr = NULL;

  /* open the file */
  fd = fopen(mergefile, "r");
  if (!fd) {
    fprintf(stdout, "ERROR: merge data file not found!\n");
    return 0;
  }

  /* skip first line */
  pPtr = fgets(pBuffer, 131, fd);

  /* read the number of merged data points */
  fgets(pBuffer, 131, fd);
  iTest = sscanf(pBuffer, "%d", &nMerged);
  if ((iTest != 1) || (nMerged <= 0)) {
    fprintf(stdout, "ERROR: Invalid or corrupt merge dat file");
    return 0;
  }

  /* allocate space */
  mergedTheta = (float *) malloc(nMerged * sizeof(float));
  mergeUp = (int *) malloc(nMerged * sizeof(int));
  mergeMed = (int *) malloc(nMerged * sizeof(int));
  mergeLow = (int *) malloc(nMerged * sizeof(int));
  if (mergedTheta == NULL || mergeUp == NULL || mergeMed == NULL
      || mergeLow == NULL) {
    fprintf(stdout, "ERROR: out of memory in fomerge\n");
    return 0;
  }

  /* read the merge table */
  for (i = 0; i < nMerged; i++) {
    fgets(pBuffer, 131, fd);
    iTest = sscanf(pBuffer, "%d %f %d %d %d", &j, &mergedTheta[i],
                   &mergeUp[i], &mergeMed[i], &mergeLow[i]);
    if (iTest != 5) {
      fprintf(stdout, "ERROR: Invalid or corrupt merge dat file");
      return 0;
    }
  }

  /* read upper bank two theta */
  fgets(pBuffer, 131, fd);
  fgets(pBuffer, 131, fd);
  iTest = sscanf(pBuffer, "%d", &nUpper);
  if ((iTest != 1) || (nUpper <= 0)) {
    fprintf(stdout, "ERROR: Invalid or corrupt merge data  file");
    return 0;
  }
  upperTheta = (float *) malloc(nUpper * sizeof(float));
  if (upperTheta == NULL) {
    fprintf(stdout, "ERROR: out of memory in fomerge\n");
    return 0;
  }
  for (i = 0; i < nUpper; i++) {
    fgets(pBuffer, 131, fd);
    sscanf(pBuffer, "%f", &upperTheta[i]);
  }

  /* read middle bank two theta */
  fgets(pBuffer, 131, fd);
  fgets(pBuffer, 131, fd);
  iTest = sscanf(pBuffer, "%d", &nMedium);
  if ((iTest != 1) || (nMedium <= 0)) {
    fprintf(stdout, "ERROR: Invalid or corrupt merge data file");
    return 0;
  }
  mediumTheta = (float *) malloc(nMedium * sizeof(float));
  if (mediumTheta == NULL) {
    fprintf(stdout, "ERROR: out of memory in fomerge\n");
    return 0;
  }
  for (i = 0; i < nMedium; i++) {
    fgets(pBuffer, 131, fd);
    sscanf(pBuffer, "%f", &mediumTheta[i]);
  }

  /* read lower bank two theta */
  fgets(pBuffer, 131, fd);
  fgets(pBuffer, 131, fd);
  iTest = sscanf(pBuffer, "%d", &nLower);
  if ((iTest != 1) || (nLower <= 0)) {
    fprintf(stdout, "ERROR: Invalid or corrupt merge data file");
    return 0;
  }
  lowerTheta = (float *) malloc(nLower * sizeof(float));
  if (lowerTheta == NULL) {
    fprintf(stdout, "ERROR: out of memory in fomerge\n");
    return 0;
  }
  for (i = 0; i < nLower; i++) {
    fgets(pBuffer, 131, fd);
    sscanf(pBuffer, "%f", &lowerTheta[i]);
  }

  /* done */
  timeBin = 0;
  iMerged = 0;
  medium = 1;
  upper = lower = 1;

  fclose(fd);
  return 1;
}

/*--------------------------------------------------------------------------
   sets the data pointer and finds the starting points for the different 
   banks. This routine contains the knowledge how the detector banks are 
   laid out in the histogram returned from the histogram memory.
   -------------------------------------------------------------------------*/
int setFMDataPointer(HistInt * lData, int mytimeBins, int bank)
{
  HistInt *dataPtr;

  if (lData == NULL)
    return 0;

  /* the first set is the medium bank */
  masterData = dataPtr = lData;
  if (bank == 2) {
    if (medium) {
      mediumData = masterData;
    } else {
      mediumData = NULL;
    }
  }
  /* the next set is the upper bank */
  if (bank == 1) {
    if (upper) {
      upperData = masterData;
    } else {
      upperData = NULL;
    }
  }
  /* the last is the lower bank */
  if (bank == 3) {
    if (lower) {
      lowerData = masterData;
    } else {
      lowerData = NULL;
    }
  }
  /* deal with data allocation for merged data */
  if (mytimeBins != timeBin) {
    if (mergedData != NULL) {
      free(mergedData);
    }
    timeBin = mytimeBins;
    mergedData = (HistInt *) malloc(nMerged * timeBin * sizeof(HistInt));
    if (mergedData == NULL) {
      fprintf(stdout,
              "ERROR: out of memory in fomerge.setFMdataPointer\n");
      return 0;
    }
  }

  /*
     clear mergedData to 0 
   */
  memset(mergedData, 0, nMerged * timeBin * sizeof(HistInt));
  iMerged = 0;
  return 1;
}

/*-----------------------------------------------------------------------*/
void setFMconfiguration(int up, int med, int low)
{
  upper = up;
  medium = med;
  lower = low;
}

/*--------------------------------------------------------------------------
  The algorithm is suboptimal because we have to take care of the case of
  missing detector banks due to missing or broken electronics or whatever.
*/
static void mergeData(void)
{
  int i, j, nDiv;
  HistInt *startMerge, *startData;

  if (!mergedData)
    return;

  memset(mergedData, 0, nMerged * sizeof(HistInt));
  for (i = 0; i < nMerged; i++) {
    startMerge = mergedData + i * timeBin;
    nDiv = 0;
    /* upper bank contribution */
    if (mergeUp[i] != 0 && upperData != NULL) {
      startData = upperData + (mergeUp[i] - 1) * timeBin;
      for (j = 0; j < timeBin; j++) {
        startMerge[j] += startData[j];
      }
      nDiv++;
    }
    /* medium bank contribution */
    if (mergeMed[i] != 0 && mediumData != NULL) {
      startData = mediumData + (mergeMed[i] - 1) * timeBin;
      for (j = 0; j < timeBin; j++) {
        startMerge[j] += startData[j];
      }
      nDiv++;
    }
    /* lower bank contribution */
    if (mergeLow[i] != 0 && lowerData != NULL) {
      startData = lowerData + (mergeLow[i] - 1) * timeBin;
      for (j = 0; j < timeBin; j++) {
        startMerge[j] += startData[j];
      }
      nDiv++;
    }
    /* do we need to run a division? */
    if (nDiv > 1) {
      for (j = 0; j < timeBin; j++) {
        startMerge[j] /= nDiv;
      }
    }
  }                             /* end of for */
  iMerged = 1;
}

/*-----------------------------------------------------------------------*/
void killFM(void)
{
  if (mergedData != NULL)
    free(mergedData);

  if (lowerTheta != NULL)
    free(lowerTheta);

  if (mediumTheta != NULL)
    free(mediumTheta);

  if (upperTheta != NULL)
    free(upperTheta);

  if (mergedTheta != NULL)
    free(mergedTheta);

  if (mergeUp != NULL)
    free(mergeUp);

  if (mergeMed != NULL)
    free(mergeMed);

  if (mergeLow != NULL)
    free(mergeLow);

}

/*--------------------------------------------------------------------------
  Below are some not very interesting data access routines 
  ---------------------------------------------------------------------------*/
HistInt *getFMBankPointer(int which)
{
  switch (which) {
  case UPPER:
    return upperData;
    break;
  case MIDDLE:
    return mediumData;
    break;
  case LOWER:
    return lowerData;
    break;
  case MERGED:
    if (!iMerged)
      mergeData();
    return mergedData;
    break;
  default:
    assert(0);
  }
}

/*-------------------------------------------------------------------------*/
float *getFMBankTheta(int which)
{
  switch (which) {
  case UPPER:
    return upperTheta;
    break;
  case MIDDLE:
    return mediumTheta;
    break;
  case LOWER:
    return lowerTheta;
    break;
  case MERGED:
    return mergedTheta;
    break;
  default:
    assert(0);
  }
}

/*-------------------------------------------------------------------------*/
int getFMdim(int which)
{
  switch (which) {
  case UPPER:
    return nUpper;
    break;
  case MIDDLE:
    return nMedium;
    break;
  case LOWER:
    return nLower;
    break;
  case MERGED:
    return nMerged;
    break;
  case TIMEBIN:
    return timeBin;
    break;
  default:
    assert(0);
  }
}

/*---------------------------------------------------------------------*/
int InstallFocusMerge(SConnection * pCon, SicsInterp * pSics, void *pData,
                      int argc, char *argv[])
{
  int status;
  char pBueffel[256];

  if (argc < 2) {
    SCWrite(pCon, "ERROR: Insufficient arguments to InstallFocusMerge",
            eError);
    return 0;
  }

  status = initializeFM(argv[1]);
  if (!status) {
    snprintf(pBueffel, 255, "ERROR: failed to read mergefile %s", argv[1]);
    SCWrite(pCon, pBueffel, eError);
    return 0;
  }

  /* Install command */
  AddCommand(pSics, "focusmerge", FocusMergeAction, NULL, NULL);
  return 1;
}

/*---------------------------------------------------------------------*/
static pNXScript checkNXScript(SicsInterp * pSics, char *name)
{
  pNXScript result = NULL;

  result = FindCommandData(pSics, name, "NXScript");
  if (result == NULL) {
    return NULL;
  }
  if (result->fileHandle == NULL || result->dictHandle == NULL) {
    return NULL;
  }
  return result;
}

/*---------------------------------------------------------------------*/
static int updateHMFMData(SicsInterp * pSics, SConnection * pCon)
{
  int status, iTime;
  const float *fTimeBin = NULL;
  HistInt *data = NULL;
  pHistMem pMem = NULL;

  pMem = (pHistMem) FindHM(pSics, "hm2");
  if (pMem == NULL) {
    return 0;
  }
  fTimeBin = GetHistTimeBin(pMem, &iTime);
  setFMDataPointer(GetHistogramPointer(pMem, pCon), iTime, MIDDLE);

  pMem = (pHistMem) FindHM(pSics, "hm1");
  if (pMem == NULL) {
    return 0;
  }
  setFMDataPointer(GetHistogramPointer(pMem, pCon), iTime, LOWER);

  pMem = (pHistMem) FindHM(pSics, "hm3");
  if (pMem == NULL) {
    return 0;
  }
  setFMDataPointer(GetHistogramPointer(pMem, pCon), iTime, UPPER);
  setFMconfiguration(1, 1, 1);
  return 1;
}

/*-------------------------------------------------------------------*/
static int *calculateDetSum(HistInt * data, int iDet, int iTime)
{
  int i, j, iIndex;
  int *sum = NULL;

  sum = (int *) malloc(iDet * sizeof(int));
  if (!sum) {
    return NULL;
  }
  memset(sum, 0, iDet * sizeof(int));

  for (i = 0; i < iDet; i++) {
    iIndex = i * iTime;
    for (j = 0; j < iTime; j++) {
      sum[i] += data[iIndex + j];
    }
  }
  return sum;
}

/*-------------------------------------------------------------------*/
static int *calculateTimeSum(HistInt * data, int iDet, int iTime)
{
  int i, j;
  int *sum = NULL;

  if(data == NULL){
	  return NULL;
  }
  
  sum = (int *) malloc(iTime * sizeof(int));
  if (!sum) {
    return NULL;
  }
  memset(sum, 0, iTime * sizeof(int));

  for (i = 0; i < iTime; i++) {
    for (j = 0; j < iDet; j++) {
      sum[i] += data[j * iTime + i];
    }
  }
  return sum;
}

/*--------------------------------------------------------------------*/
static void checkSum(HistInt * sum, int iDet, char *name,
                     SConnection * pCon)
{
  int i, count;
  char pBueffel[256];

  for (i = 0, count = 0; i < iDet; i++) {
    if (sum[i] == 0) {
      count++;
    }
  }
  if (count == iDet || count == 0) {
    return;
  }
  snprintf(pBueffel, 255,
           "WARNING: %d of %d detectors in bank %s are empty", count, iDet,
           name);
  SCWrite(pCon, pBueffel, eWarning);
}

/*---------------------------------------------------------------------*/
static int putSum(SicsInterp * pSics, SConnection * pCon,
                  pNXScript nxscript, char *name, char *alias)
{
  HistInt *data = NULL;
  HistInt *sum = NULL;
  int iDet, iTime, i, j, iIndex, status;

  iTime = getFMdim(TIMEBIN);
  if (strcmp(name, "upper") == 0) {
    iDet = getFMdim(UPPER);
    data = getFMBankPointer(UPPER);
  } else if (strcmp(name, "middle") == 0) {
    iDet = getFMdim(MIDDLE);
    data = getFMBankPointer(MIDDLE);
  } else if (strcmp(name, "lower") == 0) {
    iDet = getFMdim(LOWER);
    data = getFMBankPointer(LOWER);
  } else if (strcmp(name, "merged") == 0) {
    iDet = getFMdim(MERGED);
    data = getFMBankPointer(MERGED);
  } else {
    SCWrite(pCon, "ERROR: detector bank to sum  not recognised", eError);
    return NX_ERROR;
  }

  sum = calculateDetSum(data, iDet, iTime);
  if (!sum) {
    SCWrite(pCon, "ERROR: out of memory summing bank", eError);
    return NX_ERROR;
  }
  checkSum(sum, iDet, name, pCon);

  status = NXDputalias(nxscript->fileHandle, nxscript->dictHandle,
                       alias, sum);
  free(sum);
  return status;
}

/*---------------------------------------------------------------------*/
static int TOFLambda(SicsInterp * pSics, SConnection * pCon,
                     int argc, char *argv[])
{
  int status, iTime, iDet, i;
  const float *fTimeBin = NULL;
  int *sum = NULL;
  long *lSum = NULL;
  pHistMem pMem = NULL;
  float fCenter, fFWHM, fStdDev, fVal;
  float fMon, fData, distMonoDet, distFermiDet, tdiff, lambda;

  pMem = (pHistMem) FindHM(pSics, "hm1");
  if (pMem == NULL) {
    SCWrite(pCon,
            "ERROR: need lower detector bank for lambda calculation",
            eError);
    return 0;
  }

  /**
   * locate elastic position in data
   */
  fTimeBin = GetHistTimeBin(pMem, &iTime);
  iDet = getFMdim(LOWER);
  sum = calculateTimeSum(GetHistogramPointer(pMem, pCon), iDet, iTime);
  if (!sum) {
    SCWrite(pCon, "ERROR: out of memory calculating lambda", eError);
    return 0;
  }
  if (fitter == NULL) {
    fitter = CreateFitCenter(NULL);
    if (!fitter) {
      SCWrite(pCon, "ERROR: cannot allocate fitting structure", eError);
      return 0;
    }
  }
  /*
     copy sum to make compiler happy
   */
  lSum = (long *) malloc(iTime * sizeof(long));
  if (lSum == NULL) {
    SCWrite(pCon, "ERROR: out of memory in TOFLambda", eError);
    free(sum);
    return 0;
  }
  for (i = 0; i < iTime; i++) {
    lSum[i] = sum[i];
  }
  status = CalculateFitFromData(fitter, (float *) fTimeBin, lSum, iTime);
  if (status < 0) {
    SCWrite(pCon, "ERROR: no peak in data", eError);
    free(sum);
    free(lSum);
    return 0;
  }
  GetFitResults(fitter, &fCenter, &fStdDev, &fFWHM, &fVal);
  fData = fCenter;

  /*
   * locate elastic position in tofmon
   */
  GetHistogram(pMem, pCon, 0, iTime * iDet, iTime * (iDet + 1),
               sum, iTime * sizeof(HistInt));
  for (i = 0; i < iTime; i++) {
    lSum[i] = sum[i];
  }
  status = CalculateFitFromData(fitter, (float *) fTimeBin, lSum, iTime);
  GetFitResults(fitter, &fCenter, &fStdDev, &fFWHM, &fVal);
  fMon = fCenter;
  free(sum);
  free(lSum);

  /*
   * calculate
   */
  distFermiDet = 3000.;
  distMonoDet = distFermiDet - 215.7;
  tdiff = fData - fMon;
  lambda = tdiff / (252.78 * distMonoDet * .001);
  SCPrintf(pCon, eValue, "toflambda = %f", lambda);

  return 1;
}

/*---------------------------------------------------------------------*/
static float calcElastic(SicsInterp * pSics, SConnection * pCon)
{
  int status, iTime, iDet, i;
  const float *fTimeBin = NULL;
  int *sum = NULL;
  long *lSum = NULL;
  pHistMem pMem = NULL;
  float fCenter, fFWHM, fStdDev, fVal;

  pMem = (pHistMem) FindHM(pSics, "hm2");
  if (pMem == NULL) {
    SCWrite(pCon,
            "ERROR: need middle detector bank for elastic peak calculation",
            eError);
    return -1.;
  }
  fTimeBin = GetHistTimeBin(pMem, &iTime);
  iDet = getFMdim(MIDDLE);
  sum = calculateTimeSum(GetHistogramPointer(pMem, pCon), iDet, iTime);
  if (!sum) {
    SCWrite(pCon, "ERROR: out of memory calculating elastic peak position",
            eError);
    return -1;
  }
  if (fitter == NULL) {
    fitter = CreateFitCenter(NULL);
    if (!fitter) {
      SCWrite(pCon, "ERROR: cannot allocate fitting structure", eError);
      return -1.;
    }
  }
  /*
     copy sum to make compiler happy
   */
  lSum = (long *) malloc(iTime * sizeof(long));
  if (lSum == NULL) {
    SCWrite(pCon, "ERROR: out of memory in putElastic", eError);
    free(sum);
    return -1.;
  }
  for (i = 0; i < iTime; i++) {
    lSum[i] = sum[i];
  }
  status = CalculateFitFromData(fitter, (float *) fTimeBin, lSum, iTime);
  free(lSum);
  GetFitResults(fitter, &fCenter, &fStdDev, &fFWHM, &fVal);
  free(sum);
  return fCenter;
}

/*---------------------------------------------------------------------*/
#define ABS(x) (x < 0 ? -(x) : (x))

static int putElastic(SicsInterp * pSics, SConnection * pCon,
                      pNXScript pNexus, char *alias, float fElastic)
{
  float fCalc;
  int status;

  fCalc = calcElastic(pSics, pCon);
  if (ABS(fElastic - fCalc) < 20) {
    fElastic = fCalc;
  }

  status = NXDputalias(pNexus->fileHandle, pNexus->dictHandle, alias,
                       &fElastic);
  return status;
}

/*----------------------------------------------------------------------*/
static int FMputTTH(SConnection * pCon, int argc, char *argv[])
{
  pSICSData data = NULL;
  int length = -1, i;
  float *tthData = NULL;

  if (argc < 4) {
    SCWrite(pCon, "ERROR: insufficient no of arguments to FMputTTH",
            eError);
    return 0;
  }

  data = (pSICSData) FindCommandData(pServ->pSics, argv[3], "SICSData");
  if (data == NULL) {
    SCWrite(pCon, "ERROR: SICSData object not found", eError);
    return 0;
  }
  if (strcmp(argv[2], "upper") == 0) {
    length = getFMdim(UPPER);
    tthData = getFMBankTheta(UPPER);
  } else if (strcmp(argv[2], "middle") == 0) {
    length = getFMdim(MIDDLE);
    tthData = getFMBankTheta(MIDDLE);
  } else if (strcmp(argv[2], "lower") == 0) {
    length = getFMdim(LOWER);
    tthData = getFMBankTheta(LOWER);
  } else if (strcmp(argv[2], "merged") == 0) {
    length = getFMdim(MERGED);
    tthData = getFMBankTheta(MERGED);
  } else {
    SCWrite(pCon, "ERROR: requested two_theta for invalid detector bank",
            eError);
    return 0;
  }
  if (length < 0 || tthData == NULL) {
    SCWrite(pCon, "ERROR: requested two_theta for invalid detector bank",
            eError);
    return 0;
  }
  clearSICSData(data);
  for (i = 0; i < length; i++) {
    setSICSDataFloat(data, i, tthData[i]);
  }
  SCSendOK(pCon);
  return 1;
}

/*---------------------------------------------------------------------*/
static int FMcopyMerged(SConnection * pCon, int argc, char *argv[])
{
  pSICSData data = NULL;
  int i, length;
  HistInt *hmData = NULL;

  if (argc < 3) {
    SCWrite(pCon, "ERROR: insufficient no of arguments to FMcopyMerged",
            eError);
    return 0;
  }

  data = (pSICSData) FindCommandData(pServ->pSics, argv[2], "SICSData");
  if (data == NULL) {
    SCWrite(pCon, "ERROR: SICSData object not found", eError);
    return 0;
  }
  if (!updateHMFMData(pServ->pSics, pCon)) {
    SCWrite(pCon,
            "ERROR: not enough HM's to merge or bad names in fomerge.c",
            eError);
    return 0;
  }

  clearSICSData(data);
  length = getFMdim(MERGED) * getFMdim(TIMEBIN);
  hmData = getFMBankPointer(MERGED);
  if (hmData == NULL) {
    SCWrite(pCon, "ERROR: merged data not available", eError);
    return 0;
  }
  for (i = 0; i < length; i++) {
    setSICSDataInt(data, i, hmData[i]);
  }
  SCSendOK(pCon);
  return 1;
}

/*---------------------------------------------------------------------*/
static int FMcopyMergedSum(SConnection * pCon, int argc, char *argv[])
{
  pSICSData data = NULL;
  int i, length, tbin, j, row;
  HistInt *hmData = NULL, *sumData = NULL;

  if (argc < 3) {
    SCWrite(pCon, "ERROR: insufficient no of arguments to FMcopyMerged",
            eError);
    return 0;
  }

  data = (pSICSData) FindCommandData(pServ->pSics, argv[2], "SICSData");
  if (data == NULL) {
    SCWrite(pCon, "ERROR: SICSData object not found", eError);
    return 0;
  }
  if (!updateHMFMData(pServ->pSics, pCon)) {
    SCWrite(pCon,
            "ERROR: not enough HM's to merge or bad names in fomerge.c",
            eError);
    return 0;
  }

  clearSICSData(data);
  length = getFMdim(MERGED);
  tbin = getFMdim(TIMEBIN);
  hmData = getFMBankPointer(MERGED);
  if (hmData == NULL) {
    SCWrite(pCon, "ERROR: merged data not available", eError);
    return 0;
  }
  sumData = malloc(tbin * sizeof(int));
  if (sumData == NULL) {
    SCWrite(pCon, "ERROR: out-of-memory in FMcopyMergedSum", eError);
    return 0;
  }
  memset(sumData, 0, tbin * sizeof(int));
  for (j = 0; j < length; j++) {
    row = j * tbin;
    for (i = 0; i < tbin; i++) {
      sumData[i] += hmData[row + i];
    }
  }
  for (i = 0; i < tbin; i++) {
    setSICSDataInt(data, i, sumData[i]);
  }
  free(sumData);
  SCSendOK(pCon);
  return 1;
}

/*-----------------------------------------------------------------------
 Usage:
       focusmerge puttwotheta nxscriptmod bankname alias
       focusmerge puttth bankname sicsdataname
       focusmerge copymerged sicsdataname
       focusmerge copymergedsum sicsdataname
       focusmerge putmerged nxscriptmod alias
       focusmerge putsum nxscriptmod bankname alias
       focusmerge putelastic nxscriptmod alias theoelastic
       focusmerge toflambda
       focusmerge elastic 

nxscriptmod = name of the nxscript module used for writing, must be open
alias       = The alias under which to write the data item
theoelastic = theoretical elastic peak position
------------------------------------------------------------------------*/
int FocusMergeAction(SConnection * pCon, SicsInterp * pSics, void *pData,
                     int argc, char *argv[])
{
  int status;
  pNXScript pNexus = NULL;
  float fElastic;
  char pNum[20];

  if (argc < 2) {
    SCWrite(pCon, "ERROR: Insufficient arguments to focusmerge", eError);
    return 0;
  }

  strtolower(argv[1]);

  if (strcmp(argv[1], "puttth") == 0) {
    return FMputTTH(pCon, argc, argv);
  }

  if (strcmp(argv[1], "copymerged") == 0) {
    return FMcopyMerged(pCon, argc, argv);
  }
  if (strcmp(argv[1], "copymergedsum") == 0) {
    return FMcopyMergedSum(pCon, argc, argv);
  }
  if (strcmp(argv[1], "toflambda") == 0) {
    return TOFLambda(pSics, pCon, argc, argv);
  }
  if (strcmp(argv[1], "elastic") == 0) {
    fElastic = calcElastic(pSics, pCon);
    SCPrintf(pCon, eValue, "tofelastic = %f", fElastic);
    return 1;
  }

  if (strcmp(argv[1], "puttwotheta") == 0) {
    if (argc < 4) {
      SCWrite(pCon,
              "ERROR: Insufficient arguments to focusmerge puttwotheta",
              eError);
      return 0;
    }
    pNexus = checkNXScript(pSics, argv[2]);
    if (pNexus == NULL) {
      SCWrite(pCon, "ERROR: bad nxscript name or NeXus file not open",
              eError);
      return 0;
    }
    strtolower(argv[3]);
    if (strcmp(argv[3], "upper") == 0) {
      status = NXDputalias(pNexus->fileHandle, pNexus->dictHandle,
                           argv[4], getFMBankTheta(UPPER));
    } else if (strcmp(argv[3], "middle") == 0) {
      status = NXDputalias(pNexus->fileHandle, pNexus->dictHandle,
                           argv[4], getFMBankTheta(MIDDLE));
    } else if (strcmp(argv[3], "lower") == 0) {
      status = NXDputalias(pNexus->fileHandle, pNexus->dictHandle,
                           argv[4], getFMBankTheta(LOWER));
    } else if (strcmp(argv[3], "merged") == 0) {
      status = NXDputalias(pNexus->fileHandle, pNexus->dictHandle,
                           argv[4], getFMBankTheta(MERGED));
    } else {
      SCWrite(pCon, "ERROR: requested two_theta for invalid detector bank",
              eError);
      return 0;
    }
    if (status == NX_OK) {
      SCSendOK(pCon);
      return 1;
    } else {
      SCWrite(pCon, "ERROR: failed to write two theta array to file",
              eError);
      return 0;
    }
  } else if (strcmp(argv[1], "putmerged") == 0) {
    if (argc < 4) {
      SCWrite(pCon, "ERROR: Insufficient arguments to focusmerge", eError);
      return 0;
    }
    pNexus = checkNXScript(pSics, argv[2]);
    if (pNexus == NULL) {
      SCWrite(pCon, "ERROR: bad nxscript name or NeXus file not open",
              eError);
      return 0;
    }
    if (!updateHMFMData(pSics, pCon)) {
      SCWrite(pCon,
              "ERROR: not enough HM's to merge or bad names in fomerge.c",
              eError);
      return 0;
    }
    snprintf(pNum, 19, "%d", getFMdim(MERGED));
    NXDupdate(pNexus->dictHandle, "noofdetectors", pNum);
    snprintf(pNum, 19, "%d", getFMdim(TIMEBIN));
    NXDupdate(pNexus->dictHandle, "timebin", pNum);
    status = NXDputalias(pNexus->fileHandle, pNexus->dictHandle,
                         argv[3], getFMBankPointer(MERGED));
    if (status == NX_OK) {
      SCSendOK(pCon);
      return 1;
    } else {
      SCWrite(pCon, "ERROR: failed to write merged data to file", eError);
      return 0;
    }

  } else if (strcmp(argv[1], "putsum") == 0) {
    if (argc < 4) {
      SCWrite(pCon, "ERROR: Insufficient arguments to focusmerge putsum",
              eError);
      return 0;
    }
    pNexus = checkNXScript(pSics, argv[2]);
    if (pNexus == NULL) {
      SCWrite(pCon, "ERROR: bad nxscript name or NeXus file not open",
              eError);
      return 0;
    }
    updateHMFMData(pSics, pCon);
    status = putSum(pSics, pCon, pNexus, argv[3], argv[4]);
    if (status == NX_OK) {
      SCSendOK(pCon);
      return 1;
    } else {
      SCWrite(pCon, "ERROR: failed to write summed data to file", eError);
      return 0;
    }
  } else if (strcmp(argv[1], "putelastic") == 0) {
    if (argc < 4) {
      SCWrite(pCon,
              "ERROR: Insufficient arguments to focusmerge putelastic",
              eError);
      return 0;
    }
    pNexus = checkNXScript(pSics, argv[2]);
    if (pNexus == NULL) {
      SCWrite(pCon, "ERROR: bad nxscript name or NeXus file not open",
              eError);
      return 0;
    }
    fElastic = atof(argv[4]);
    status = putElastic(pSics, pCon, pNexus, argv[3], fElastic);
    if (status == NX_OK) {
      SCSendOK(pCon);
      return 1;
    } else {
      SCWrite(pCon, "ERROR: failed to write elastic peak position",
              eError);
      return 0;
    }
  } else {
    SCWrite(pCon, "ERROR: subcommand to focusmerge not understood",
            eError);
    return 0;
  }
  return 0;
}