Gatestrain release

lalapps/src/pulsar/MakeData/ComputePSDFunction.c

+/*
+*  Copyright (C) 2007, 2010 Karl Wette
+*  Copyright (C) 2007 Badri Krishnan, Iraj Gholami, Reinhard Prix, Alicia Sintes
+*  Copyright (C) 2017-2021 Benjamin Steltner
+*
+*  This program is free software; you can redistribute it and/or modify
+*  it under the terms of the GNU General Public License as published by
+*  the Free Software Foundation; either version 2 of the License, or
+*  (at your option) any later version.
+*
+*  This program is distributed in the hope that it will be useful,
+*  but WITHOUT ANY WARRANTY; without even the implied warranty of
+*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*  GNU General Public License for more details.
+*
+*  You should have received a copy of the GNU General Public License
+*  along with with program; see the file COPYING. If not, write to the
+*  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
+*  MA  02111-1307  USA
+*/
+
+/**
+ * \file
+ * \ingroup lalapps_pulsar_SFTTools
+ * \author Badri Krishnan, Iraj Gholami, Reinhard Prix, Alicia Sintes, Karl Wette, Benjamin Steltner
+ * \brief Compute power spectral densities
+ */
+
+/**
+ * This file is a copy of ComputePSD.c ( https://lscsoft.docs.ligo.org/lalsuite/lalapps/_compute_p_s_d_8c_source.html ).
+ * Most file I/O is removed, the main method takes an in-memory SFT and computes & returns an in-memory PSD to the caller.
+ */
+
+#include <glob.h>
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <gsl/gsl_sort.h>
+#include <gsl/gsl_math.h>
+#include <lal/LALStdlib.h>
+#include <lal/LALConstants.h>
+#include <lal/AVFactories.h>
+#include <lal/SeqFactories.h>
+#include <lal/SFTfileIO.h>
+#include <lal/Random.h>
+#include <lal/PulsarDataTypes.h>
+#include <lal/UserInput.h>
+#include <lal/NormalizeSFTRngMed.h>
+#include <lal/LALInitBarycenter.h>
+#include <lal/SFTClean.h>
+#include <lal/Segments.h>
+#include <lal/FrequencySeries.h>
+#include <lal/Units.h>
+
+#include <lal/LogPrintf.h>
+
+#include <LALAppsVCSInfo.h>
+#include "ComputePSDFunction.h"
+
+/* ---------- Error codes and messages ---------- */
+#define COMPUTEPSDC_ENORM 0
+#define COMPUTEPSDC_ESUB  1
+#define COMPUTEPSDC_EARG  2
+#define COMPUTEPSDC_EBAD  3
+#define COMPUTEPSDC_EFILE 4
+#define COMPUTEPSDC_ENULL 5
+#define COMPUTEPSDC_EMEM  6
+
+#define COMPUTEPSDC_MSGENORM "Normal exit"
+#define COMPUTEPSDC_MSGESUB  "Subroutine failed"
+#define COMPUTEPSDC_MSGEARG  "Error parsing arguments"
+#define COMPUTEPSDC_MSGEBAD  "Bad argument values"
+#define COMPUTEPSDC_MSGEFILE "Could not create output file"
+#define COMPUTEPSDC_MSGENULL "Null Pointer"
+#define COMPUTEPSDC_MSGEMEM  "Out of memory"
+
+/*---------- local defines ---------- */
+
+#define TRUE (1==1)
+#define FALSE (1==0)
+
+/* ----- Macros ----- */
+
+/* ---------- local types ---------- */
+
+/** types of mathematical operations */
+
+
+/**
+ * Config variables 'derived' from user-input
+ */
+typedef struct
+{
+  UINT4 firstBin;		/**< first PSD bin for output */
+  UINT4 lastBin;		/**< last PSD bin for output */
+  LALSeg dataSegment;		/**< the data-segment for which PSD was computed */
+} ConfigVariables_t;
+
+/* ---------- global variables ----------*/
+extern int vrbflg;
+
+/* ---------- local prototypes ---------- */
+static REAL8 math_op(REAL8*, size_t, INT4);
+int XLALDumpMultiPSDVector ( const CHAR *outbname, const MultiPSDVector *multiPSDVect );
+
+int XLALCropMultiPSDandSFTVectors ( MultiPSDVector *multiPSDVect, MultiSFTVector *multiSFTVect, UINT4 firstBin, UINT4 lastBin );
+
+/*============================================================
+ * FUNCTION definitions
+ *============================================================*/
+int
+computePSD(PSDUserVarsTag uvar, SFTtype *sft, REAL8FrequencySeries **Fseries)  {
+  ConfigVariables_t XLAL_INIT_DECL(cfg);
+
+  uvar.PSDmthopSFTs = MATH_OP_HARMONIC_MEAN;
+  uvar.PSDmthopIFOs = MATH_OP_HARMONIC_SUM;
+
+  uvar.nSFTmthopSFTs = MATH_OP_ARITHMETIC_MEAN;
+  uvar.nSFTmthopIFOs = MATH_OP_MAXIMUM;
+
+  uvar.PSDmthopBins  = MATH_OP_ARITHMETIC_MEDIAN;
+  uvar.nSFTmthopBins = MATH_OP_MAXIMUM;
+
+  UINT4 k, numBins, numIFOs, maxNumSFTs, X, alpha;
+  REAL8 Freq0, dFreq, normPSD;
+  UINT4 finalBinSize, finalBinStep, finalNumBins;
+
+  REAL8Vector *overSFTs = NULL; /* one frequency bin over SFTs */
+  REAL8Vector *overIFOs = NULL; /* one frequency bin over IFOs */
+  REAL8Vector *finalPSD = NULL; /* math. operation PSD over SFTs and IFOs */
+  REAL8Vector *finalNormSFT = NULL; /* normalised SFT power */
+
+  vrbflg = 1;	/* verbose error-messages */
+
+  /* set LAL error-handler */
+  lal_errhandler = LAL_ERR_EXIT;
+
+  MultiSFTVector *inputSFTs = NULL;
+  UINT4Vector *uint4vector = NULL;
+  XLAL_CHECK ( (uint4vector = XLALCreateUINT4Vector(1)) != NULL, XLAL_ENOMEM);
+  uint4vector->length = 1;
+  uint4vector->data[0] = 1;
+
+  XLAL_CHECK ( (inputSFTs = XLALCreateMultiSFTVector(sft->data->length, uint4vector)) != NULL, XLAL_ENOMEM );
+
+  //Input SFT into the MultiSFTVector
+  //MultiSFTVector -> SFTVector ** -> SFTtype * -> COMPLEX8Sequence * -> COMPLEX8 *
+  COMPLEX8Sequence *tmp = inputSFTs->data[0]->data->data;     //Temporarily safe the pointer to the data
+  memcpy(inputSFTs->data[0]->data, sft, sizeof(*sft));        //Copy all fields of SFT into the SFTVector
+  inputSFTs->data[0]->data->data = tmp;                       //Restore the correct data pointer
+  memcpy(inputSFTs->data[0]->data->data->data, sft->data->data, sizeof(sft->data->data[0]) * sft->data->length);    //Copy the actual sft data
+
+  cfg.firstBin = (UINT4) 0;
+  cfg.lastBin = (UINT4) inputSFTs->data[0]->data[0].data->length - 1;
+
+  LogPrintf (LOG_DEBUG, "Computing spectrogram and PSD ... ");
+
+  /* get power running-median rngmed[ |data|^2 ] from SFTs */
+  MultiPSDVector *multiPSD = NULL;
+  XLAL_CHECK_MAIN( ( multiPSD = XLALNormalizeMultiSFTVect ( inputSFTs, uvar.blocksRngMed, NULL ) ) != NULL, XLAL_EFUNC);
+  /* restrict this PSD to just the "physical" band if requested using {--Freq, --FreqBand} */
+  if ( ( XLALCropMultiPSDandSFTVectors ( multiPSD, inputSFTs, cfg.firstBin, cfg.lastBin )) != XLAL_SUCCESS ) {
+    XLALPrintError ("%s: XLALCropMultiPSDandSFTVectors (inputPSD, inputSFTs, %d, %d) failed with xlalErrno = %d\n", __func__, cfg.firstBin, cfg.lastBin, xlalErrno );
+    return EXIT_FAILURE;
+  }
+
+  /* start frequency and frequency spacing */
+  Freq0 = multiPSD->data[0]->data[0].f0;
+  dFreq = multiPSD->data[0]->data[0].deltaF;
+
+  /* number of raw bins in final PSD */
+  numBins = multiPSD->data[0]->data[0].data->length;
+  if ( (finalPSD = XLALCreateREAL8Vector ( numBins )) == NULL ) {
+    LogPrintf (LOG_CRITICAL, "Out of memory!\n");
+    return EXIT_FAILURE;
+  }
+
+  /* number of IFOs */
+  numIFOs = multiPSD->length;
+  if ( (overIFOs = XLALCreateREAL8Vector ( numIFOs )) == NULL ) {
+    LogPrintf (LOG_CRITICAL, "Out of memory!\n");
+    return EXIT_FAILURE;
+  }
+
+  /* maximum number of SFTs */
+  maxNumSFTs = 0;
+  for (X = 0; X < numIFOs; ++X) {
+    maxNumSFTs = GSL_MAX(maxNumSFTs, multiPSD->data[X]->length);
+  }
+  if ( (overSFTs = XLALCreateREAL8Vector ( maxNumSFTs )) == NULL ) {
+    LogPrintf (LOG_CRITICAL, "Out of memory!\n");
+    return EXIT_FAILURE;
+  }
+
+  /* normalize rngmd(power) to get proper *single-sided* PSD: Sn = (2/Tsft) rngmed[|data|^2]] */
+  normPSD = 2.0 * dFreq;
+
+  /* loop over frequency bins in final PSD */
+  for (k = 0; k < numBins; ++k) {
+
+    /* loop over IFOs */
+    for (X = 0; X < numIFOs; ++X) {
+
+      /* number of SFTs for this IFO */
+      UINT4 numSFTs = multiPSD->data[X]->length;
+
+      /* copy PSD frequency bins and normalise multiPSD for later use */
+      for (alpha = 0; alpha < numSFTs; ++alpha) {
+      	multiPSD->data[X]->data[alpha].data->data[k] *= normPSD;
+       	overSFTs->data[alpha] = multiPSD->data[X]->data[alpha].data->data[k];
+      }
+
+      /* compute math. operation over SFTs for this IFO */
+      overIFOs->data[X] = math_op(overSFTs->data, numSFTs, uvar.PSDmthopSFTs);
+      if ( isnan( overIFOs->data[X] ) )
+        XLAL_ERROR ( EXIT_FAILURE, "Found Not-A-Number in overIFOs->data[X=%d] = NAN ... exiting\n", X );
+
+    } /* for IFOs X */
+
+    /* compute math. operation over IFOs for this frequency */
+    finalPSD->data[k] = math_op(overIFOs->data, numIFOs, uvar.PSDmthopIFOs);
+    if ( isnan ( finalPSD->data[k] ) )
+      XLAL_ERROR ( EXIT_FAILURE, "Found Not-A-Number in finalPSD->data[k=%d] = NAN ... exiting\n", k );
+
+  } /* for freq bins k */
+  LogPrintfVerbatim ( LOG_DEBUG, "done.\n");
+
+  /* compute normalised SFT power */
+  if (uvar.outputNormSFT) {
+    LogPrintf (LOG_DEBUG, "Computing normalised SFT power ... ");
+
+    if ( (finalNormSFT = XLALCreateREAL8Vector ( numBins )) == NULL ) {
+      LogPrintf (LOG_CRITICAL, "Out of memory!\n");
+      return EXIT_FAILURE;
+    }
+
+    /* loop over frequency bins in SFTs */
+    for (k = 0; k < numBins; ++k) {
+
+      /* loop over IFOs */
+      for (X = 0; X < numIFOs; ++X) {
+
+	/* number of SFTs for this IFO */
+	UINT4 numSFTs = inputSFTs->data[X]->length;
+
+	/* compute SFT power */
+	for (alpha = 0; alpha < numSFTs; ++alpha) {
+	  COMPLEX8 bin = inputSFTs->data[X]->data[alpha].data->data[k];
+	  overSFTs->data[alpha] = crealf(bin)*crealf(bin) + cimagf(bin)*cimagf(bin);
+	}
+
+	/* compute math. operation over SFTs for this IFO */
+	overIFOs->data[X] = math_op(overSFTs->data, numSFTs, uvar.nSFTmthopSFTs);
+	if ( isnan ( overIFOs->data[X] ))
+          XLAL_ERROR ( EXIT_FAILURE, "Found Not-A-Number in overIFOs->data[X=%d] = NAN ... exiting\n", X );
+
+      } /* over IFOs */
+
+      /* compute math. operation over IFOs for this frequency */
+      finalNormSFT->data[k] = math_op(overIFOs->data, numIFOs, uvar.nSFTmthopIFOs);
+      if ( isnan( finalNormSFT->data[k] ) )
+        XLAL_ERROR ( EXIT_FAILURE, "Found Not-A-Number in bin finalNormSFT->data[k=%d] = NAN ... exiting\n", k );
+
+    } /* over freq bins */
+    LogPrintfVerbatim ( LOG_DEBUG, "done.\n");
+  }
+
+  /* ---------- if user requested it, output complete MultiPSDVector over IFOs X, timestamps and freq-bins into ASCI file(s) */
+  if ( uvar.dumpMultiPSDVector ) {
+    if ( XLALDumpMultiPSDVector ( uvar.outputPSD, multiPSD ) != XLAL_SUCCESS ) {
+      XLALPrintError ("%s: XLALDumpMultiPSDVector() failed, xlalErrnor = %d\n", __func__, xlalErrno );
+      return EXIT_FAILURE;
+    }
+  } /* if uvar.dumpMultiPSDVector */
+
+  finalBinSize = 1;
+  finalBinStep = finalBinSize;
+
+  /* work out total number of bins */
+  finalNumBins = (UINT4)floor((numBins - finalBinSize) / finalBinStep) + 1;
+
+  /* write final PSD to file */
+  if( 1==0) { // DONT WRITE PSD TO FILE
+    FILE *fpOut = NULL;
+
+    if ((fpOut = fopen(uvar.outputPSD, "wb")) == NULL) {
+      LogPrintf ( LOG_CRITICAL, "Unable to open output file %s for writing...exiting \n", uvar.outputPSD );
+      return EXIT_FAILURE;
+    }
+
+    /* write header info in comments */
+
+    if ( XLAL_SUCCESS != XLALOutputVCSInfo(fpOut, lalAppsVCSInfoList, 0, "%% "))
+      XLAL_ERROR ( XLAL_EFUNC );
+
+    /* write column headings */
+    fprintf(fpOut,"%%%% columns:\n%%%% FreqBinStart");
+    if (uvar.outFreqBinEnd)
+      fprintf(fpOut," FreqBinEnd");
+    fprintf(fpOut," PSD");
+    if (uvar.outputNormSFT)
+      fprintf(fpOut," normSFTpower");
+    fprintf(fpOut,"\n");
+
+    LogPrintf(LOG_DEBUG, "Printing PSD to file ... ");
+    for (k = 0; k < finalNumBins; ++k) {
+      UINT4 b = k * finalBinStep;
+
+      REAL8 f0 = Freq0 + b * dFreq;
+      REAL8 f1 = f0 + finalBinStep * dFreq;
+      fprintf(fpOut, "%f", f0);
+      if (uvar.outFreqBinEnd)
+	    fprintf(fpOut, "   %f", f1);
+
+      REAL8 psd = math_op(&(finalPSD->data[b]), finalBinSize, uvar.PSDmthopBins);
+      if ( isnan ( psd ))
+        XLAL_ERROR ( EXIT_FAILURE, "Found Not-A-Number in psd[k=%d] = NAN ... exiting\n", k );
+
+      fprintf(fpOut, "   %e", psd);
+
+      if (uvar.outputNormSFT) {
+	      REAL8 nsft = math_op(&(finalNormSFT->data[b]), finalBinSize, uvar.nSFTmthopBins);
+    	if ( isnan ( nsft ))
+          XLAL_ERROR ( EXIT_FAILURE, "Found Not-A-Number in nsft[k=%d] = NAN ... exiting\n", k );
+
+      	fprintf(fpOut, "   %f", nsft);
+      }
+
+      fprintf(fpOut, "\n");
+    } // k < finalNumBins
+    LogPrintfVerbatim ( LOG_DEBUG, "done.\n");
+
+    fclose(fpOut);
+
+  }
+//  Fseries = Fseries; //Else warning set but not used
+  for (k = 0; k < finalNumBins; ++k) {
+    UINT4 b = k * finalBinStep;
+//    REAL8 f0 = Freq0 + b * dFreq;
+    REAL8 psd = math_op(&(finalPSD->data[b]), finalBinSize, uvar.PSDmthopBins);
+    if ( isnan ( psd )) {
+        XLAL_ERROR ( EXIT_FAILURE, "Found Not-A-Number in psd[k=%d] = NAN ... exiting\n", k );
+    } else  {
+    (*Fseries)->data->data[k] = math_op(&(finalPSD->data[b]), finalBinSize, uvar.PSDmthopBins);
+    }
+  }
+
+  /* we are now done with the psd */
+  XLALDestroyMultiPSDVector  ( multiPSD);
+  XLALDestroyMultiSFTVector  ( inputSFTs);
+  XLALDestroyUINT4Vector     ( uint4vector);
+
+  XLALDestroyREAL8Vector ( overSFTs );
+  XLALDestroyREAL8Vector ( overIFOs );
+  XLALDestroyREAL8Vector ( finalPSD );
+  XLALDestroyREAL8Vector ( finalNormSFT );
+
+  return EXIT_SUCCESS;
+
+} /* main() */
+
+/** compute the various kinds of math. operation */
+REAL8 math_op(REAL8* data, size_t length, INT4 type) {
+
+  UINT4 i;
+  REAL8 res = 0.0;
+
+  switch (type) {
+
+  case MATH_OP_ARITHMETIC_SUM: /* sum(data) */
+
+    for (i = 0; i < length; ++i) res += *(data++);
+
+    break;
+
+  case MATH_OP_ARITHMETIC_MEAN: /* sum(data)/length  */
+
+    for (i = 0; i < length; ++i) res += *(data++);
+    res /= (REAL8)length;
+
+    break;
+
+  case MATH_OP_ARITHMETIC_MEDIAN: /* middle element of sort(data) */
+
+    gsl_sort(data, 1, length);
+    if (length/2 == (length+1)/2) /* length is even */ {
+      res = (data[length/2] + data[length/2+1])/2;
+    }
+    else /* length is odd */ {
+      res = data[length/2];
+    }
+
+    break;
+
+  case MATH_OP_HARMONIC_SUM: /* 1 / sum(1 / data) */
+
+    for (i = 0; i < length; ++i) res += 1.0 / *(data++);
+    res = 1.0 / res;
+
+    break;
+
+  case MATH_OP_HARMONIC_MEAN: /* length / sum(1 / data) */
+
+    for (i = 0; i < length; ++i) res += 1.0 / *(data++);
+    res = (REAL8)length / res;
+
+    break;
+
+  case MATH_OP_POWERMINUS2_SUM: /*   1 / sqrt ( sum(1 / data/data) )*/
+
+    for (i = 0; i < length; ++i) res += 1.0 / (data[i]*data[i]);
+    res = 1.0 / sqrt(res);
+
+    break;
+
+   case MATH_OP_POWERMINUS2_MEAN: /*   1 / sqrt ( sum(1/data/data) / length )*/
+
+    for (i = 0; i < length; ++i) res += 1.0 / (data[i]*data[i]);
+    res = 1.0 / sqrt(res / (REAL8)length);
+
+    break;
+
+  case MATH_OP_MINIMUM: /* first element of sort(data) */
+
+    gsl_sort(data, 1, length);
+    res = data[0];
+    break;
+
+  case MATH_OP_MAXIMUM: /* first element of sort(data) */
+
+    gsl_sort(data, 1, length);
+    res = data[length-1];
+    break;
+
+  default:
+
+    XLALPrintError("'%i' is not a valid math. operation", type);
+    XLAL_ERROR_REAL8(XLAL_EINVAL);
+
+  }
+
+  return res;
+
+}
+
+/**
+ * Dump complete multi-PSDVector over IFOs, timestamps and frequency-bins into
+ * per-IFO ASCII output-files 'outbname-IFO'
+ *
+ */
+int
+XLALDumpMultiPSDVector ( const CHAR *outbname,			/**< output basename 'outbname' */
+                         const MultiPSDVector *multiPSDVect	/**< multi-psd vector to output */
+                  )
+{
+  /* check input consistency */
+  if ( outbname == NULL ) {
+    XLALPrintError ("%s: NULL input 'outbname'\n", __func__ );
+    XLAL_ERROR ( XLAL_EINVAL );
+  }
+  if ( multiPSDVect == NULL ) {
+    XLALPrintError ("%s: NULL input 'multiPSDVect'\n", __func__ );
+    XLAL_ERROR ( XLAL_EINVAL );
+  }
+  if ( multiPSDVect->length == 0 || multiPSDVect->data==0 ) {
+    XLALPrintError ("%s: invalid multiPSDVect input (length=0 or data=NULL)\n", __func__ );
+    XLAL_ERROR ( XLAL_EINVAL );
+  }
+
+  CHAR *fname;
+  FILE *fp;
+
+  UINT4 len = strlen ( outbname ) + 4;
+  if ( ( fname = XLALMalloc ( len * sizeof(CHAR) )) == NULL ) {
+    XLALPrintError ("%s: XLALMalloc(%d) failed.\n", __func__, len );
+    XLAL_ERROR ( XLAL_ENOMEM);
+  }
+
+  UINT4 numIFOs = multiPSDVect->length;
+  UINT4 X;
+  for ( X = 0; X < numIFOs; X ++ )
+    {
+      PSDVector *thisPSDVect = multiPSDVect->data[X];
+      char buf[100];
+
+      sprintf ( fname, "%s-%c%c", outbname, thisPSDVect->data[0].name[0], thisPSDVect->data[0].name[1] );
+
+      if ( ( fp = fopen( fname, "wb" ))  == NULL ) {
+        XLALPrintError ("%s: Failed to open PSDperSFT file '%s' for writing!\n", __func__, fname );
+        XLALFree ( fname );
+        XLAL_ERROR ( XLAL_ESYS );
+      }
+
+      REAL8 f0       = thisPSDVect->data[0].f0;
+      REAL8 dFreq    = thisPSDVect->data[0].deltaF;
+      UINT4 numFreqs = thisPSDVect->data[0].data->length;
+      UINT4 iFreq;
+
+      /* write comment header line into this output file */
+      /* FIXME: output code-version/cmdline/history info */
+      fprintf(fp,"%%%% first line holds frequencies [Hz] of PSD-Columns\n");
+      /* loop over frequency and output comnment-header markers */
+      fprintf(fp,"%%%%%-17s", "dummy");
+      for ( iFreq = 0; iFreq < numFreqs; iFreq ++ )
+        {
+          sprintf (buf, "f%d [Hz]", iFreq + 1 );
+          fprintf (fp, " %-21s", buf );
+        }
+      fprintf (fp, "\n");
+
+      /* write parseable header-line giving bin frequencies for PSDs */
+      fprintf (fp, "%-19d", -1 );
+      for (iFreq = 0; iFreq < numFreqs; iFreq++ )
+        fprintf (fp, " %-21.16g", f0 + iFreq * dFreq );
+      fprintf (fp, "\n\n\n");
+
+      /* output another header line describing the following format "ti[GPS] PSD(f1) ... " */
+      fprintf(fp,"%%%%%-17s", "ti[GPS]");
+      for ( iFreq = 0; iFreq < numFreqs; iFreq ++ )
+        {
+          sprintf (buf, "PSD(f%d)", iFreq + 1 );
+          fprintf (fp, " %-21s", buf );
+        }
+      fprintf (fp, "\n");
+
+      /* loop over timestamps: dump all PSDs over frequencies into one line */
+      UINT4 numTS = thisPSDVect->length;
+      UINT4 iTS;
+      for ( iTS = 0; iTS < numTS; iTS++ )
+        {
+          REAL8FrequencySeries *thisPSD = &thisPSDVect->data[iTS];
+
+          /* first output timestamp GPS time for this line */
+          REAL8 tGPS = XLALGPSGetREAL8( &thisPSD->epoch );
+          fprintf (fp, "%-19.18g", tGPS );
+
+          /* some internal consistency/paranoia checks */
+          if ( ( f0 != thisPSD->f0) || ( dFreq != thisPSD->deltaF ) || (numFreqs != thisPSD->data->length ) ) {
+            XLALPrintError ("%s: %d-th timestamp %f: inconsistent PSDVector: f0 = %g : %g,  dFreq = %g : %g, numFreqs = %d : %d \n",
+                            __func__, iTS, tGPS, f0, thisPSD->f0, dFreq, thisPSD->deltaF, numFreqs, thisPSD->data->length );
+            XLALFree ( fname );
+            fclose ( fp );
+            XLAL_ERROR ( XLAL_EDOM );
+          }
+
+          /* loop over all frequencies and dump PSD-value */
+          for ( iFreq = 0; iFreq < numFreqs; iFreq ++ )
+            fprintf (fp, " %-21.16g", thisPSD->data->data[iFreq] );
+
+          fprintf (fp, "\n");
+
+        } /* for iTS < numTS */
+
+      fclose ( fp );
+
+    } /* for X < numIFOs */
+
+  XLALFree ( fname );
+
+  return XLAL_SUCCESS;
+
+} /* XLALDumpMultiPSDVector() */
+
+/**
+ * Function that *truncates the PSD in place* to the requested frequency-bin interval [firstBin, lastBin] for the given multiPSDVector.
+ * Now also truncates the original SFT vector, as necessary for correct computation of normalized SFT power.
+ */
+int
+XLALCropMultiPSDandSFTVectors ( MultiPSDVector *multiPSDVect,
+                         MultiSFTVector *multiSFTVect,
+                         UINT4 firstBin,
+                         UINT4 lastBin
+                         )
+{
+  /* check user input */
+  if ( !multiPSDVect ) {
+    XLALPrintError ("%s: invalid NULL input 'multiPSDVect'\n", __func__ );
+    XLAL_ERROR ( XLAL_EINVAL );
+  }
+  if ( !multiSFTVect ) {
+    XLALPrintError ("%s: invalid NULL input 'multiSFTVect'\n", __func__ );
+    XLAL_ERROR ( XLAL_EINVAL );
+  }
+  if ( lastBin < firstBin ) {
+    XLALPrintError ("%s: empty bin interval requested [%d, %d]\n", __func__, firstBin, lastBin );
+    XLAL_ERROR ( XLAL_EDOM );
+  }
+
+  UINT4 numIFOs = multiPSDVect->length;
+  UINT4 numBins = multiPSDVect->data[0]->data[0].data->length;
+
+  if ( numIFOs != multiSFTVect->length ) {
+    XLALPrintError ("%s: inconsistent number of IFOs between PSD (%d) and SFT (%d) vectors.\n", __func__, numIFOs, multiSFTVect->length );
+    XLAL_ERROR ( XLAL_EDOM );
+  }
+
+  if ( numBins != multiSFTVect->data[0]->data[0].data->length ) {
+    XLALPrintError ("%s: inconsistent number of bins between PSD (%d bins) and SFT (%d bins) vectors.\n", __func__, numBins, multiSFTVect->data[0]->data[0].data->length );
+    XLAL_ERROR ( XLAL_EDOM );
+  }
+
+  if ( (firstBin >= numBins) || (lastBin >= numBins ) ) {
+    XLALPrintError ("%s: requested bin-interval [%d, %d] outside of PSD bins [0, %d]\n", __func__, firstBin, lastBin, numBins - 1 );
+    XLAL_ERROR ( XLAL_EDOM );
+  }
+
+  /* ----- check if there's anything to do at all? ----- */
+  if ( (firstBin == 0)  && (lastBin == numBins - 1) )
+    return XLAL_SUCCESS;
+
+  /* ----- loop over detectors, timestamps, then crop each PSD ----- */
+  UINT4 X;
+  for ( X=0; X < numIFOs; X ++ )
+    {
+      PSDVector *thisPSDVect = multiPSDVect->data[X];
+      SFTVector *thisSFTVect = multiSFTVect->data[X];
+      UINT4 numTS   = thisPSDVect->length;
+
+      UINT4 iTS;
+      for ( iTS = 0; iTS < numTS; iTS ++ )
+        {
+          REAL8FrequencySeries *thisPSD = &thisPSDVect->data[iTS];
+          COMPLEX8FrequencySeries *thisSFT = &thisSFTVect->data[iTS];
+
+          if ( numBins != thisPSD->data->length ) {
+            XLALPrintError ("%s: inconsistent number of frequency-bins across multiPSDVector: X=%d, iTS=%d: numBins = %d != %d\n",
+                            __func__, X, iTS, numBins, thisPSD->data->length );
+            XLAL_ERROR ( XLAL_EDOM );
+          }
+
+          if ( numBins != thisSFT->data->length ) {
+            XLALPrintError ("%s: inconsistent number of frequency-bins across multiSFTVector: X=%d, iTS=%d: numBins = %d != %d\n",
+                            __func__, X, iTS, numBins, thisSFT->data->length );
+            XLAL_ERROR ( XLAL_EDOM );
+          }
+
+          UINT4 numNewBins = lastBin - firstBin + 1;
+
+          /* crop PSD */
+          XLAL_CHECK(XLALResizeREAL8FrequencySeries(thisPSD, firstBin, numNewBins) != NULL, XLAL_EFUNC);
+
+          /* crop SFT */
+          XLAL_CHECK(XLALResizeCOMPLEX8FrequencySeries(thisSFT, firstBin, numNewBins) != NULL, XLAL_EFUNC);
+
+        } /* for iTS < numTS */
+
+    } /* for X < numIFOs */
+
+  /* that should be all ... */
+  return XLAL_SUCCESS;
+
+} /* XLALCropMultiPSDandSFTVectors() */

lalapps/src/pulsar/MakeData/ComputePSDFunction.h

+/*
+*  Copyright (C) 2017 Benjamin Steltner
+*
+*  This program is free software; you can redistribute it and/or modify
+*  it under the terms of the GNU General Public License as published by
+*  the Free Software Foundation; either version 2 of the License, or
+*  (at your option) any later version.
+*
+*  This program is distributed in the hope that it will be useful,
+*  but WITHOUT ANY WARRANTY; without even the implied warranty of
+*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*  GNU General Public License for more details.
+*
+*  You should have received a copy of the GNU General Public License
+*  along with with program; see the file COPYING. If not, write to the
+*  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
+*  MA  02111-1307  USA
+*/
+
+#ifndef _COMPUTEPSDFUNCTION_H
+#define _COMPUTEPSDFUNCTION_H
+
+#include <lal/LALDatatypes.h>
+#include "PSDConfigVars.h"
+
+#if defined(__cplusplus)
+extern "C" {
+#elif 0
+} /* so that editors will match preceding brace */
+#endif
+
+/**
+  * \defgroup ComputePSDFunction_h Header ComputePSDFunction_h
+  * \ingroup lal_fft
+  *
+  * \brief Computes the PSD of SFT
+  *
+  * ### Synopsis ###
+  *
+  * \code
+  * #include <lal/ComputePSDFunction.h>
+  * \endcode
+  *
+  * Computes PSD
+  *
+  * ### Description ###
+**/
+
+int computePSD(PSDUserVarsTag PSDUserVars, SFTtype *sft, REAL8FrequencySeries **Fseries);
+
+#endif

lalapps/src/pulsar/MakeData/MakeSFTsFunction.h

+/*
+*  Copyright (C) 2017 Benjamin Steltner
+*
+*  This program is free software; you can redistribute it and/or modify
+*  it under the terms of the GNU General Public License as published by
+*  the Free Software Foundation; either version 2 of the License, or
+*  (at your option) any later version.
+*
+*  This program is distributed in the hope that it will be useful,
+*  but WITHOUT ANY WARRANTY; without even the implied warranty of
+*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*  GNU General Public License for more details.
+*
+*  You should have received a copy of the GNU General Public License
+*  along with with program; see the file COPYING. If not, write to the
+*  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
+*  MA  02111-1307  USA
+*/
+
+#ifndef _MAKESFTSFUNCTION_H
+#define _MAKESFTSFUNCTION_H
+
+#include <lal/LALDatatypes.h>
+#include <lal/ComplexFFT.h>
+#include <lal/RealFFT.h>
+#include <lal/Window.h>
+#include "SFTConfigVars.h"
+
+#if defined(__cplusplus)
+extern "C" {
+#elif 0
+} /* so that editors will match preceding brace */
+#endif
+
+/**
+  * \defgroup MakeSFTs_h Header MakeSFTs_h
+  * \ingroup lal_fft
+  *
+  * \brief Performs short fourier transforms
+  *
+  * ### Synopsis ###
+  *
+  * \code
+  * #include <lal/MakeSFTsFunction.h>
+  * \endcode
+  *
+  * Performs SFTs
+  *
+  * ### Description ###
+**/
+
+int makeSFT(SFTConfigVarsTag SFTConfigVars, REAL8TimeSeries *Tseries, REAL8 fmin, REAL8 band, SFTtype **sft);
+
+
+#endif /* _MAKESFTSFUNCTION_H */

lalapps/src/pulsar/MakeData/Makefile.am

@@ -15,11 +15,16 @@ EXTRA_PROGRAMS =
 
 if LALFRAME
 
-bin_PROGRAMS += lalapps_MakeSFTs
+bin_PROGRAMS += lalapps_MakeSFTs lalapps_gateStrain_v1
 
 lalapps_MakeSFTs_SOURCES = MakeSFTs.c
 lalapps_MakeSFTs_CPPFLAGS = $(AM_CPPFLAGS)
 
+lalapps_gateStrain_v1_SOURCES = gateStrain_v1.c MakeSFTsFunction.c ComputePSDFunction.c
+
+
+
+
 if PSS
 lalapps_MakeSFTs_SOURCES += XLALPSSInterface.c
 lalapps_MakeSFTs_CPPFLAGS += -DPSS_ENABLED

lalapps/src/pulsar/MakeData/PSDConfigVars.h

+#ifndef _PSDCONFIGVARS_H
+#define _PSDCONFIGVARS_H
+
+struct PSDUserVarsTag {
+  CHAR *inputData;      /**< directory for input sfts */
+  CHAR *outputPSD;      /**< directory for output sfts */
+  CHAR *outputSpectBname;
+
+  REAL8 Freq;   /**< *physical* start frequency to compute PSD for (excluding rngmed wings) */
+  REAL8 FreqBand; /**< *physical* frequency band to compute PSD for (excluding rngmed wings) */
+
+  REAL8 startTime;  /**< earliest SFT-timestamp to include */
+  REAL8 endTime;  /**< last SFT-timestamp to include */
+  CHAR *IFO;
+  CHAR  *timeStampsFile;
+  LALStringVector *linefiles;
+  INT4 blocksRngMed;  /**< number of running-median bins to use */
+  INT4 maxBinsClean;
+
+  INT4 PSDmthopSFTs;     /**< for PSD, type of math. operation over SFTs */
+  INT4 PSDmthopIFOs;     /**< for PSD, type of math. operation over IFOs */
+  BOOLEAN outputNormSFT; /**< output normalised SFT power? */
+  INT4 nSFTmthopSFTs;    /**< for norm. SFT, type of math. operation over SFTs */
+  INT4 nSFTmthopIFOs;    /**< for norm. SFT, type of math. operation over IFOs */
+
+  REAL8 binSizeHz;       /**< output PSD bin size in Hz */
+  INT4  binSize;         /**< output PSD bin size in no. of bins */
+  INT4  PSDmthopBins;    /**< for PSD, type of math. operation over bins */
+  INT4  nSFTmthopBins;   /**< for norm. SFT, type of math. operation over bins */
+  REAL8 binStepHz;       /**< output PSD bin step in Hz */
+  INT4  binStep;         /**< output PSD bin step in no. of bins */
+  BOOLEAN outFreqBinEnd; /**< output the end frequency of each bin? */
+
+  BOOLEAN dumpMultiPSDVector; /**< output multi-PSD vector over IFOs, timestamps, and frequencies into file(s) */
+  CHAR *outputQ;  /**< output the 'data-quality factor' Q(f) into this file */
+
+  REAL8 fStart;   /**< Start Frequency to load from SFT and compute PSD, including wings (it is RECOMMENDED to use --Freq instead) */
+  REAL8 fBand;    /**< Frequency Band to load from SFT and compute PSD, including wings (it is RECOMMENDED to use --FreqBand instead) */
+
+};
+
+typedef struct PSDUserVarsTag PSDUserVarsTag;
+
+#endif

lalapps/src/pulsar/MakeData/SFTConfigVars.h

+#ifndef _SFTCONFIGVARS_H
+#define _SFTCONFIGVARS_H
+
+struct SFTConfigVarsTag { /* Was "CommandLineArgsTag" */
+  REAL8 HPf;              /* High pass filtering frequency */
+  INT4 T;                 /* SFT duration */
+  char *stringT;          /* 12/27/05 gam; string with SFT duration */
+  INT4 GPSStart;
+  INT4 GPSEnd;
+  INT4 makeGPSDirs;        /* 12/27/05 gam; add option to make directories based on gps time */
+  INT4 sftVersion;         /* 12/28/05 gam; output SFT version */
+  char *commentField;      /* 12/28/05 gam; string comment for version 2 SFTs */
+  BOOLEAN htdata;          /* flag that indicates we're doing h(t) data */
+  BOOLEAN makeTmpFile;     /* 01/09/06 gam */
+  char *FrCacheFile;       /* Frame cache file */
+  char *ChannelName;
+  char *IFO;               /* 01/14/07 gam */
+  char *SFTpath;           /* path to SFT file location */
+  char *miscDesc;          /* 12/28/05 gam; string giving misc. part of the SFT description field in the filename */
+  INT4 PSSCleaning;    /* 1=YES and 0=NO*/
+  REAL8 PSSCleanHPf;       /* Cut frequency for the bilateral highpass filter. It has to be used only if PSSCleaning is YES.*/
+  INT4 PSSCleanExt;        /* Extend the timeseries at the beginning before calculating the autoregressive mean */
+  INT4 windowOption;       /* 12/28/05 gam; window options; 0 = no window, 1 = default = Matlab style Tukey window; 2 = make_sfts.c Tukey window; 3 = Hann window */
+  REAL8 windowR;
+  REAL8 overlapFraction;   /* 12/28/05 gam; overlap fraction (for use with windows; e.g., use -P 0.5 with -w 3 Hann windows; default is 1.0). */
+  BOOLEAN useSingle;       /* 11/19/05 gam; use single rather than double precision */
+  char *frameStructType;   /* 01/10/07 gam */
+
+  BOOLEAN writeFile;       /* Write file or not */
+}; /* Was CommandLineArgs */
+
+//extern struct SFTConfigVarsTag SFTConfigVars;
+
+struct SFTheadertag { /* Was "headertag" */
+  REAL8 endian;
+  INT4  gps_sec;
+  INT4  gps_nsec;
+  REAL8 tbase;
+  INT4  firstfreqindex;
+  INT4  nsamples;
+}; /* Was header */
+
+//extern struct SFTheadertag SFTheader;
+typedef struct SFTConfigVarsTag SFTConfigVarsTag;
+typedef struct SFTheadertag SFTheadertag;
+
+
+
+#endif /* _SFTCONFIGVARS_H */

lalapps/src/pulsar/MakeData/testMakeSFTs.sh

@@ -39,10 +39,122 @@ for file in $MSFTsft; do
 done
 
 ## compare SFTs produced by MFDv5 and MakeSFTs, should be very nearly identical
-tol=1e-10
+tol=1.5e-10
 cmdline="lalapps_compareSFTs -V -e $tol -1 $MFDv5sft -2 $MSFTsft"
 echo "Comparing SFTs produced by MFDv5 and MakeSFTs, allowed tolerance=$tol:"
 if ! eval "$cmdline"; then
     echo "ERROR: something failed when running '$cmdline'"
     exit 1
 fi
+
+#Make SFT analogue to gatestrain
+Band3=200
+highpass=0
+echo "Make SFT analogue to gatestrain"
+# high-pass-freq (-f) sft-duration (-t) sft-write-path (-p) frame-cache (-C) gps-start-time (-s) gps-end-time (-e) channel-name (-N)
+# window-type (-w) window-radius (-r) overlap-fraction (-P) frame-struct-type (-u) start-freq (-F) band (-B)
+cmdline="lalapps_MakeSFTs -f ${highpass} -t ${Tsft} -p . -C $framecache -s ${tstart} -e ${tend} -N H1:mfdv5 -v 2 -i H1 -u PROC_REAL8 -w 1 -r 0.001 -F ${highpass} -B ${Band3} -D 4 -X MSFTGate"
+if ! eval "$cmdline"; then
+    echo "ERROR: something failed when running '$cmdline'"
+    exit 1
+fi
+MSFTsft="./H-1_H1_${Tsft}SFT_MSFTGate-1257/H-1_H1_${Tsft}SFT_MSFTGate-${tstart}-${Tsft}.sft"
+
+../../../frametools/lalapps_frread H-H1_mfdv5-1257741529-1800.gwf H1:mfdv5
+mv H1:mfdv5-1257741529.dat Mfd.gwf.dat
+
+## run gatestrain
+echo "Run gateStrain_v1"
+cmdline="lalapps_gateStrain_v1 --startTime ${tstart} --duration ${Tsft} --SFTFreqStart ${highpass} --SFTFreqBand ${Band3} --outDir . --outputGWF 1 --inFrames $framecache --inFrChannel H1:mfdv5 --outFrChannels H1:mfdv5 --padData 0.0 --Highpass ${highpass}"
+if ! eval "$cmdline"; then
+    echo "ERROR: something failed when running '$cmdline'"
+    exit 1
+fi
+GSFTsft="H-1_H1_1800SFT_GATED_WIN001-1257741529-1800.sft"
+
+../../../frametools/lalapps_frread H-H1_gating-1257741529-1800.gwf H1:mfdv5
+mv H1:mfdv5-1257741529.dat Gate.gwf.dat
+
+#Compare gwf:
+echo "Compare gwf produced by MFD and gatestrain"
+cmdline="diff -q Mfd.gwf.dat Gate.gwf.dat"
+if ! eval "$cmdline"; then
+    echo "ERROR: something failed when running '$cmdline'"
+    exit 1
+fi
+
+cmdline="lalapps_compareSFTs -V -e $tol -1 $MSFTsft -2 $GSFTsft"
+echo "Comparing SFTs produced by MFDv5 and gateStrain_v1, allowed tolerance=$tol:"
+if ! eval "$cmdline"; then
+    echo "ERROR: something failed when running '$cmdline'"
+    exit 1
+fi
+
+rm ${GSFTsft}
+
+## run gatestrain in a few different configurations
+gateStart=`echo "${tstart} + ${Tsft} / 3" | bc`
+gateEnd=`echo "${gateStart} + 3" | bc`
+gateFile=gateFile.txt
+echo "${gateStart} ${gateEnd}" > ${gateFile}
+
+echo "Run gateStrain_v1"
+cmdline="lalapps_gateStrain_v1 --startTime ${tstart} --duration ${Tsft} --SFTFreqStart ${highpass} --SFTFreqBand ${Band3} --outDir . --outputGWF 0 --outputUngatedSFT=1 --outputpyCBCSFT=1 --inFrames $framecache --inFrChannel H1:mfdv5 --outFrChannels H1:mfdv5 --padData 0.0 --Highpass ${highpass} --externalGateFile ${gateFile} --additionalGateFile ${gateFile}"
+if ! eval "$cmdline"; then
+    echo "ERROR: something failed when running '$cmdline'"
+    exit 1
+fi
+#Different gatestrain output formats
+GSFTGatsft="H-1_H1_1800SFT_GATED_WIN001-1257741529-1800.sft"
+GSFTUngsft="H-1_H1_1800SFT_UNGATED_WIN001-1257741529-1800.sft"
+GSFTCBCsft="H-1_H1_1800SFT_pyCBCGATED_WIN001-1257741529-1800.sft"
+GSFTExtsft="H-1_H1_1800SFT_ExternalGateFile_WIN001-1257741529-1800.sft"
+
+cmdline="lalapps_compareSFTs -V -e $tol -1 ${GSFTGatsft} -2 ${GSFTExtsft}"
+echo "Comparing gated SFTs produced by additional gates and external gates (same file), allowed tolerance=$tol:"
+if ! eval "$cmdline"; then
+    echo "ERROR: something failed when running '$cmdline'"
+    exit 1
+fi
+
+cmdline="lalapps_compareSFTs -V -e $tol -1 ${GSFTUngsft} -2 ${GSFTCBCsft}"
+echo "Comparing ungated SFT to pyCBC gated SFT (with no gate found), allowed tolerance=$tol:"
+if ! eval "$cmdline"; then
+    echo "ERROR: something failed when running '$cmdline'"
+    exit 1
+fi
+
+#exit 1
+
+#Test with and without padding
+## run gatestrain
+echo "Run gateStrain_v1 with and without padding"
+pad=8
+padStart=`echo "${tstart} + 8" | bc`
+padDur=`echo "${Tsft} - 16" | bc`
+
+cmdline="lalapps_gateStrain_v1 --startTime ${padStart} --duration ${padDur} --SFTFreqStart ${highpass} --SFTFreqBand ${Band3} --outDir . --outputGWF 1 --inFrames $framecache --inFrChannel H1:mfdv5 --outFrChannels H1:mfdv5 --padData 0.0 --Highpass ${highpass}"
+if ! eval "$cmdline"; then
+    echo "ERROR: something failed when running '$cmdline'"
+    exit 1
+fi
+GSFTsft="H-1_H1_1784SFT_GATED_WIN001-1257741537-1784.sft"
+mv ${GSFTsft} ${GSFTsft}_NoPad
+
+## run gatestrain
+echo "Run gateStrain_v1"
+cmdline="lalapps_gateStrain_v1 --startTime ${tstart} --duration ${Tsft} --SFTFreqStart ${highpass} --SFTFreqBand ${Band3} --outDir . --outputGWF 1 --inFrames $framecache --inFrChannel H1:mfdv5 --outFrChannels H1:mfdv5 --padData ${pad} --Highpass ${highpass}"
+if ! eval "$cmdline"; then
+    echo "ERROR: something failed when running '$cmdline'"
+    exit 1
+fi
+GSFTsft="H-1_H1_1784SFT_GATED_WIN001-1257741537-1784.sft"
+mv ${GSFTsft} ${GSFTsft}_WithPad
+
+#Compare SFT with and without padding
+cmdline="lalapps_compareSFTs -V -e $tol -1 ${GSFTsft}_NoPad -2 ${GSFTsft}_WithPad"
+echo "Comparing ungated SFT to pyCBC gated SFT (with no gate found), allowed tolerance=$tol:"
+if ! eval "$cmdline"; then
+    echo "ERROR: something failed when running '$cmdline'"
+    exit 1
+fi