diff --git a/lalapps/src/pulsar/MakeData/ComputePSDFunction.c b/lalapps/src/pulsar/MakeData/ComputePSDFunction.c
new file mode 100644
index 0000000000000000000000000000000000000000..bd494f7ef76e10f75726bd1a605bc315bbe974db
--- /dev/null
+++ b/lalapps/src/pulsar/MakeData/ComputePSDFunction.c
@@ -0,0 +1,654 @@
+/*
+* 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() */
diff --git a/lalapps/src/pulsar/MakeData/ComputePSDFunction.h b/lalapps/src/pulsar/MakeData/ComputePSDFunction.h
new file mode 100644
index 0000000000000000000000000000000000000000..c4c4586af3d2906a0597fa6524e1eacb47f5d7d2
--- /dev/null
+++ b/lalapps/src/pulsar/MakeData/ComputePSDFunction.h
@@ -0,0 +1,51 @@
+/*
+* 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
diff --git a/lalapps/src/pulsar/MakeData/MakeSFTsFunction.c b/lalapps/src/pulsar/MakeData/MakeSFTsFunction.c
new file mode 100644
index 0000000000000000000000000000000000000000..cbe4ab153942aa3f236958b2551507d4835fd2c9
--- /dev/null
+++ b/lalapps/src/pulsar/MakeData/MakeSFTsFunction.c
@@ -0,0 +1,1251 @@
+/*
+* Copyright (C) 2007 Gregory Mendell
+* Copyright (C) 2010,2011,2016 Bernd Machenschalk
+* 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 Gregory Mendell, Xavier Siemens, Bruce Allen, Bernd Machenschalk, Benjamin Steltner
+ * \brief generate SFTs
+ */
+
+/**
+ * This file is a copy of MakeSFTs.c ( https://lscsoft.docs.ligo.org/lalsuite/lalapps/_make_s_f_ts_8c_source.html ).
+ * Most file I/O is removed, the main method takes an in-memory time series and computes & returns an in-memory SFT to the caller.
+ */
+
+#include <config.h>
+#if !defined HAVE_LIBGSL || !defined HAVE_LIBLALFRAME
+#include <stdio.h>
+int main(void) {fputs("disabled, no gsl or no lal frame library support.\n", stderr);return 1;}
+#else
+
+#ifdef HAVE_UNISTD_H
+#include <unistd.h>
+#endif
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include <glob.h>
+#include <errno.h>
+#include <stdarg.h>
+
+#include <lal/LALDatatypes.h>
+#include <lal/LALgetopt.h>
+#include <lal/LALStdlib.h>
+#include <lal/LALStdio.h>
+#include <lal/FileIO.h>
+#include <lal/AVFactories.h>
+#include <lal/LALCache.h>
+#include <lal/LALFrStream.h>
+#include <lal/Window.h>
+#include <lal/Calibration.h>
+#include <lal/LALConstants.h>
+#include <lal/TimeSeries.h>
+#include <lal/BandPassTimeSeries.h>
+#include <lal/LALStdlib.h>
+#include <lal/LALConstants.h>
+#include <lal/AVFactories.h>
+#include <lal/TimeFreqFFT.h>
+#include <lal/RealFFT.h>
+#include <lal/ComplexFFT.h>
+#include <lal/SFTfileIO.h>
+#include <lal/SFTutils.h>
+#include <lal/LALVCSInfo.h>
+#include <LALAppsVCSInfo.h>
+#include "MakeSFTsFunction.h"
+#include <lal/Sequence.h>
+
+/* track memory usage under linux */
+#define TRACKMEMUSE 0
+
+/* print the first NUMTOPRINT, middle NUMTOPRINT, and last NUMTOPRINT input/ouput data at various stages */
+#define PRINTEXAMPLEDATA 0
+#define NUMTOPRINT 2
+
+/* 06/26/07 gam; use finite to check that data does not contains a non-FINITE (+/- Inf, NaN) values */
+#define CHECKFORINFINITEANDNANS 1
+
+#define TESTSTATUS( pstat ) \
+ if ( (pstat)->statusCode ) { REPORTSTATUS(pstat); return 100; } else ((void)0)
+
+/***************************************************************************/
+
+/***************************************************************************/
+
+/* GLOBAL VARIABLES */
+REAL8 FMIN = 48.0; /* default start frequency */
+REAL8 DF = 2000.0; /* default band */
+
+REAL8 winFncRMS = 1.0; /* 10/05/12 gam; global variable with the RMS of the window function; default value is 1.0 */
+
+static LALStatus status;
+LALCache *framecache; /* frame reading variables */
+LALFrStream *framestream=NULL;
+
+REAL8TimeSeries resetdataDouble;
+REAL8TimeSeries dataDouble;
+REAL4TimeSeries dataSingle;
+INT2TimeSeries dataINT2;
+INT4TimeSeries dataINT4;
+INT8TimeSeries dataINT8;
+INT4 SegmentDuration;
+LIGOTimeGPS gpsepoch;
+
+/* REAL8Vector *window; */ /* 11/02/05 gam */
+
+REAL8FFTPlan *fftPlanDouble; /* fft plan and data container, double precision case */
+COMPLEX16Vector *fftDataDouble = NULL;
+
+REAL4FFTPlan *fftPlanSingle; /* 11/19/05 gam; fft plan and data container, single precision case */
+COMPLEX8Vector *fftDataSingle = NULL;
+
+CHAR allargs[16384]; /* 06/26/07 gam; copy all command line args into commentField, based on /lalapps/src/calibration/ComputeStrainDriver.c */
+/***************************************************************************/
+
+/* FUNCTION PROTOTYPES */
+/* Reads the command line */
+int ReadCommandLine(int argc,char *argv[],struct SFTConfigVarsTag *CLA);
+
+/* Allocates space for data time series */
+int AllocateDataNoFileIO(struct SFTConfigVarsTag CLA);
+
+/* High passes data */
+int DoHighPass(struct SFTConfigVarsTag CLA);
+
+/* Windows data */
+int WindowData(struct SFTConfigVarsTag CLA);
+int WindowDataTukey2(struct SFTConfigVarsTag CLA);
+int WindowDataHann(struct SFTConfigVarsTag CLA);
+
+#ifdef PSS_ENABLED
+/* Time Domain Cleaning with PSS functions */
+int PSSTDCleaningDouble(struct SFTConfigVarsTag CLA);
+int PSSTDCleaningREAL8(REAL8TimeSeries *LALTS, REAL4 highpassFrequency, INT4 extendTimeseries);
+#endif
+
+/* create an SFT */
+int CreateSFT(struct SFTConfigVarsTag CLA);
+
+/* writes out an SFT */
+int WriteSFT(struct SFTConfigVarsTag CLA);
+int WriteVersion2SFT(struct SFTConfigVarsTag CLA);
+int WriteVersion2SFTNoFileIO(struct SFTConfigVarsTag CLA, SFTtype **sft);
+
+/* Frees the memory */
+int FreeMem(struct SFTConfigVarsTag CLA);
+
+/* prototypes */
+FILE* tryopen(char *name, const char *mode);
+void getSFTDescField(CHAR *sftDescField, CHAR *numSFTs, CHAR *ifo, CHAR *stringT, CHAR *typeMisc);
+void mkSFTDir(CHAR *sftPath, CHAR *site, CHAR *numSFTs, CHAR *ifo, CHAR *stringT, CHAR *typeMisc,CHAR *gpstime, INT4 numGPSdigits);
+void mkSFTFilename(CHAR *sftFilename, CHAR *site, CHAR *numSFTs, CHAR *ifo, CHAR *stringT, CHAR *typeMisc,CHAR *gpstime);
+void mvFilenames(CHAR *filename1, CHAR *filename2);
+
+
+/* -------------------- function definitions -------------------- */
+
+
+/* This repeatedly tries to re-open a file. This is useful on a
+ cluster because automount may fail or time out. This allows a
+ number of tries with a bit of rest in between.
+*/
+FILE* tryopen(char *name, const char *mode){
+ int count=0;
+ FILE *fp;
+
+ while (!(fp=fopen(name, mode))){
+ fprintf(stderr,"Unable to open file %s in mode %s. Will retry...\n", name, mode);
+ if (count++<10)
+ sleep(10);
+ else
+ exit(3);
+ }
+ return fp;
+}
+
+
+/* 12/27/05 gam; function to create sftDescField for output directory or filename. */
+void getSFTDescField(CHAR *sftDescField, CHAR *numSFTs, CHAR *ifo, CHAR *stringT, CHAR *typeMisc) {
+ strcpy(sftDescField,numSFTs);
+ strcat(sftDescField, "_");
+ strcat(sftDescField,ifo);
+ strcat(sftDescField, "_");
+ strcat(sftDescField,stringT);
+ strcat(sftDescField, "SFT");
+ if (typeMisc != NULL) {
+ strcat(sftDescField, "_");
+ strcat(sftDescField, typeMisc);
+ }
+}
+
+/* 12/27/05 gam; concat to the sftPath the directory name based on GPS time; make this directory if it does not already exist */
+void mkSFTDir(CHAR *sftPath, CHAR *site, CHAR *numSFTs, CHAR *ifo, CHAR *stringT, CHAR *typeMisc,CHAR *gpstime, INT4 numGPSdigits) {
+ CHAR sftDescField[256];
+ CHAR mkdirCommand[256];
+ strcat(sftPath,"/");
+ strcat(sftPath,site);
+ strcat(sftPath,"-");
+ getSFTDescField(sftDescField, numSFTs, ifo, stringT, typeMisc);
+ strcat(sftPath,sftDescField);
+ strcat(sftPath,"-");
+ strncat(sftPath,gpstime,numGPSdigits);
+ sprintf(mkdirCommand,"mkdir -p %s",sftPath);
+ if ( system(mkdirCommand) ) XLALPrintError ("system() returned non-zero status\n");
+}
+
+/* 12/27/05 gam; make SFT file name according to LIGO T040164-01 specification */
+void mkSFTFilename(CHAR *sftFilename, CHAR *site, CHAR *numSFTs, CHAR *ifo, CHAR *stringT, CHAR *typeMisc,CHAR *gpstime) {
+ CHAR sftDescField[256];
+ strcpy(sftFilename,site);
+ strcat(sftFilename,"-");
+ getSFTDescField(sftDescField, numSFTs, ifo, stringT, typeMisc);
+ strcat(sftFilename,sftDescField);
+ strcat(sftFilename,"-");
+ strcat(sftFilename,gpstime);
+ strcat(sftFilename,"-");
+ strcat(sftFilename,stringT);
+ strcat(sftFilename,".sft");
+}
+
+/* 01/09/06 gam; move filename1 to filename2 */
+void mvFilenames(CHAR *filename1, CHAR *filename2) {
+ CHAR mvFilenamesCommand[512+4];
+ sprintf(mvFilenamesCommand,"mv %s %s",filename1,filename2);
+ if ( system(mvFilenamesCommand) ) XLALPrintError ("system() returned non-zero status\n");
+}
+
+#if TRACKMEMUSE
+void printmemuse() {
+ pid_t mypid=getpid();
+ char commandline[256];
+ fflush(NULL);
+ sprintf(commandline,"cat /proc/%d/status | /bin/grep Vm | /usr/bin/fmt -140 -u", (int)mypid);
+ if ( system(commandline) ) XLALPrintError ("system() returned non-zero status\n");
+ fflush(NULL);
+ }
+#endif
+
+#if PRINTEXAMPLEDATA
+void printExampleDataSingle() {
+ INT4 i,dataLength;
+ dataLength = dataSingle.data->length;
+ fprintf(stdout,"dataSingle.deltaT, 1.0/dataSingle.deltaT = %23.16e, %23.16e\n",dataSingle.deltaT,1.0/dataSingle.deltaT);
+ fprintf(stdout,"dataSingle.epoch.gpsSeconds,dataSingle.epoch.gpsNanoSeconds = %i, %i\n",dataSingle.epoch.gpsSeconds,dataSingle.epoch.gpsNanoSeconds);
+ for (i=0;i<NUMTOPRINT;i++) {
+ fprintf(stdout,"%i %23.16e\n",i,dataSingle.data->data[i]);
+ }
+ for (i=dataLength/2;i<dataLength/2+NUMTOPRINT;i++) {
+ fprintf(stdout,"%i %23.16e\n",i,dataSingle.data->data[i]);
+ }
+ for (i=dataLength-NUMTOPRINT;i<dataLength;i++) {
+ fprintf(stdout,"%i %23.16e\n",i,dataSingle.data->data[i]);
+ }
+ fflush(stdout);
+}
+
+void printExampleDataDouble() {
+ INT4 i,dataLength;
+ dataLength = dataDouble.data->length;
+ fprintf(stdout,"dataDouble.deltaT, 1.0/dataDouble.deltaT = %23.16e, %23.16e\n",dataDouble.deltaT,1.0/dataDouble.deltaT);
+ fprintf(stdout,"dataDouble.epoch.gpsSeconds,dataDouble.epoch.gpsNanoSeconds = %i, %i\n",dataDouble.epoch.gpsSeconds,dataDouble.epoch.gpsNanoSeconds);
+ for (i=0;i<NUMTOPRINT;i++) {
+ fprintf(stdout,"%i %23.16e\n",i,dataDouble.data->data[i]);
+ }
+ for (i=dataLength/2;i<dataLength/2+NUMTOPRINT;i++) {
+ fprintf(stdout,"%i %23.16e\n",i,dataDouble.data->data[i]);
+ }
+ for (i=dataLength-NUMTOPRINT;i<dataLength;i++) {
+ fprintf(stdout,"%i %23.16e\n",i,dataDouble.data->data[i]);
+ }
+ fflush(stdout);
+}
+
+void printExampleFFTData(struct SFTConfigVarsTag CLA)
+{
+ int firstbin=(INT4)(FMIN*CLA.T+0.5);
+ int k;
+ INT4 nsamples=(INT4)(DF*CLA.T+0.5);
+
+ if(CLA.useSingle) {
+ for (k=0; k<NUMTOPRINT; k++) {
+ fprintf(stdout,"%i %23.16e %23.16e\n",k,fftDataSingle->data[k+firstbin].re,fftDataSingle->data[k+firstbin].im);
+ }
+ for (k=nsamples/2; k<nsamples/2+NUMTOPRINT; k++) {
+ fprintf(stdout,"%i %23.16e %23.16e\n",k,fftDataSingle->data[k+firstbin].re,fftDataSingle->data[k+firstbin].im);
+ }
+ for (k=nsamples-NUMTOPRINT; k<nsamples; k++) {
+ fprintf(stdout,"%i %23.16e %23.16e\n",k,fftDataSingle->data[k+firstbin].re,fftDataSingle->data[k+firstbin].im);
+ }
+ } else {
+ for (k=0; k<NUMTOPRINT; k++) {
+ fprintf(stdout,"%i %23.16e %23.16e\n",k,fftDataDouble->data[k+firstbin].re,fftDataDouble->data[k+firstbin].im);
+ }
+ for (k=nsamples/2; k<nsamples/2+NUMTOPRINT; k++) {
+ fprintf(stdout,"%i %23.16e %23.16e\n",k,fftDataDouble->data[k+firstbin].re,fftDataDouble->data[k+firstbin].im);
+ }
+ for (k=nsamples-NUMTOPRINT; k<nsamples; k++) {
+ fprintf(stdout,"%i %23.16e %23.16e\n",k,fftDataDouble->data[k+firstbin].re,fftDataDouble->data[k+firstbin].im);
+ }
+ }
+ fflush(stdout);
+}
+
+void printExampleSFTDataGoingToFile(struct SFTConfigVarsTag CLA)
+{
+ int firstbin=(INT4)(FMIN*CLA.T+0.5);
+ int k;
+ INT4 nsamples=(INT4)(DF*CLA.T+0.5);
+ REAL4 rpw,ipw;
+
+ if(CLA.useSingle) {
+ for (k=0; k<NUMTOPRINT; k++) {
+ rpw=((REAL4)(((REAL8)DF)/(0.5*(REAL8)(1/dataSingle.deltaT))))
+ * fftDataSingle->data[k+firstbin].re;
+ ipw=((REAL4)(((REAL8)DF)/(0.5*(REAL8)(1/dataSingle.deltaT))))
+ * fftDataSingle->data[k+firstbin].im;
+ fprintf(stdout,"%i %23.16e %23.16e\n",k,rpw,ipw);
+ }
+ for (k=nsamples/2; k<nsamples/2+NUMTOPRINT; k++) {
+ rpw=((REAL4)(((REAL8)DF)/(0.5*(REAL8)(1/dataSingle.deltaT))))
+ * fftDataSingle->data[k+firstbin].re;
+ ipw=((REAL4)(((REAL8)DF)/(0.5*(REAL8)(1/dataSingle.deltaT))))
+ * fftDataSingle->data[k+firstbin].im;
+ fprintf(stdout,"%i %23.16e %23.16e\n",k,rpw,ipw);
+ }
+ for (k=nsamples-NUMTOPRINT; k<nsamples; k++) {
+ rpw=((REAL4)(((REAL8)DF)/(0.5*(REAL8)(1/dataSingle.deltaT))))
+ * fftDataSingle->data[k+firstbin].re;
+ ipw=((REAL4)(((REAL8)DF)/(0.5*(REAL8)(1/dataSingle.deltaT))))
+ * fftDataSingle->data[k+firstbin].im;
+ fprintf(stdout,"%i %23.16e %23.16e\n",k,rpw,ipw);
+ }
+ } else {
+ for (k=0; k<NUMTOPRINT; k++) {
+ rpw=(((REAL8)DF)/(0.5*(REAL8)(1/dataDouble.deltaT)))
+ * fftDataDouble->data[k+firstbin].re;
+ ipw=(((REAL8)DF)/(0.5*(REAL8)(1/dataDouble.deltaT)))
+ * fftDataDouble->data[k+firstbin].im;
+ fprintf(stdout,"%i %23.16e %23.16e\n",k,rpw,ipw);
+ }
+ for (k=nsamples/2; k<nsamples/2+NUMTOPRINT; k++) {
+ rpw=(((REAL8)DF)/(0.5*(REAL8)(1/dataDouble.deltaT)))
+ * fftDataDouble->data[k+firstbin].re;
+ ipw=(((REAL8)DF)/(0.5*(REAL8)(1/dataDouble.deltaT)))
+ * fftDataDouble->data[k+firstbin].im;
+ fprintf(stdout,"%i %23.16e %23.16e\n",k,rpw,ipw);
+ }
+ for (k=nsamples-NUMTOPRINT; k<nsamples; k++) {
+ rpw=(((REAL8)DF)/(0.5*(REAL8)(1/dataDouble.deltaT)))
+ * fftDataDouble->data[k+firstbin].re;
+ ipw=(((REAL8)DF)/(0.5*(REAL8)(1/dataDouble.deltaT)))
+ * fftDataDouble->data[k+firstbin].im;
+ fprintf(stdout,"%i %23.16e %23.16e\n",k,rpw,ipw);
+ }
+ }
+ fflush(stdout);
+}
+
+void printExampleVersion2SFTDataGoingToFile(struct SFTConfigVarsTag CLA, SFTtype *oneSFT)
+{
+ int k;
+ INT4 nsamples=(INT4)(DF*CLA.T+0.5);
+
+ for (k=0; k<NUMTOPRINT; k++) {
+ fprintf(stdout,"%i %23.16e %23.16e\n",k,oneSFT->data->data[k].re,oneSFT->data->data[k].im);
+ }
+ for (k=nsamples/2; k<nsamples/2+NUMTOPRINT; k++) {
+ fprintf(stdout,"%i %23.16e %23.16e\n",k,oneSFT->data->data[k].re,oneSFT->data->data[k].im);
+ }
+ for (k=nsamples-NUMTOPRINT; k<nsamples; k++) {
+ fprintf(stdout,"%i %23.16e %23.16e\n",k,oneSFT->data->data[k].re,oneSFT->data->data[k].im);
+ }
+}
+#endif
+
+
+#ifdef PSS_ENABLED
+/* debug function */
+int PrintREAL4ArrayToFile(const char*name, REAL4*array, UINT4 length);
+int PrintREAL4ArrayToFile(const char*name, REAL4*array, UINT4 length) {
+ UINT4 i;
+ FILE*fp=fopen(name,"w");
+ if(!fp) {
+ fprintf(stderr,"Could not open file '%s' for writing\n",name);
+ return -1;
+ } else {
+ for(i=0;i<length;i++)
+ fprintf(fp,"%23.16e\n",array[i]);
+ fclose(fp);
+ }
+ return 0;
+}
+#endif
+
+
+/************************************* MAIN PROGRAM *************************************/
+/**
+ * This method is a copy of the main method, but with all file reading and writing disabled for processing in memory
+ */
+int makeSFT(SFTConfigVarsTag NewSFTConfigVars, REAL8TimeSeries *Tseries, REAL8 fmin, REAL8 band, SFTtype **sft) {
+ SFTConfigVarsTag SFTConfigVars = NewSFTConfigVars;
+ FMIN = fmin;
+ DF = band;
+/** Reset everything for calling it multiple times **/
+
+ dataDouble = resetdataDouble;
+ fftPlanDouble = NULL;
+ fftDataDouble = NULL;
+
+ #if TRACKMEMUSE
+ printf("Memory use at startup is:\n"); printmemuse();
+ #endif
+
+ SegmentDuration = SFTConfigVars.GPSEnd - SFTConfigVars.GPSStart ;
+
+ if( SegmentDuration < SFTConfigVars.T)
+ {
+ fprintf(stderr,"Cannot fit an SFT of duration %d between %d and %d\n",
+ SFTConfigVars.T, SFTConfigVars.GPSStart, SFTConfigVars.GPSEnd );
+ return 0;;
+ }
+
+ gpsepoch.gpsSeconds = SFTConfigVars.GPSStart;
+ gpsepoch.gpsNanoSeconds = 0;
+
+ dataDouble = *Tseries;
+
+ /* Allocates space for data */
+ if (AllocateDataNoFileIO(SFTConfigVars)) return 2;
+
+
+ /* for(j=0; j < SegmentDuration/SFTConfigVars.T; j++) */ /* 12/28/05 gam */
+ while(gpsepoch.gpsSeconds + SFTConfigVars.T <= SFTConfigVars.GPSEnd)
+ {
+
+ /* gpsepoch.gpsSeconds = SFTConfigVars.GPSStart + j*SFTConfigVars.T;
+ gpsepoch.gpsNanoSeconds = 0; */ /* 12/28/05 gam; now updated at the bottom of while loop */
+
+ /* Reads T seconds of data */
+// if (ReadDataNoFileIO(SFTConfigVars)) return 3;
+
+ /* High-pass data with Butterworth filter */
+ if(DoHighPass(SFTConfigVars)) return 4;
+
+ /* Window data; 12/28/05 gam; add options */
+ if (SFTConfigVars.windowOption==1) {
+ if(WindowData(SFTConfigVars)) return 5; /* SFTConfigVars.windowOption==1 is the default */
+ } else if (SFTConfigVars.windowOption==2) {
+ if(WindowDataTukey2(SFTConfigVars)) return 5;
+ } else if (SFTConfigVars.windowOption==3) {
+ if(WindowDataHann(SFTConfigVars)) return 5;
+ } else {
+ /* Continue with no windowing; parsing of command line args makes sure options are one of the above or 0 for now windowing. */
+ }
+
+#ifdef PSS_ENABLED
+ /* Time Domain cleaning procedure */
+ if (SFTConfigVars.PSSCleaning)
+ if(PSSTDCleaningDouble(SFTConfigVars))
+ return 9;
+#endif
+
+ /* create an SFT */
+ if(CreateSFT(SFTConfigVars)) return 6;
+
+ /* write out sft - dummy write in this case */
+ if (SFTConfigVars.sftVersion==1) {
+ if(WriteSFT(SFTConfigVars)) return 7;
+ } else if (SFTConfigVars.sftVersion==2) {
+ if(WriteVersion2SFTNoFileIO(SFTConfigVars, sft)) return 7; /* default now is to output version 2 SFTs */
+ } else {
+ fprintf( stderr, "Invalid input for 'sft-version = %d', must be either '1' or '2'\n", SFTConfigVars.sftVersion );
+ return 0;;
+ }
+
+ gpsepoch.gpsSeconds = gpsepoch.gpsSeconds + (INT4)((1.0 - SFTConfigVars.overlapFraction)*((REAL8)SFTConfigVars.T));
+ gpsepoch.gpsNanoSeconds = 0;
+ }
+ XLALDestroyCOMPLEX16Vector(fftDataDouble);
+ XLALDestroyREAL8FFTPlan(fftPlanDouble);
+
+ #if TRACKMEMUSE
+ printf("Memory use after freeing all allocated memory:\n"); printmemuse();
+ #endif
+
+ return 0;
+} /* makeSFT(args) */
+
+
+/*** FUNCTIONS ***/
+
+/*******************************************************************************/
+/**
+ * This is a copy of AllocateData, but with no file reading and/or writing
+ */
+
+int AllocateDataNoFileIO(struct SFTConfigVarsTag CLA)
+{
+ dataSingle.deltaT=dataDouble.deltaT;
+
+ /* 11/19/05 gam; will keep dataDouble or dataSingle in memory at all times, depending on whether using single or double precision */
+ if(CLA.useSingle) {
+ LALCreateVector(&status,&dataSingle.data,(UINT4)(CLA.T/dataSingle.deltaT +0.5));
+ TESTSTATUS( &status );
+
+ #if TRACKMEMUSE
+ printf("Memory use after creating dataSingle and before calling LALCreateForwardRealFFTPlan.\n"); printmemuse();
+ #endif
+
+ fftPlanSingle = XLALCreateForwardREAL4FFTPlan( dataSingle.data->length, 0 );
+ XLAL_CHECK( fftPlanSingle != NULL, XLAL_EFUNC );
+
+
+ #if TRACKMEMUSE
+ printf("Memory use after creating dataSingle and after calling XLALCreateForwardREAL4FFTPlan.\n"); printmemuse();
+ #endif
+
+ } else {
+ // LALDCreateVector(&status,&dataDouble.data,(UINT4)(CLA.T/dataDouble.deltaT +0.5));
+ // TESTSTATUS( &status );
+
+ #if TRACKMEMUSE
+ printf("Memory use after creating dataDouble and before calling LALCreateForwardREAL8FFTPlan.\n"); printmemuse();
+ #endif
+
+ fftPlanDouble = XLALCreateForwardREAL8FFTPlan( dataDouble.data->length, 0 );
+ XLAL_CHECK( fftPlanDouble != NULL, XLAL_EFUNC );
+
+ #if TRACKMEMUSE
+ printf("Memory use after creating dataDouble and after calling XLALCreateForwardREAL8FFTPlan.\n"); printmemuse();
+ #endif
+ }
+
+ return 0;
+}
+/*******************************************************************************/
+
+/*******************************************************************************/
+int DoHighPass(struct SFTConfigVarsTag CLA)
+{
+ PassBandParamStruc filterpar;
+ char tmpname[] = "Butterworth High Pass";
+
+ filterpar.name = tmpname;
+ filterpar.nMax = 10;
+ filterpar.f2 = CLA.HPf;
+ filterpar.a2 = 0.5;
+ filterpar.f1 = -1.0;
+ filterpar.a1 = -1.0;
+
+ if (CLA.HPf > 0.0 )
+ {
+ if(CLA.useSingle) {
+ #if TRACKMEMUSE
+ printf("Memory use before calling LALDButterworthREAL4TimeSeries:\n"); printmemuse();
+ #endif
+
+ /* High pass the time series */
+ LALDButterworthREAL4TimeSeries(&status,&dataSingle,&filterpar);
+ TESTSTATUS( &status );
+
+ #if TRACKMEMUSE
+ printf("Memory use after calling LALDButterworthREAL4TimeSeries:\n"); printmemuse();
+ #endif
+
+ #if PRINTEXAMPLEDATA
+ printf("\nExample dataSingle values after filtering data in HighPass:\n"); printExampleDataSingle();
+ #endif
+ } else {
+ #if TRACKMEMUSE
+ printf("Memory use before calling LALButterworthREAL8TimeSeries:\n"); printmemuse();
+ #endif
+
+ /* High pass the time series */
+ LALButterworthREAL8TimeSeries(&status,&dataDouble,&filterpar);
+ TESTSTATUS( &status );
+
+ #if TRACKMEMUSE
+ printf("Memory use after calling LALButterworthREAL8TimeSeries:\n"); printmemuse();
+ #endif
+
+ #if PRINTEXAMPLEDATA
+ printf("\nExample dataDouble values after filtering data in HighPass:\n"); printExampleDataDouble();
+ #endif
+ }
+ }
+
+ return 0;
+}
+/*******************************************************************************/
+
+/*******************************************************************************/
+int WindowData(struct SFTConfigVarsTag CLA)
+{
+ REAL8 r = CLA.windowR;
+ INT4 k, N, kl, kh;
+ /* 10/05/12 gam */
+ REAL8 win;
+ /* This implementation of a Tukey window is describes
+ in the Matlab tukey window documentation */
+
+ /* initialize the RMS of the window function */
+ winFncRMS = 0.0;
+
+ /* 11/19/05 gam */
+ if(CLA.useSingle) {
+ N=dataSingle.data->length;
+ kl=r/2*(N-1)+1;
+ kh=N-r/2*(N-1)+1;
+ for(k = 1; k < kl; k++)
+ {
+ /* dataSingle.data->data[k-1]=dataSingle.data->data[k-1]*0.5*( (REAL4)( 1.0 + cos(LAL_TWOPI/r*(k-1)/(N-1) - LAL_PI) ) ); */
+ win = 0.5*( 1.0 + cos(LAL_TWOPI/r*(k-1)/(N-1) - LAL_PI) );
+ dataSingle.data->data[k-1] *= ((REAL4) win);
+ winFncRMS += win*win;
+ }
+ for(k = kh; k <= N; k++)
+ {
+ /*dataSingle.data->data[k-1]=dataSingle.data->data[k-1]*0.5*( (REAL4)( 1.0 + cos(LAL_TWOPI/r - LAL_TWOPI/r*(k-1)/(N-1) - LAL_PI) ) );*/
+ win = 0.5*( 1.0 + cos(LAL_TWOPI/r - LAL_TWOPI/r*(k-1)/(N-1) - LAL_PI) );
+ dataSingle.data->data[k-1] *= ((REAL4) win);
+ winFncRMS += win*win;
+ }
+
+ #if PRINTEXAMPLEDATA
+ printf("\nExample dataSingle values after windowing data in WindowData:\n"); printExampleDataSingle();
+ #endif
+ } else {
+ N=dataDouble.data->length;
+ kl=r/2*(N-1)+1;
+ kh=N-r/2*(N-1)+1;
+ for(k = 1; k < kl; k++)
+ {
+ /* dataDouble.data->data[k-1]=dataDouble.data->data[k-1]*0.5*( 1.0 + cos(LAL_TWOPI/r*(k-1)/(N-1) - LAL_PI) ); */
+ win = 0.5*( 1.0 + cos(LAL_TWOPI/r*(k-1)/(N-1) - LAL_PI) );
+ dataDouble.data->data[k-1] *= win;
+ winFncRMS += win*win;
+ }
+ for(k = kh; k <= N; k++)
+ {
+ /* dataDouble.data->data[k-1]=dataDouble.data->data[k-1]*0.5*( 1.0 + cos(LAL_TWOPI/r - LAL_TWOPI/r*(k-1)/(N-1) - LAL_PI) ); */
+ win = 0.5*( 1.0 + cos(LAL_TWOPI/r - LAL_TWOPI/r*(k-1)/(N-1) - LAL_PI) );
+ dataDouble.data->data[k-1] *= win;
+ winFncRMS += win*win;
+ }
+ #if PRINTEXAMPLEDATA
+ printf("\nExample dataDouble values after windowing data in WindowData:\n"); printExampleDataDouble();
+ #endif
+ }
+
+ /* Add to sum of squares of the window function the parts of window which are equal to 1, and then find RMS value*/
+ winFncRMS += (REAL8) (kh - kl);
+ winFncRMS = sqrt( ( winFncRMS/((REAL8) N) ) );
+
+ return 0;
+}
+/*******************************************************************************/
+
+/*******************************************************************************/
+/* Same as window function given in lalapps/src/pulsar/make_sfts.c */
+int WindowDataTukey2(struct SFTConfigVarsTag CLA)
+{
+ /* Define the parameters to make the window */
+ INT4 WINSTART = 4096;
+ INT4 WINEND = 8192;
+ INT4 WINLEN = (WINEND-WINSTART);
+ /* INT4 i; */
+ INT4 i, N;
+ REAL8 win;
+
+
+ /* initialize the RMS of the window function */
+ winFncRMS = 0.0;
+
+ if(CLA.useSingle) {
+ N=dataSingle.data->length;
+ /* window data. Off portion */
+ for (i=0; i<WINSTART; i++) {
+ dataSingle.data->data[i] = 0.0;
+ dataSingle.data->data[dataSingle.data->length - 1 - i] = 0.0;
+ }
+ /* window data, smooth turn-on portion */
+ for (i=WINSTART; i<WINEND; i++) {
+ win=((sin((i - WINSTART)*LAL_PI/(WINLEN)-LAL_PI_2)+1.0)/2.0);
+ dataSingle.data->data[i] *= win;
+ dataSingle.data->data[dataSingle.data->length - 1 - i] *= win;
+ winFncRMS += 2.0*win*win;
+ }
+ #if PRINTEXAMPLEDATA
+ printf("\nExample dataSingle values after windowing data in WindowDataTukey2:\n"); printExampleDataSingle();
+ #endif
+ } else {
+ N=dataDouble.data->length;
+ /* window data. Off portion */
+ for (i=0; i<WINSTART; i++) {
+ dataDouble.data->data[i] = 0.0;
+ dataDouble.data->data[dataDouble.data->length - 1 - i] = 0.0;
+ }
+ /* window data, smooth turn-on portion */
+ for (i=WINSTART; i<WINEND; i++) {
+ win=((sin((i - WINSTART)*LAL_PI/(WINLEN)-LAL_PI_2)+1.0)/2.0);
+ dataDouble.data->data[i] *= win;
+ dataDouble.data->data[dataDouble.data->length - 1 - i] *= win;
+ winFncRMS += 2.0*win*win;
+ }
+ #if PRINTEXAMPLEDATA
+ printf("\nExample dataDouble values after windowing data in WindowDataTukey2:\n"); printExampleDataDouble();
+ #endif
+ }
+
+ /* Add to sum of squares of the window function the parts of window which are equal to 1, and then find RMS value*/
+ winFncRMS += (REAL8) (N - 2*WINEND);
+ winFncRMS = sqrt( ( winFncRMS/((REAL8) N) ) );
+
+ return 0;
+}
+/*******************************************************************************/
+
+/*******************************************************************************/
+/* Hann window based on Matlab, but with C indexing: w[k] = 0.5*( 1 - cos(2*pi*k/(N-1)) ) k = 0, 1, 2,...N-1 */
+int WindowDataHann(struct SFTConfigVarsTag CLA)
+{
+ INT4 k;
+ REAL8 win,N,Nm1;
+ REAL8 real8TwoPi = 2.0*((REAL8)(LAL_PI));
+
+ /* initialize the RMS of the window function */
+ winFncRMS = 0.0;
+
+ if(CLA.useSingle) {
+ N = ((REAL8)dataSingle.data->length);
+ Nm1 = N - 1;
+ for (k=0; k<N; k++) {
+ win=0.5*( 1.0 - cos(real8TwoPi*((REAL8)(k))/Nm1) );
+ dataSingle.data->data[k] *= win;
+ winFncRMS += win*win;
+ }
+ #if PRINTEXAMPLEDATA
+ printf("\nExample dataSingle values after windowing data in WindowDataHann:\n"); printExampleDataSingle();
+ #endif
+ } else {
+ N = ((REAL8)dataDouble.data->length);
+ Nm1 = N - 1;
+ for (k=0; k<N; k++) {
+ win=0.5*( 1.0 - cos(real8TwoPi*((REAL8)(k))/Nm1) );
+ dataDouble.data->data[k] *= win;
+ winFncRMS += win*win;
+ }
+ #if PRINTEXAMPLEDATA
+ printf("\nExample dataDouble values after windowing data in WindowDataHann:\n"); printExampleDataDouble();
+ #endif
+ }
+
+ /* Find RMS value; note that N is REAL8 in this function */
+ winFncRMS = sqrt( (winFncRMS/N) );
+
+ return 0;
+}
+/*******************************************************************************/
+
+/*******************************************************************************/
+int CreateSFT(struct SFTConfigVarsTag CLA)
+{
+
+ /* 11/19/05 gam */
+ if(CLA.useSingle) {
+ #if TRACKMEMUSE
+ printf("Memory use before creating output vector fftDataSingle and calling XLALREAL4ForwardFFT:\n"); printmemuse();
+ #endif
+
+ /* 11/02/05 gam; allocate container for SFT data */
+ LALCCreateVector( &status, &fftDataSingle, dataSingle.data->length / 2 + 1 );
+ TESTSTATUS( &status );
+
+ /* compute sft */
+ XLAL_CHECK( XLALREAL4ForwardFFT( fftDataSingle, dataSingle.data, fftPlanSingle ) == XLAL_SUCCESS, XLAL_EFUNC );
+
+ #if TRACKMEMUSE
+ printf("Memory use after creating output vector fftDataSingle and calling XLALREAL4ForwardFFT:\n"); printmemuse();
+ #endif
+
+ #if PRINTEXAMPLEDATA
+ printf("\nExample real and imaginary value of fftDataSingle from CreateSFT:\n"); printExampleFFTData(CLA);
+ #endif
+ } else {
+ #if TRACKMEMUSE
+ printf("Memory use before creating output vector fftDataDouble and calling XLALREAL8ForwardFFT:\n"); printmemuse();
+ #endif
+
+ /* 11/02/05 gam; allocate container for SFT data */
+ LALZCreateVector( &status, &fftDataDouble, dataDouble.data->length / 2 + 1 );
+ TESTSTATUS( &status );
+
+ /* compute sft */
+ XLAL_CHECK( XLALREAL8ForwardFFT( fftDataDouble, dataDouble.data, fftPlanDouble ) == XLAL_SUCCESS, XLAL_EFUNC );
+
+ #if TRACKMEMUSE
+ printf("Memory use after creating output vector fftDataDouble and calling XLALREAL8ForwardFFT:\n"); printmemuse();
+ #endif
+
+ #if PRINTEXAMPLEDATA
+ printf("\nExample real and imaginary value of fftDataDouble from CreateSFT:\n"); printExampleFFTData(CLA);
+ #endif
+ }
+ return 0;
+}
+/*******************************************************************************/
+
+/*******************************************************************************/
+int WriteSFT(struct SFTConfigVarsTag CLA)
+{
+ char sftname[256];
+ char sftnameFinal[256]; /* 01/09/06 gam */
+ char numSFTs[2]; /* 12/27/05 gam */
+ char site[2]; /* 12/27/05 gam */
+ char ifo[3]; /* 12/27/05 gam; allow 3rd charactor for null termination */
+ int firstbin=(INT4)(FMIN*CLA.T+0.5), k;
+ FILE *fpsft;
+ int errorcode1=0,errorcode2=0;
+ char gpstime[11]; /* 12/27/05 gam; allow for 10 digit GPS times and null termination */
+ REAL4 rpw,ipw;
+
+ /* 12/27/05 gam; set up the number of SFTs, site, and ifo as null terminated strings */
+ numSFTs[0] = '1';
+ numSFTs[1] = '\0'; /* null terminate */
+ strncpy( site, CLA.ChannelName, 1 );
+ site[1] = '\0'; /* null terminate */
+ if (CLA.IFO != NULL) {
+ strncpy( ifo, CLA.IFO, 2 );
+ } else {
+ strncpy( ifo, CLA.ChannelName, 2 );
+ }
+ ifo[2] = '\0'; /* null terminate */
+ sprintf(gpstime,"%09d",gpsepoch.gpsSeconds);
+
+ strcpy( sftname, CLA.SFTpath );
+ /* 12/27/05 gam; add option to make directories based on gps time */
+ if (CLA.makeGPSDirs > 0) {
+ /* 12/27/05 gam; concat to the sftname the directory name based on GPS time; make this directory if it does not already exist */
+ mkSFTDir(sftname, site, numSFTs, ifo, CLA.stringT, CLA.miscDesc, gpstime, CLA.makeGPSDirs);
+ }
+
+ /* 01/09/06 gam; sftname will be temporary; will move to sftnameFinal. */
+ if(CLA.makeTmpFile) {
+ /* set up sftnameFinal with usual SFT name */
+ strcpy(sftnameFinal,sftname);
+ strcat(sftnameFinal,"/SFT_");
+ strcat(sftnameFinal,ifo);
+ strcat(sftnameFinal,".");
+ strcat(sftnameFinal,gpstime);
+ /* sftname begins with . and ends in .tmp */
+ strcat(sftname,"/.SFT_");
+ strcat(sftname,ifo);
+ strcat(sftname,".");
+ strcat(sftname,gpstime);
+ strcat(sftname,".tmp");
+ } else {
+ strcat(sftname,"/SFT_");
+ strcat(sftname,ifo);
+ strcat(sftname,".");
+ strcat(sftname,gpstime);
+ }
+
+ /* open SFT file for writing */
+ fpsft=tryopen(sftname,"w");
+
+ /* write header */
+ SFTheadertag header;
+ header.endian=1.0;
+ header.gps_sec=gpsepoch.gpsSeconds;
+ header.gps_nsec=gpsepoch.gpsNanoSeconds;
+ header.tbase=CLA.T;
+ header.firstfreqindex=firstbin;
+ header.nsamples=(INT4)(DF*CLA.T+0.5);
+ if (1!=fwrite((void*)&header,sizeof(header),1,fpsft)){
+ fprintf(stderr,"Error in writing header into file %s!\n",sftname);
+ return 7;
+ }
+
+ /* 11/19/05 gam */
+ if(CLA.useSingle) {
+ /* Write SFT */
+ for (k=0; k<header.nsamples; k++)
+ {
+ rpw=((REAL4)(((REAL8)DF)/(0.5*(REAL8)(1/dataSingle.deltaT))))
+ * crealf(fftDataSingle->data[k+firstbin]);
+ ipw=((REAL4)(((REAL8)DF)/(0.5*(REAL8)(1/dataSingle.deltaT))))
+ * cimagf(fftDataSingle->data[k+firstbin]);
+ /* 06/26/07 gam; use finite to check that data does not contains a non-FINITE (+/- Inf, NaN) values */
+ #if CHECKFORINFINITEANDNANS
+ if (!isfinite(rpw) || !isfinite(ipw)) {
+ fprintf(stderr, "Infinite or NaN data at freq bin %d.\n", k);
+ return 7;
+ }
+ #endif
+ errorcode1=fwrite((void*)&rpw, sizeof(REAL4),1,fpsft);
+ errorcode2=fwrite((void*)&ipw, sizeof(REAL4),1,fpsft);
+ }
+ #if PRINTEXAMPLEDATA
+ printf("\nExample real and imaginary SFT values going to file from fftDataSingle in WriteSFT:\n"); printExampleSFTDataGoingToFile(CLA);
+ #endif
+ LALCDestroyVector( &status, &fftDataSingle );
+ TESTSTATUS( &status );
+ #if TRACKMEMUSE
+ printf("Memory use after destroying fftDataSingle in WriteSFT:\n"); printmemuse();
+ #endif
+ } else {
+ /* Write SFT */
+ for (k=0; k<header.nsamples; k++)
+ {
+ rpw=(((REAL8)DF)/(0.5*(REAL8)(1/dataDouble.deltaT)))
+ * creal(fftDataDouble->data[k+firstbin]);
+ ipw=(((REAL8)DF)/(0.5*(REAL8)(1/dataDouble.deltaT)))
+ * cimag(fftDataDouble->data[k+firstbin]);
+ /* 06/26/07 gam; use finite to check that data does not contains a non-FINITE (+/- Inf, NaN) values */
+ #if CHECKFORINFINITEANDNANS
+ if (!isfinite(rpw) || !isfinite(ipw)) {
+ fprintf(stderr, "Infinite or NaN data at freq bin %d.\n", k);
+ return 7;
+ }
+ #endif
+ errorcode1=fwrite((void*)&rpw, sizeof(REAL4),1,fpsft);
+ errorcode2=fwrite((void*)&ipw, sizeof(REAL4),1,fpsft);
+ }
+ #if PRINTEXAMPLEDATA
+ printf("\nExample real and imaginary SFT values going to file from fftDataDouble in WriteSFT:\n"); printExampleSFTDataGoingToFile(CLA);
+ #endif
+ LALZDestroyVector( &status, &fftDataDouble );
+ TESTSTATUS( &status );
+ #if TRACKMEMUSE
+ printf("Memory use after destroying fftDataDouble in WriteSFT:\n"); printmemuse();
+ #endif
+ }
+
+ /* Check that there were no errors while writing SFTS */
+ if (errorcode1-1 || errorcode2-1)
+ {
+ fprintf(stderr,"Error in writing data into SFT file %s!\n",sftname);
+ return 8;
+ }
+
+ fclose(fpsft);
+
+ /* 01/09/06 gam; sftname is temporary; move to sftnameFinal. */
+ if(CLA.makeTmpFile) {
+ mvFilenames(sftname,sftnameFinal);
+ }
+
+ return 0;
+}
+/*******************************************************************************/
+
+/*******************************************************************************/
+/* 12/28/05 gam; write out version 2 SFT */
+int WriteVersion2SFT(struct SFTConfigVarsTag CLA)
+{
+ char sftname[256];
+ char sftFilename[256];
+ char sftnameFinal[256]; /* 01/09/06 gam */
+ char numSFTs[2]; /* 12/27/05 gam */
+ char site[2]; /* 12/27/05 gam */
+ char ifo[3]; /* 12/27/05 gam; allow 3rd charactor for null termination */
+ int firstbin=(INT4)(FMIN*CLA.T+0.5), k;
+ char gpstime[11]; /* 12/27/05 gam; allow for 10 digit GPS times and null termination */
+ SFTtype *oneSFT = NULL;
+ INT4 nBins = (INT4)(DF*CLA.T+0.5);
+ REAL4 singleDeltaT = 0.0; /* 01/05/06 gam */
+ REAL8 doubleDeltaT = 0.0; /* 01/05/06 gam */
+
+ /* 12/27/05 gam; set up the number of SFTs, site, and ifo as null terminated strings */
+ numSFTs[0] = '1';
+ numSFTs[1] = '\0'; /* null terminate */
+ strncpy( site, CLA.ChannelName, 1 );
+ site[1] = '\0'; /* null terminate */
+ if (CLA.IFO != NULL) {
+ strncpy( ifo, CLA.IFO, 2 );
+ } else {
+ strncpy( ifo, CLA.ChannelName, 2 );
+ }
+ ifo[2] = '\0'; /* null terminate */
+ sprintf(gpstime,"%09d",gpsepoch.gpsSeconds);
+
+ strcpy( sftname, CLA.SFTpath );
+ /* 12/27/05 gam; add option to make directories based on gps time */
+ if (CLA.makeGPSDirs > 0) {
+ /* 12/27/05 gam; concat to the sftname the directory name based on GPS time; make this directory if it does not already exist */
+ mkSFTDir(sftname, site, numSFTs, ifo, CLA.stringT, CLA.miscDesc, gpstime, CLA.makeGPSDirs);
+ }
+
+ strcat(sftname,"/");
+ mkSFTFilename(sftFilename, site, numSFTs, ifo, CLA.stringT, CLA.miscDesc, gpstime);
+ /* 01/09/06 gam; sftname will be temporary; will move to sftnameFinal. */
+ if(CLA.makeTmpFile) {
+ /* set up sftnameFinal with usual SFT name */
+ strcpy(sftnameFinal,sftname);
+ strcat(sftnameFinal,sftFilename);
+ /* sftname begins with . and ends in .tmp */
+ strcat(sftname,".");
+ strcat(sftname,sftFilename);
+ strcat(sftname,".tmp");
+ } else {
+ strcat(sftname,sftFilename);
+ }
+
+ /* make container to store the SFT data */
+ XLAL_CHECK( ( oneSFT = XLALCreateSFT ( ((UINT4)nBins)) ) != NULL, XLAL_EFUNC );
+ #if TRACKMEMUSE
+ printf("Memory use after creating oneSFT and calling XLALCreateSFT in WriteVersion2SFT:\n"); printmemuse();
+ #endif
+
+ /* copy the data to oneSFT */
+ strcpy(oneSFT->name,ifo);
+ oneSFT->epoch.gpsSeconds=gpsepoch.gpsSeconds;
+ oneSFT->epoch.gpsNanoSeconds=gpsepoch.gpsNanoSeconds;
+ oneSFT->f0 = FMIN;
+ oneSFT->deltaF = 1.0/((REAL8)CLA.T);
+ oneSFT->data->length=nBins;
+
+ if(CLA.useSingle) {
+ /* singleDeltaT = ((REAL4)dataSingle.deltaT); */ /* 01/05/06 gam; and normalize SFTs using this below */
+ singleDeltaT = (REAL4)(dataSingle.deltaT/winFncRMS); /* include 1 over window function RMS */
+ for (k=0; k<nBins; k++)
+ {
+ oneSFT->data->data[k] = (((REAL4) singleDeltaT) * fftDataSingle->data[k+firstbin]);
+ /* 06/26/07 gam; use finite to check that data does not contains a non-FINITE (+/- Inf, NaN) values */
+ #if CHECKFORINFINITEANDNANS
+ if (!isfinite(crealf(oneSFT->data->data[k])) || !isfinite(cimagf(oneSFT->data->data[k]))) {
+ fprintf(stderr, "Infinite or NaN data at freq bin %d.\n", k);
+ return 7;
+ }
+ #endif
+ }
+ #if PRINTEXAMPLEDATA
+ printf("\nExample real and imaginary SFT values going to file from fftDataSingle in WriteVersion2SFT:\n"); printExampleVersion2SFTDataGoingToFile(CLA,oneSFT);
+ #endif
+ LALCDestroyVector( &status, &fftDataSingle );
+ TESTSTATUS( &status );
+ #if TRACKMEMUSE
+ printf("Memory use after destroying fftDataSingle in WriteVersion2SFT:\n"); printmemuse();
+ #endif
+ } else {
+ /*doubleDeltaT = ((REAL8)dataDouble.deltaT); */ /* 01/05/06 gam; and normalize SFTs using this below */
+ doubleDeltaT = (REAL8)(dataDouble.deltaT/winFncRMS); /* include 1 over window function RMS */
+ for (k=0; k<nBins; k++)
+ {
+ oneSFT->data->data[k] = crectf( doubleDeltaT*creal(fftDataDouble->data[k+firstbin]), doubleDeltaT*cimag(fftDataDouble->data[k+firstbin]) );
+ /* 06/26/07 gam; use finite to check that data does not contains a non-FINITE (+/- Inf, NaN) values */
+ #if CHECKFORINFINITEANDNANS
+ if (!isfinite(crealf(oneSFT->data->data[k])) || !isfinite(cimagf(oneSFT->data->data[k]))) {
+ fprintf(stderr, "Infinite or NaN data at freq bin %d.\n", k);
+ return 7;
+ }
+ #endif
+ }
+ #if PRINTEXAMPLEDATA
+ printf("\nExample real and imaginary SFT values going to file from fftDataDouble in WriteVersion2SFT:\n"); printExampleVersion2SFTDataGoingToFile(CLA,oneSFT);
+ #endif
+ LALZDestroyVector( &status, &fftDataDouble );
+ TESTSTATUS( &status );
+ #if TRACKMEMUSE
+ printf("Memory use after destroying fftDataDouble in WriteVersion2SFT:\n"); printmemuse();
+ #endif
+ }
+
+ /* write the SFT */
+ XLAL_CHECK( XLALWriteSFT2file(oneSFT, sftname, CLA.commentField) == XLAL_SUCCESS, XLAL_EFUNC );
+
+ /* 01/09/06 gam; sftname is temporary; move to sftnameFinal. */
+ if(CLA.makeTmpFile) {
+ mvFilenames(sftname,sftnameFinal);
+ }
+
+ XLALDestroySFT (oneSFT);
+ #if TRACKMEMUSE
+ printf("Memory use after destroying oneSFT and calling XLALDestroySFT in WriteVersion2SFT:\n"); printmemuse();
+ #endif
+
+ return 0;
+}
+/*******************************************************************************/
+
+/*******************************************************************************/
+/**
+ * This is a copy of WriteVersion2 method, but with no File Reading / Writing,
+ * instead saving the sft into the pointer parameter sft
+ */
+int WriteVersion2SFTNoFileIO(struct SFTConfigVarsTag CLA, SFTtype **sft)
+{
+ char sftname[256];
+ char sftFilename[256];
+ char sftnameFinal[256]; /* 01/09/06 gam */
+ char numSFTs[2]; /* 12/27/05 gam */
+ char site[2]; /* 12/27/05 gam */
+ char ifo[3]; /* 12/27/05 gam; allow 3rd charactor for null termination */
+ int firstbin=(INT4)(FMIN*CLA.T+0.5), k;
+ char gpstime[11]; /* 12/27/05 gam; allow for 10 digit GPS times and null termination */
+ SFTtype *oneSFT = NULL;
+ INT4 nBins = (INT4)(DF*CLA.T+0.5);
+ REAL4 singleDeltaT = 0.0; /* 01/05/06 gam */
+ REAL8 doubleDeltaT = 0.0; /* 01/05/06 gam */
+
+ /* 12/27/05 gam; set up the number of SFTs, site, and ifo as null terminated strings */
+ numSFTs[0] = '1';
+ numSFTs[1] = '\0'; /* null terminate */
+ strncpy( site, CLA.ChannelName, 1 );
+ site[1] = '\0'; /* null terminate */
+ if (CLA.IFO != NULL) {
+ strncpy( ifo, CLA.IFO, 2 );
+ } else {
+ strncpy( ifo, CLA.ChannelName, 2 );
+ }
+ ifo[2] = '\0'; /* null terminate */
+ sprintf(gpstime,"%09d",gpsepoch.gpsSeconds);
+
+ strcpy( sftname, CLA.SFTpath );
+ /* 12/27/05 gam; add option to make directories based on gps time */
+ if (CLA.makeGPSDirs > 0) {
+ /* 12/27/05 gam; concat to the sftname the directory name based on GPS time; make this directory if it does not already exist */
+ mkSFTDir(sftname, site, numSFTs, ifo, CLA.stringT, CLA.miscDesc, gpstime, CLA.makeGPSDirs);
+ }
+
+ strcat(sftname,"/");
+ mkSFTFilename(sftFilename, site, numSFTs, ifo, CLA.stringT, CLA.miscDesc, gpstime);
+ /* 01/09/06 gam; sftname will be temporary; will move to sftnameFinal. */
+ if(CLA.makeTmpFile) {
+ /* set up sftnameFinal with usual SFT name */
+ strcpy(sftnameFinal,sftname);
+ strcat(sftnameFinal,sftFilename);
+ /* sftname begins with . and ends in .tmp */
+ strcat(sftname,".");
+ strcat(sftname,sftFilename);
+ strcat(sftname,".tmp");
+ } else {
+ strcat(sftname,sftFilename);
+ }
+
+ /* make container to store the SFT data */
+ XLAL_CHECK( ( oneSFT = XLALCreateSFT ( ((UINT4)nBins)) ) != NULL, XLAL_EFUNC );
+ #if TRACKMEMUSE
+ printf("Memory use after creating oneSFT and calling XLALCreateSFT in WriteVersion2SFT:\n"); printmemuse();
+ #endif
+
+ /* copy the data to oneSFT */
+ strcpy(oneSFT->name,ifo);
+ oneSFT->epoch.gpsSeconds=gpsepoch.gpsSeconds;
+ oneSFT->epoch.gpsNanoSeconds=gpsepoch.gpsNanoSeconds;
+ oneSFT->f0 = FMIN;
+ oneSFT->deltaF = 1.0/((REAL8)CLA.T);
+ oneSFT->data->length=nBins;
+
+ if(CLA.useSingle) {
+ /* singleDeltaT = ((REAL4)dataSingle.deltaT); */ /* 01/05/06 gam; and normalize SFTs using this below */
+ singleDeltaT = (REAL4)(dataSingle.deltaT/winFncRMS); /* include 1 over window function RMS */
+ for (k=0; k<nBins; k++)
+ {
+ oneSFT->data->data[k] = (((REAL4) singleDeltaT) * fftDataSingle->data[k+firstbin]);
+ /* 06/26/07 gam; use finite to check that data does not contains a non-FINITE (+/- Inf, NaN) values */
+ #if CHECKFORINFINITEANDNANS
+ if (!isfinite(crealf(oneSFT->data->data[k])) || !isfinite(cimagf(oneSFT->data->data[k]))) {
+ fprintf(stderr, "Infinite or NaN data at freq bin %d.\n", k);
+ return 7;
+ }
+ #endif
+ }
+ #if PRINTEXAMPLEDATA
+ printf("\nExample real and imaginary SFT values going to file from fftDataSingle in WriteVersion2SFT:\n"); printExampleVersion2SFTDataGoingToFile(CLA,oneSFT);
+ #endif
+ LALCDestroyVector( &status, &fftDataSingle );
+ TESTSTATUS( &status );
+ #if TRACKMEMUSE
+ printf("Memory use after destroying fftDataSingle in WriteVersion2SFT:\n"); printmemuse();
+ #endif
+ } else {
+ /*doubleDeltaT = ((REAL8)dataDouble.deltaT); */ /* 01/05/06 gam; and normalize SFTs using this below */
+ doubleDeltaT = (REAL8)(dataDouble.deltaT/winFncRMS); /* include 1 over window function RMS */
+ for (k=0; k<nBins; k++)
+ {
+ oneSFT->data->data[k] = crectf( doubleDeltaT*creal(fftDataDouble->data[k+firstbin]), doubleDeltaT*cimag(fftDataDouble->data[k+firstbin]) );
+ /* 06/26/07 gam; use finite to check that data does not contains a non-FINITE (+/- Inf, NaN) values */
+ #if CHECKFORINFINITEANDNANS
+ if (!isfinite(crealf(oneSFT->data->data[k])) || !isfinite(cimagf(oneSFT->data->data[k]))) {
+ fprintf(stderr, "Infinite or NaN data at freq bin %d.\n", k);
+ return 7;
+ }
+ #endif
+ }
+ #if PRINTEXAMPLEDATA
+ printf("\nExample real and imaginary SFT values going to file from fftDataDouble in WriteVersion2SFT:\n"); printExampleVersion2SFTDataGoingToFile(CLA,oneSFT);
+ #endif
+ LALZDestroyVector( &status, &fftDataDouble );
+ TESTSTATUS( &status );
+ #if TRACKMEMUSE
+ printf("Memory use after destroying fftDataDouble in WriteVersion2SFT:\n"); printmemuse();
+ #endif
+
+ }
+
+ /* write the SFT if SFTVars says so */
+ if(CLA.writeFile == 1) {
+ XLAL_CHECK( XLALWriteSFT2file(oneSFT, sftname, CLA.commentField) == XLAL_SUCCESS, XLAL_EFUNC );
+ } else {
+ printf("Omitting writing SFT to file, user requested no file IO in MakeSFTsFunction");
+ }
+
+ /* 01/09/06 gam; sftname is temporary; move to sftnameFinal. */
+ if(CLA.makeTmpFile) {
+ mvFilenames(sftname,sftnameFinal);
+ }
+
+ *sft = oneSFT;
+ return 0;
+}
+/*******************************************************************************/
+
+/*******************************************************************************/
+int FreeMem(struct SFTConfigVarsTag CLA)
+{
+
+ LALFrClose(&status,&framestream);
+ TESTSTATUS( &status );
+
+ /* 11/19/05 gam */
+ if(CLA.useSingle) {
+ LALDestroyVector(&status,&dataSingle.data);
+ TESTSTATUS( &status );
+ XLALDestroyREAL4FFTPlan( fftPlanSingle );
+ } else {
+ LALDDestroyVector(&status,&dataDouble.data);
+ TESTSTATUS( &status );
+ XLALDestroyREAL8FFTPlan( fftPlanDouble );
+ }
+
+ LALCheckMemoryLeaks();
+
+ return 0;
+}
+/*******************************************************************************/
+
+#endif
diff --git a/lalapps/src/pulsar/MakeData/MakeSFTsFunction.h b/lalapps/src/pulsar/MakeData/MakeSFTsFunction.h
new file mode 100644
index 0000000000000000000000000000000000000000..c0c9c81db9c4bb9a6cb01e3e0adf82c30850cc26
--- /dev/null
+++ b/lalapps/src/pulsar/MakeData/MakeSFTsFunction.h
@@ -0,0 +1,55 @@
+/*
+* 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 */
diff --git a/lalapps/src/pulsar/MakeData/Makefile.am b/lalapps/src/pulsar/MakeData/Makefile.am
index bd03f30b3309973b550fbc1d2229595bc7d92ab7..eb336f2c3b126dca7c6920af66994901791b4e5d 100644
--- a/lalapps/src/pulsar/MakeData/Makefile.am
+++ b/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
diff --git a/lalapps/src/pulsar/MakeData/PSDConfigVars.h b/lalapps/src/pulsar/MakeData/PSDConfigVars.h
new file mode 100644
index 0000000000000000000000000000000000000000..a4c04946acc6e2ee54c9d1b6f43a64c26a40d15d
--- /dev/null
+++ b/lalapps/src/pulsar/MakeData/PSDConfigVars.h
@@ -0,0 +1,44 @@
+#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
diff --git a/lalapps/src/pulsar/MakeData/SFTConfigVars.h b/lalapps/src/pulsar/MakeData/SFTConfigVars.h
new file mode 100644
index 0000000000000000000000000000000000000000..9a3fe6fbd4965110f1bf139ee472581f8a70809a
--- /dev/null
+++ b/lalapps/src/pulsar/MakeData/SFTConfigVars.h
@@ -0,0 +1,49 @@
+#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 */
diff --git a/lalapps/src/pulsar/MakeData/gateStrain_v1.c b/lalapps/src/pulsar/MakeData/gateStrain_v1.c
new file mode 100644
index 0000000000000000000000000000000000000000..39a72d844a3f9eafad370fe73166ba368d3d9faf
--- /dev/null
+++ b/lalapps/src/pulsar/MakeData/gateStrain_v1.c
@@ -0,0 +1,2163 @@
+/*
+* 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 Benjamin Steltner, Maria Alessandra Papa, Heinz-Bernd Eggenstein
+ * \brief Remove short-lived high-powered non-Gaussian noise transients from files. Generate 'gated' Short Fourier Transforms (SFTs).
+ */
+
+/*-----------------------------------------------------------------------
+ *
+ * File Name: gateStrain_v1.c
+ * Purpose: Remove short-lived high-powered non-Gaussian noise transients from files. Generate 'gated' Short Fourier Transforms (SFTs)
+ * Authors: Benjamin Steltner, Maria Alessandra Papa, Heinz-Bernd Eggenstein
+ * based on pycbc.strain.strain.detect_loud_glitches of the pyCBC python package by Alex Nitz (2013)
+ *
+ *-----------------------------------------------------------------------
+ */
+
+/* ---------- includes ---------- */
+
+#include "config.h"
+
+#include <sys/stat.h>
+
+#include <lal/FrequencySeries.h>
+#include <lal/TimeSeries.h>
+#include <lal/ResampleTimeSeries.h>
+#include <lal/UserInput.h>
+#include <lal/SFTfileIO.h>
+#include <lal/BandPassTimeSeries.h>
+
+#include <lal/Random.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_interp.h>
+
+#include <lal/Window.h>
+#include <lal/LALString.h>
+#include <lal/LALCache.h>
+
+#ifdef HAVE_LIBLALFRAME
+#include <lal/LALFrameIO.h>
+#include <lal/LALFrStream.h>
+#endif
+
+#include <lal/LALMalloc.h>
+
+#include <lal/TimeFreqFFT.h>
+#include <lal/Interpolate.h>
+#include <lal/Units.h>
+#include <lal/GeneratePulsarSignal.h> //Needed for time to be included
+#include <LALAppsVCSInfo.h>
+
+#include "MakeSFTsFunction.h"
+#include "ComputePSDFunction.h"
+
+#include <time.h>
+
+/** configuration-variables derived from user-variables */
+typedef struct{
+ REAL8TimeSeries *inputTS; /**< Input TimeSeries read from frame file */
+ CHAR *VCSInfoString; /**< LAL + LALapps Git version string */
+
+ REAL8TimeSeries *originalStrain; /**< Copy of strain to gate at the end */
+ REAL8TimeSeries *strainPadded; /**< Zero-padded strain */
+ REAL8FrequencySeries *welchPSD; /**< Average PSD computed with Welchs Method */
+ COMPLEX16FrequencySeries *strainWhitened; /**< Fourier transformed strain */
+ REAL8TimeSeries *magnitude; /**< Will be used and in the end to hold the magnitude */
+ REAL8FrequencySeries *referencePSD; /**< Harmonic Mean PSD of whole Obs run */
+
+ BOOLEAN have_PSDRefFile; /**< Was the PSD Reference file given as a parameter? */
+ CHAR *logFileName; /**< Filename of the logfile to be written */
+ LALStringVector *fileNames; /**< Currently only holding the logFileName */
+ LALStringVector *identifier; /**< Filename identifier for multi-processing */
+ UINT4 strainPaddedLength; /**< Todo: Remove. */
+ UINT4 padStart; /**< Before this array index strain is zero*/
+ UINT4 padEnd; /**< After this array index strain is zero*/
+ /* Gating constants: */
+ REAL8 threshHigh; /**< High threshold to identify glitch as gate */
+ REAL8 threshLow; /**< Lower threshold to identify length of glitch */
+ UINT4 nDead; /**< Minimal distance (in samples) for two glitches to be considered seperated */
+ CHAR* SFTTukeyWindowNaming; /**< Naming constant for tukey window radius in SFT file name */
+
+} ConfigVars_t;
+
+/* ----- User variables */
+typedef struct{
+ CHAR *outLabel;
+ /* Time */
+ LIGOTimeGPS startTime; /**< Start-time of requested signal in detector-frame (GPS seconds) */
+ INT4 duration; /**< Duration of requested signal in seconds */
+
+ /* Frame input and output */
+ CHAR *outDir; /**< directory for writing output timeseries in frame files */
+ LALStringVector *inFrames; /**< At the moment only ONE timeseries supported: frame glob-patterns or cache files of time-series data to be added to output (one per IFO) */
+ LALStringVector *inFrChannels;/**< At the moment only ONE timeseries supported: list of frame channels to read time-series data from (one per IFO) */
+ LALStringVector *outFrChannels;/**< At the moment only ONE timeseries supported: list of output frame channel names to write (one per IFO) */
+ LALStringVector* IFOs; /**< At the moment only ONE timeseries supported: list of detector-names "H1,H2,L1,.." or single detector*/
+ BOOLEAN outputGWF; /**< Determines if the gated timeseries will be written into a *.gwf file if value is not zero. Default is zero and therefore no file is written. */
+ CHAR *PSDReferenceFile; /**< Reference PSD for comparing 'quality' of SFT after gating */
+
+ /* SFT params */
+ REAL8 Tsft; /**< SFT time baseline Tsft */
+ REAL8 SFToverlap; /**< overlap SFTs by this many seconds */
+ REAL8 tukeyWindowRadius; /**< Tukey window radius for (only) the output SFT */
+ REAL8 SFTFreqStart; /**< Start frequency of output SFT */
+ REAL8 SFTFreqBand; /**< Frequency band of output SFT */
+ REAL8 Highpass; /** */
+
+ /* Gating constants: Named mostly like pycbc/strain.py detect_loud_glitch(...) method, minus naming conventions */
+ REAL8 psdDuration; /**< Duration of the segments for PSD estimation in seconds */
+ REAL8 psdStride; /**< Separation between PSD estimation segments in seconds */
+ REAL8 refFreqStart; /**< Frequency start for internally calculated reference PSD */
+ REAL8 refFreqBand; /**< Frequency band for internally calculated reference PSD */
+ REAL8 refPSDRatio; /**< Ratio of reference PSD versus intermediate PSD in each frequency bin. Default 1.05 */
+ REAL8 lowFreqCutoff; /**< Minimum frequency to include in the whitened strain */
+ REAL8 highFreqCutoff; /**< Maximum frequency to include in the whitened strain. ((If given, the input series is downsampled accordingly. If omitted, the Nyquist frequency is used. */
+ REAL8 upThreshold; /**< Minimum magnitude of whitened strain for considering a transient to be present */
+ REAL8 sigmaFactor; /**< Points are considered to belong to the transient if they exceed this factor times calculated sigma */
+ REAL8 secondsDeadTime; /**< Time in seconds outliers must be apart to be considered two discrete gates */
+ REAL8 padData; /**< Amount of time padded at the beginning and end of the input time series */
+ UINT4 taperLength; /**< Length of timeseries tapered which is between zeropadding and actual timeseries */
+ UINT4 corruptLength; /**< TODO: This is a derived variable. Get rid of it. Describes the amount of non-useable timeseries, e.g. taper. */
+ REAL8 autogatingTaper; /**< Taper of the Tukey Window to gate data */
+
+ UINT4 maxIterations; /**< Maximum number of iterations comparing intermediate PSD to best-guess PSD. Default 3 */
+ REAL8 iterReduceHighThresh; /**< Reduce the high threshold by this amount each iteration */
+ REAL8 iterReduceLowThreshFactor; /**< Reduce the lower threshold factor by this amount each iteration */
+ REAL8 iterIncreaseDeadtime; /**< Increase the dead time by this amount each iteration */
+
+ /* Or provide a file with all gates and gate that */
+ CHAR *externalGateFile; /**< Use this file of externally provided gates instead of finding gates. Format: t_start t_end. t_start and t_end are allowed to lie outside given time window. */
+ CHAR *additionalGateFile; /**< Additionally provide times of gates which are used together with found gates. Format: t_start t_end. t_start and t_end are allowed to lie outside given time window. */
+
+ /* Further analysis output options */
+ BOOLEAN outputPSD; /**< Output PSD of SFT */
+ BOOLEAN outputMagnitude; /**< Output magnitude TS*/
+ BOOLEAN outputHighpassedStrain; /**< Output highpassed input TS */
+ BOOLEAN outputWelchPSD; /**< Output Welch PSD estimate */
+ BOOLEAN outputUngatedSFT; /**< Output ungated SFT */
+ BOOLEAN outputpyCBCSFT; /**< Output SFT gated with pyCBC gating */
+ REAL8 pyCBCFixedDurationGate; /**< Fixed gate duration when using pyCBC gating */
+ BOOLEAN outputIntermediateSFT; /**< Output intermediate SFTs */
+
+// ---------- DEPRECATED and DEFUNCT options ----------
+} UserVariables_t;
+
+typedef struct {
+ UINT4 gateAmount; /**< Amount of gates. */
+ REAL8Vector* gateStart; /**< Vector of start times for gates (REAL8 repr. of GPS Time) */
+ REAL8Vector* gateEnd; /**<Vector of end times for gates */
+ REAL8Vector* gateDuration; /**< Vector of duration for each gate */
+} GateParams_t;
+
+// ---------- local prototypes ----------
+int XLALWriteREAL8TimeSeries2fp ( FILE *fp, const REAL8TimeSeries *TS );
+int XLALWriteREAL8FrequencySeries2fp ( FILE *fp, const REAL8FrequencySeries *FS );
+int XLALInitUserVars ( UserVariables_t *uvar, int argc, char *argv[] );
+int XLALFreeMem (UserVariables_t *uvar, ConfigVars_t *cfg, GateParams_t *gtParams );
+BOOLEAN is_directory ( const CHAR *fname );
+UINT4 nextPowerOf2(UINT4);
+
+int XLALInitgateStrain ( ConfigVars_t *cfg, UserVariables_t *uvar );
+int outputTimeSeriesToFile(UserVariables_t *uvar, ConfigVars_t *cfg, REAL8TimeSeries *Tseries);
+int taperStrain(UserVariables_t *uvar, REAL8TimeSeries *Tseries);
+int zeroPadStrain(ConfigVars_t *cfg, REAL8TimeSeries *mseries, REAL8TimeSeries **strain);
+int estimateWelchPSD(UserVariables_t *uvar, ConfigVars_t *cfg, REAL8TimeSeries *unPaddedTseries, REAL8TimeSeries *paddedTseries, REAL8FrequencySeries **psd);
+int copyReferencePSDFromWelchsPSD(UserVariables_t *uvar, ConfigVars_t *cfg, REAL8FrequencySeries **referencePSD );
+int fftAndWhitenStrain(UserVariables_t *uvar, REAL8TimeSeries *Tseries, COMPLEX16FrequencySeries **CFseries, REAL8FrequencySeries *psdFseries);
+int ifft(REAL8TimeSeries **Tseries, COMPLEX16FrequencySeries *CFseries);
+COMPLEX16FrequencySeries *WhitenCOMPLEX16FrequencySeries(UserVariables_t *uvar, REAL8 strainSampleRate, COMPLEX16FrequencySeries *fseries, const REAL8FrequencySeries *psd);
+int HighPass(ConfigVars_t *cfg, REAL8TimeSeries *Tseries, REAL8 frequency);
+int findGlitchesAndDuration(UserVariables_t *uvar, ConfigVars_t *cfg, REAL8TimeSeries *magnitude, REAL8TimeSeries *Tseries, GateParams_t *gtParams);
+int compMagnitudeStatistics(ConfigVars_t *cfg, REAL8TimeSeries *Tseries, REAL8 *mean, REAL8 *std, REAL8 *max, UINT2 segments);
+int getGlitchDuration(ConfigVars_t *cfg, REAL8TimeSeries *magnitude, GateParams_t *gtParams);
+int gateData(UserVariables_t *uvar, REAL8TimeSeries *Tseries, GateParams_t *gtParams);
+int callMakeSFT(UserVariables_t *uvar, ConfigVars_t *cfg, REAL8 freqStart, REAL8 freqBand, REAL8TimeSeries *Tseries, SFTtype **sft, BOOLEAN write2File, char* description, REAL8 tukeyWindowRadius);
+int callComputePSD(REAL8 freqStart, REAL8 freqBand, UINT2 blocksRngMed, SFTtype *sft, REAL8FrequencySeries **Fseries);
+int gateWithExternalGatesFile(UserVariables_t *uvar, ConfigVars_t *cfg, REAL8TimeSeries *Tseries);
+int loadAndMergeAdditionalGatesFile(UserVariables_t *uvar, GateParams_t *gtParams);
+int loadHMeanPSD(UserVariables_t *uvar, ConfigVars_t *cfg);
+int compareHMeanToThis(ConfigVars_t *cfg, REAL8FrequencySeries *sftPsd, REAL8 *sftRatio);
+int SpectrumInvertTruncate(REAL8FrequencySeries *spectrum, REAL8 lowfreq, UINT4 seglen, UINT4 trunclen, REAL8FFTPlan *fwdplan, REAL8FFTPlan *revplan);
+
+// ---------- debug and comparison prototypes ----------
+int producepyCBCGatedSFT(UserVariables_t *uvar, ConfigVars_t *cfg, REAL8 freqStart, REAL8 freqBand, REAL8TimeSeries *Tseries, REAL8TimeSeries *magnitude);
+int produceUngatedSFT(UserVariables_t *uvar, ConfigVars_t *cfg, REAL8 freqStart, REAL8 freqBand, REAL8TimeSeries *Tseries);
+int pycbcGateDataMethod(UserVariables_t *uvar, ConfigVars_t *cfg, REAL8TimeSeries *Magnitude, GateParams_t *gtParams);
+
+#define TRUE (1==1)
+#define FALSE (1==0)
+
+/* track memory usage under linux */
+#define TRACKMEMUSE 0
+
+/* Copied from MakeSFTs.c - https://lscsoft.docs.ligo.org/lalsuite/lalapps/_make_s_f_ts_8c_source.html */
+#if TRACKMEMUSE
+#include <sys/types.h>
+#include <unistd.h>
+void printmemuse(void);
+void printmemuse() {
+ pid_t mypid=getpid();
+ char commandline[256];
+ fflush(NULL);
+ sprintf(commandline,"cat /proc/%d/status | /bin/grep Vm | /usr/bin/fmt -140 -u", (int)mypid);
+ if ( system(commandline) ) XLALPrintError ("system() returned non-zero status\n");
+ fflush(NULL);
+ }
+#endif
+
+#define MEASURETIMING 0
+#if MEASURETIMING
+void measure_time(const char* method_name, clock_t begin, clock_t* single_begin);
+void measure_time(const char* method_name, clock_t begin, clock_t* single_begin) {
+ FILE *fp = NULL;
+ XLAL_CHECK( (fp = fopen("./timing.txt","a")) != NULL, XLAL_EFUNC, "Could not open timing file");
+ fprintf(fp, "%s %.3f %.3f\n", method_name, ((double) (clock()-begin)) / CLOCKS_PER_SEC, ((double) (clock()-*single_begin)) / CLOCKS_PER_SEC);
+ fclose(fp);
+ *single_begin = clock();
+}
+#endif
+
+/*----------------------------------------------------------------------*
+ * main function
+ *----------------------------------------------------------------------*/
+int main(int argc, char *argv[]){
+ //Compile time info
+ printf("Program build on: %s %s\n", __DATE__, __TIME__);
+ XLALClobberDebugLevel(LALMEMDBG | LALMEMTRACE | LALALLDBG);
+ printf("Debug Level: %d\n", XLALGetDebugLevel());
+
+ clock_t begin = clock();
+#if TRACKMEMUSE
+ printf("\nMemory use at startup is:\n"); printmemuse();
+#endif
+
+ UserVariables_t XLAL_INIT_DECL(uvar);
+ ConfigVars_t XLAL_INIT_DECL(cfg);
+ GateParams_t XLAL_INIT_DECL(gtParams);
+
+ XLAL_CHECK ( XLALInitUserVars (&uvar, argc, argv) == XLAL_SUCCESS, XLAL_EFUNC );
+ XLAL_CHECK ( XLALInitgateStrain ( &cfg, &uvar ) == XLAL_SUCCESS, XLAL_EFUNC );
+
+#if TRACKMEMUSE
+ printf("\nMemory use after init and reading in strain is:\n"); printmemuse();
+#endif
+ FILE* fp;
+ char buffer[100];
+
+ /* Applying Highpass */
+ XLAL_CHECK(HighPass(&cfg, cfg.inputTS, uvar.Highpass) == XLAL_SUCCESS, XLAL_EFUNC);
+
+ if (uvar.outputHighpassedStrain) {
+ snprintf(buffer, sizeof buffer, "%s%s%s%s", uvar.outDir, "/C_HighpassedStrain", *cfg.identifier->data, ".txt");
+ XLAL_CHECK( (fp = fopen(buffer,"w")) != NULL, XLAL_EFUNC, "Could not open file %s", buffer);
+ XLALWriteREAL8TimeSeries2fp(fp, cfg.inputTS);
+ fclose(fp);
+ }
+
+ XLAL_CHECK( (cfg.originalStrain = XLALCutREAL8TimeSeries(cfg.inputTS, 0, \
+ cfg.inputTS->data->length)) != NULL, XLAL_ENOMEM);
+
+#if TRACKMEMUSE
+ printf("Memory use after copy to originalStrain:\n"); printmemuse();
+#endif
+
+ /* Tapering Strain */
+ XLAL_CHECK ( taperStrain(&uvar, cfg.inputTS) == XLAL_SUCCESS, XLAL_EFUNC );
+
+ /* Downsampling strain if high freq cutoff */
+ if(XLALUserVarWasSet ( &uvar.highFreqCutoff )){
+ XLAL_CHECK ( XLALResampleREAL8TimeSeries(cfg.inputTS, 0.5 / uvar.highFreqCutoff) == XLAL_SUCCESS, XLAL_EFUNC );
+ uvar.corruptLength = (UINT4) uvar.padData * 1.0 / (cfg.inputTS->deltaT);
+#if TRACKMEMUSE
+ printf("Memory use after downsample is:\n"); printmemuse();
+#endif
+ }
+
+ /* Zero pad strain to a power-of-2 length */
+ XLAL_CHECK( zeroPadStrain(&cfg, cfg.inputTS, &cfg.strainPadded) == XLAL_SUCCESS, XLAL_EFUNC);
+
+ /* Make a copy of the TimeSeries with padding */
+ XLAL_CHECK( (cfg.magnitude = XLALCutREAL8TimeSeries(cfg.strainPadded, 0, \
+ cfg.strainPadded->data->length) ) != NULL, XLAL_ENOMEM);
+
+#if TRACKMEMUSE
+ printf("Memory use after copying into magnitude is:\n"); printmemuse();
+#endif
+
+/* Estimating PSD with Welchs Method. Use non-padded TimeSeries with cut corrupt length */
+ REAL8TimeSeries *psdTSeries;
+ UINT4 psdTStart = cfg.padStart + uvar.corruptLength;
+ UINT4 psdTDur = cfg.padEnd - cfg.padStart - 2 * uvar.corruptLength;
+
+ XLAL_CHECK( (psdTSeries = XLALCutREAL8TimeSeries(cfg.strainPadded, psdTStart, \
+ psdTDur) ) != NULL, XLAL_ENOMEM);
+
+ XLAL_CHECK ( estimateWelchPSD(&uvar, &cfg, psdTSeries, cfg.strainPadded, &cfg.welchPSD) == XLAL_SUCCESS, XLAL_EFUNC );
+ XLALDestroyREAL8TimeSeries(psdTSeries);
+
+ if (uvar.outputPSD || uvar.outputWelchPSD) {
+ snprintf(buffer, sizeof buffer, "%s%s%s%s", uvar.outDir, "/C_WelchPSD", *cfg.identifier->data, ".txt");
+ XLAL_CHECK( (fp = fopen(buffer,"w")) != NULL, XLAL_EFUNC, "Could not open file %s", buffer);
+ XLALWriteREAL8FrequencySeries2fp(fp, cfg.welchPSD);
+ fclose(fp);
+ }
+
+ /* If requested copy and save a portion of the PSD to have a reference for the intermediate resulting PSD */
+ if(!cfg.have_PSDRefFile) {
+ XLAL_CHECK( copyReferencePSDFromWelchsPSD(&uvar, &cfg, &cfg.referencePSD) == XLAL_SUCCESS, XLAL_EFUNC );
+ if (uvar.outputWelchPSD) {
+ snprintf(buffer, sizeof buffer, "%s%s%s%s", uvar.outDir, "/dbg_referencePSD", *cfg.identifier->data, ".txt");
+ XLAL_CHECK( (fp = fopen(buffer,"w")) != NULL, XLAL_EFUNC, "Could not open file %s", buffer);
+ XLALWriteREAL8FrequencySeries2fp(fp, cfg.referencePSD);
+ fclose(fp);
+ }
+ } else {
+ /* Load harmonic mean of e.g. whole obs to compare against this */
+ XLAL_CHECK(loadHMeanPSD(&uvar, &cfg) == XLAL_SUCCESS, XLAL_EFUNC);
+ }
+
+ /* FFT and whiten Strain */
+ XLAL_CHECK ( fftAndWhitenStrain(&uvar, cfg.magnitude, &cfg.strainWhitened, cfg.welchPSD) == XLAL_SUCCESS, XLAL_EFUNC );
+
+ /* iFFT Whitened Strain */
+ XLAL_CHECK ( ifft(&cfg.magnitude, cfg.strainWhitened) == XLAL_SUCCESS, XLAL_EFUNC );
+
+ /* Resize: Padding of power-2 isnt necessary anymore */
+ UINT4 magStart = cfg.padStart + uvar.corruptLength;
+ UINT4 magDur = cfg.padEnd - cfg.padStart - 2 * uvar.corruptLength;
+
+ XLAL_CHECK( (cfg.magnitude = XLALResizeREAL8TimeSeries(cfg.magnitude, magStart, \
+ magDur) ) != NULL, XLAL_ENOMEM);
+
+ if (uvar.outputMagnitude) {
+ snprintf(buffer, sizeof buffer, "%s%s%s%s", uvar.outDir, "/C_WhitenedMagnitude", *cfg.identifier->data, ".txt");
+ XLAL_CHECK( (fp = fopen(buffer,"w")) != NULL, XLAL_EFUNC, "Could not open file %s", buffer);
+ XLALWriteREAL8TimeSeries2fp(fp, cfg.magnitude);
+ fclose(fp);
+ }
+
+ /* Produce ungated or pyCBC-gated SFTs */
+ if (uvar.outputUngatedSFT) {
+ XLAL_CHECK( ( produceUngatedSFT(&uvar, &cfg, uvar.SFTFreqStart, uvar.SFTFreqBand, cfg.originalStrain)) == XLAL_SUCCESS, XLAL_EFUNC);
+ }
+ if (uvar.outputpyCBCSFT) {
+ XLAL_CHECK( ( producepyCBCGatedSFT(&uvar, &cfg, uvar.SFTFreqStart, uvar.SFTFreqBand, cfg.originalStrain, cfg.magnitude)) == XLAL_SUCCESS, XLAL_EFUNC);
+ }
+ /* Gate once with externally given gating file */
+ if ( XLALUserVarWasSet(&uvar.externalGateFile) ) {
+ XLAL_CHECK( ( gateWithExternalGatesFile(&uvar, &cfg, cfg.originalStrain)) == XLAL_SUCCESS, XLAL_EFUNC);
+ }
+
+ /* Find glitches and their duration */
+ XLAL_CHECK( findGlitchesAndDuration(&uvar, &cfg, cfg.magnitude, cfg.strainPadded, >Params) == XLAL_SUCCESS, XLAL_EFUNC);
+
+ /* Load additional gates file and merge with found gates */
+ if ( XLALUserVarWasSet(&uvar.additionalGateFile) ) {
+ XLAL_CHECK( loadAndMergeAdditionalGatesFile(&uvar, >Params) == XLAL_SUCCESS, XLAL_EFUNC);
+// gtParams = loadAndMergeAdditionalGatesFile(&uvar, >Params);
+ printf("AfterM %d gates (cross check %d %d %d)\n", gtParams.gateAmount, gtParams.gateStart->length, gtParams.gateEnd->length, gtParams.gateDuration->length);
+ fflush(stdout);
+ }
+
+ /* Gate data */
+ XLAL_CHECK( gateData(&uvar, cfg.originalStrain, >Params) == XLAL_SUCCESS, XLAL_EFUNC);
+ /* Write used gates (in the end) into file */
+ FILE *fplog = fopen(cfg.logFileName, "a");
+ if( fplog == NULL) {
+ printf("Error opening file: %s", cfg.logFileName);
+ return EXIT_FAILURE;
+ }
+ fprintf(fplog, "\n\nAfterall Gate(s) found: %d\n", gtParams.gateAmount);
+ for(UINT2 index = 0; index < gtParams.gateAmount; index++) {
+ fprintf(fplog, "# Start End Dur: %.16f %.16f %.16f\n", gtParams.gateStart->data[index], gtParams.gateEnd->data[index], gtParams.gateDuration->data[index]);
+ }
+ fclose(fplog);
+
+ char newGateFileName[100];
+ snprintf(newGateFileName, 100, "%s%s%s%s", uvar.outDir, "/Gates", *cfg.identifier->data, ".txt");
+ FILE *fpNG = fopen(newGateFileName, "w");
+ if( fpNG == NULL) {
+ printf("Error opening file: %s", newGateFileName);
+ return EXIT_FAILURE;
+ }
+// fprintf(fp, "Gate(s) found: %d\n", gtParams->gateAmount);
+ for(UINT2 index = 0; index < gtParams.gateAmount; index++) {
+ fprintf(fpNG, "%.16f %.16f %.16f\n", gtParams.gateStart->data[index], gtParams.gateEnd->data[index], gtParams.gateDuration->data[index]);
+ }
+ fclose(fpNG);
+
+ if (uvar.outputHighpassedStrain) {
+ snprintf(buffer, sizeof buffer, "%s%s%s%s", uvar.outDir, "/C_B4SecHighpassStrain", *cfg.identifier->data, ".txt");
+ XLAL_CHECK( (fp = fopen(buffer,"w")) != NULL, XLAL_EFUNC, "Could not open file %s", buffer);
+ XLALWriteREAL8TimeSeries2fp(fp, cfg.originalStrain);
+ fclose(fp);
+ }
+
+ /* Highpass again */
+ XLAL_CHECK(HighPass(&cfg, cfg.originalStrain, uvar.Highpass) == XLAL_SUCCESS, XLAL_EFUNC);
+
+/* if (uvar.outputHighpassedStrain) {
+ snprintf(buffer, sizeof buffer, "%s%s%s%s", uvar.outDir, "/C_AfterSecHighpassStrain", *cfg.identifier->data, ".txt");
+ XLAL_CHECK( (fp = fopen(buffer,"w")) != NULL, XLAL_EFUNC, "Could not open file %s", buffer);
+ XLALWriteREAL8TimeSeries2fp(fp, cfg.originalStrain);
+ fclose(fp);
+ }
+*/
+ /* Remove Padding */
+ XLAL_CHECK( (cfg.originalStrain = XLALResizeREAL8TimeSeries(cfg.originalStrain,
+ uvar.padData / cfg.originalStrain->deltaT, (uvar.duration - 2*uvar.padData) / cfg.originalStrain->deltaT) ) != NULL, XLAL_ENOMEM);
+
+//#if (defined DEBUGWRITEOUT || defined DEBUGGATEDSTRAIN)
+// snprintf(buffer, sizeof buffer, "%s%s%s%s", uvar.outDir, "/C_GatedStrain", *cfg.identifier->data, ".txt");
+// XLAL_CHECK( (fp = fopen(buffer,"w")) != NULL, XLAL_EFUNC, "Could not open file %s", buffer);
+// XLALWriteREAL8TimeSeries2fp(fp, cfg.originalStrain);
+// fclose(fp);
+//#endif
+
+
+ //Output timeseries into *.gwf
+ if (uvar.outputGWF) {
+ XLAL_CHECK ( outputTimeSeriesToFile(&uvar, &cfg, cfg.originalStrain) == XLAL_SUCCESS, XLAL_EFUNC );
+ }
+
+ /* Make the gated SFT */
+ SFTtype *sft = NULL;
+ char description[32];
+ snprintf(description, 32, "GATED_%s", cfg.SFTTukeyWindowNaming);
+ XLAL_CHECK( callMakeSFT(&uvar, &cfg, uvar.SFTFreqStart, uvar.SFTFreqBand, cfg.originalStrain, &sft, 1, description, uvar.tukeyWindowRadius) == XLAL_SUCCESS, XLAL_EFUNC);
+ XLALDestroyCOMPLEX8FrequencySeries(sft);
+
+ FILE *fpend = fopen(cfg.logFileName, "a");
+ if( fpend == NULL) {
+ printf("Error opening file: %s", cfg.logFileName);
+ return EXIT_FAILURE;
+ }
+ fprintf(fpend, "\nSuccessfully leaving gateStrain program... ");
+ fclose(fpend);
+
+ XLALFreeMem (&uvar, &cfg, >Params );
+ LALCheckMemoryLeaks();
+ printf("\n------------------------------------------------\n");
+#if TRACKMEMUSE
+ printf("Memory use before program termination is:\n"); printmemuse();
+#endif
+
+ double time_spent = (double) (clock()-begin) / CLOCKS_PER_SEC;
+ printf("\n#################################################\n");
+ printf("Successfully leaving gateStrain program in %.3f ... ", time_spent);
+ printf("\n#################################################\n");
+ return 0;
+} /* Main */
+
+/**
+ * Make SFT of non-gated time series.
+ */
+int produceUngatedSFT(UserVariables_t *uvar, ConfigVars_t *cfg, REAL8 freqStart, REAL8 freqBand, REAL8TimeSeries *Tseries) {
+ REAL8TimeSeries *strainCopy = NULL;
+ SFTtype *sft = NULL;
+ char description[32];
+ snprintf(description, 32, "UNGATED_%s", cfg->SFTTukeyWindowNaming);
+
+ XLAL_CHECK( ( strainCopy = XLALCutREAL8TimeSeries(Tseries, 0, uvar->duration / Tseries->deltaT) ) != NULL, XLAL_ENOMEM);
+
+ /* Highpass again */
+ XLAL_CHECK(HighPass(cfg, strainCopy, uvar->Highpass) == XLAL_SUCCESS, XLAL_EFUNC);
+
+ /* Remove Padding */
+ XLAL_CHECK( (strainCopy = XLALResizeREAL8TimeSeries(strainCopy,
+ uvar->padData / strainCopy->deltaT, (uvar->duration - 2*uvar->padData) / strainCopy->deltaT) ) != NULL, XLAL_ENOMEM);
+
+ XLAL_CHECK( callMakeSFT(uvar, cfg, freqStart, freqBand, strainCopy, &sft, uvar->outputUngatedSFT, description, uvar->tukeyWindowRadius) == XLAL_SUCCESS, XLAL_EFUNC);
+
+ XLALDestroyREAL8TimeSeries(strainCopy);
+ XLALDestroyCOMPLEX8FrequencySeries(sft);
+ return XLAL_SUCCESS;
+}
+
+/**
+ * Make SFT with pyCBC gating methode.
+ */
+int producepyCBCGatedSFT(UserVariables_t *uvar, ConfigVars_t *cfg, REAL8 freqStart, REAL8 freqBand, REAL8TimeSeries *Tseries, REAL8TimeSeries *magnitude) {
+ GateParams_t XLAL_INIT_DECL(pyCBCgtParams);
+ REAL8TimeSeries *strainCopy = NULL;
+ REAL8TimeSeries *workCopyMag = NULL;
+ SFTtype *sft = NULL;
+ char description[32];
+
+ XLAL_CHECK( ( strainCopy = XLALCutREAL8TimeSeries(Tseries, 0, uvar->duration / Tseries->deltaT) ) != NULL, XLAL_ENOMEM);
+
+ XLAL_CHECK( (workCopyMag = XLALCutREAL8TimeSeries(magnitude, 0, \
+ magnitude->data->length) ) != NULL, XLAL_ENOMEM);
+
+ for(UINT4 indexx = 0; indexx < workCopyMag->data->length; indexx++) {
+ workCopyMag->data->data[indexx] = fabs(workCopyMag->data->data[indexx]);
+ }
+
+ XLAL_CHECK ( pycbcGateDataMethod(uvar, cfg, workCopyMag, &pyCBCgtParams)
+ == XLAL_SUCCESS, XLAL_EFUNC );
+ if(pyCBCgtParams.gateAmount != 0 ) {
+ XLAL_CHECK( gateData(uvar, strainCopy, &pyCBCgtParams) == XLAL_SUCCESS, XLAL_EFUNC);
+
+/* FILE* fpold;
+ char bufferold[100];
+ snprintf(bufferold, sizeof bufferold, "%s%s%s%s", uvar.outDir, "/C_strainCopy", *cfg->identifier->data, ".txt");
+ fpold = fopen(bufferold,"w");
+ XLALWriteREAL8TimeSeries2fp(fpold, strainCopy);
+ fclose(fpold);
+ snprintf(bufferold, sizeof bufferold, "%s%s%s%s", uvar.outDir, "/C_MagnitudeStrain", *cfg->identifier->data, ".txt");
+ fpold = fopen(bufferold,"w");
+ XLALWriteREAL8TimeSeries2fp(fpold, workCopyMag);
+ fclose(fpold);*/
+
+ } else {
+ printf("No (old) gate found");
+ }
+
+ /* Highpass again */
+ XLAL_CHECK(HighPass(cfg, strainCopy, uvar->Highpass) == XLAL_SUCCESS, XLAL_EFUNC);
+
+ /* Remove Padding */
+ XLAL_CHECK( (strainCopy = XLALResizeREAL8TimeSeries(strainCopy,
+ uvar->padData / strainCopy->deltaT, (uvar->duration - 2*uvar->padData) / strainCopy->deltaT) ) != NULL, XLAL_ENOMEM);
+
+ snprintf(description, 32, "pyCBCGATED_%s", cfg->SFTTukeyWindowNaming);
+ XLAL_CHECK( callMakeSFT(uvar, cfg, freqStart, freqBand, strainCopy, &sft, uvar->outputpyCBCSFT, description, uvar->tukeyWindowRadius) == XLAL_SUCCESS, XLAL_EFUNC);
+
+ XLALDestroyCOMPLEX8FrequencySeries(sft);
+ XLALDestroyREAL8TimeSeries(workCopyMag);
+ XLALDestroyREAL8TimeSeries(strainCopy);
+ XLALDestroyREAL8Vector(pyCBCgtParams.gateStart);
+ XLALDestroyREAL8Vector(pyCBCgtParams.gateEnd);
+ XLALDestroyREAL8Vector(pyCBCgtParams.gateDuration);
+
+ return XLAL_SUCCESS;
+}
+
+/**
+ * Register all "user-variables", and initialized them from command-line and config-files
+ */
+int XLALInitUserVars ( UserVariables_t *uvar, int argc, char *argv[] ) {
+ XLAL_CHECK ( uvar != NULL, XLAL_EINVAL, "Invalid NULL input 'uvar'\n");
+ XLAL_CHECK ( argv != NULL, XLAL_EINVAL, "Invalid NULL input 'argv'\n");
+
+#define MISC_DEFAULT "gating"
+ XLAL_CHECK ( (uvar->outLabel = XLALStringDuplicate ( MISC_DEFAULT )) != NULL, XLAL_EFUNC );
+
+ /* ---------- set a few defaults ---------- */
+ uvar->Tsft = 1800;
+ uvar->Highpass = 10.0;
+ uvar->psdDuration = 4.0;
+ uvar->psdStride = 2.0;
+ uvar->refFreqStart = 33.0;
+ uvar->refFreqBand = 5.0;
+ uvar->refPSDRatio = 1.05;
+ uvar->lowFreqCutoff = 15;
+ uvar->upThreshold = 50;
+ uvar->sigmaFactor = 6;
+ uvar->secondsDeadTime = 3.0;
+ uvar->padData = 8.0;
+ uvar->autogatingTaper = 0.25;
+ uvar->tukeyWindowRadius = 0.001;
+ uvar->SFTFreqStart = 15.0;
+ uvar->SFTFreqBand = 1990.0;
+
+ uvar->maxIterations = 3;
+ uvar->iterReduceHighThresh = 10.0;
+ uvar->iterReduceLowThreshFactor = 0.01;
+ uvar->iterIncreaseDeadtime = 0.1;
+ // uvar->fmin = 0; /* no heterodyning by default */
+ //uvar->Band = 8192; /* 1/2 LIGO sampling rate by default */
+ uvar->pyCBCFixedDurationGate = 16.0;
+
+ /* start + duration of timeseries */
+ XLALRegisterUvarMember( startTime, EPOCH, 'G', REQUIRED, "Start-time of requested signal in detector-frame (format ''xx.yy[GPS|MJD]')");
+ XLALRegisterUvarMember( duration, INT4, 0, REQUIRED, "Duration of requested signal in seconds");
+
+ /* frame input/output options */
+#ifdef HAVE_LIBLALFRAME
+ XLALRegisterUvarMember ( outDir, STRING, 'F', REQUIRED, "Output Dir: directory for output timeseries frame files and short Fourier transforms");
+ XLALRegisterUvarMember ( inFrames, STRINGVector,'C', REQUIRED, "CSV list (one per IFO) of input frame cache files");
+ XLALRegisterUvarMember ( inFrChannels, STRINGVector,'N', REQUIRED, "CSV list (one per IFO) of frame channels to read timeseries from");
+ XLALRegisterUvarMember ( outFrChannels, STRINGVector, 0, REQUIRED, "CSV list (one per IFO) of output frame channel names [default: <inFrChannels>-<outLabel> or <IFO>:<outLabel>]");
+#else
+ XLALRegisterUvarMember ( outDir, STRING, 'F', DEFUNCT, "Need to compile with lalframe support for this option to work");
+ XLALRegisterUvarMember ( inFrames, STRINGVector, 'C', DEFUNCT, "Need to compile with lalframe support for this option to work");
+ XLALRegisterUvarMember ( inFrChannels, STRINGVector, 'N', DEFUNCT, "Need to compile with lalframe support for this option to work");
+ XLALRegisterUvarMember ( outFrChannels, STRINGVector, 0, DEFUNCT, "Need to compile with lalframe support for this option to work");
+#endif
+
+ XLALRegisterUvarMember ( tukeyWindowRadius, REAL8, 0, OPTIONAL, "Tukey window radius for (only) the output SFT. Default 0.001");
+ XLALRegisterUvarMember ( SFTFreqStart, REAL8, 0, OPTIONAL, "Start frequency of output SFT");
+ XLALRegisterUvarMember ( SFTFreqBand, REAL8, 0, OPTIONAL, "Frequency band of output SFT");
+ XLALRegisterUvarMember ( Highpass, REAL8, 0, OPTIONAL, "High pass filtering frequency in Hz.");
+
+ /* Gating options */
+ XLALRegisterUvarMember ( psdDuration, REAL8, 0, OPTIONAL, "Duration of the segments for PSD estimation in seconds");
+ XLALRegisterUvarMember ( psdStride, REAL8, 0, OPTIONAL, "Separation between PSD estimation segments in seconds");
+ XLALRegisterUvarMember ( refFreqStart, REAL8, 0, OPTIONAL, "Frequency start for internally calculated reference PSD. Default 33.0 (Hz)");
+ XLALRegisterUvarMember ( refFreqBand, REAL8, 0, OPTIONAL, "Frequency band for internally calculated reference PSD. Default 5.0 (Hz)");
+ XLALRegisterUvarMember ( refPSDRatio, REAL8, 0, OPTIONAL, "Ratio of reference PSD versus intermediate PSD in each frequency bin. Default 1.05");
+ XLALRegisterUvarMember ( lowFreqCutoff, REAL8, 0, OPTIONAL, "Minimum frequency to include in the whitened strain");
+ XLALRegisterUvarMember ( highFreqCutoff, REAL8, 0, OPTIONAL, "Maximum frequency to include in the whitened strain");
+ XLALRegisterUvarMember ( upThreshold, REAL8, 0, OPTIONAL, "Minimum magnitude of whitened strain for considering a transient to be present");
+ XLALRegisterUvarMember ( sigmaFactor, REAL8, 0, OPTIONAL, "Points are considered to belong to the transient if they exceed this factor times calculated sigma");
+ XLALRegisterUvarMember ( secondsDeadTime, REAL8, 0, OPTIONAL, "Time in seconds outliers must be apart to be considered two discrete gates");
+ XLALRegisterUvarMember ( padData, REAL8, 0, OPTIONAL, "Amount of time padded at the beginning and end of the input time series. Default 8.0");
+ XLALRegisterUvarMember ( autogatingTaper, REAL8, 0, OPTIONAL, "Taper for the tukey window to gate data with. Default 0.25");
+ XLALRegisterUvarMember ( outputGWF, BOOLEAN, 0, OPTIONAL, "Write gated time series to *.gwf file if value is not zero. Default is zero");
+ XLALRegisterUvarMember ( PSDReferenceFile, STRING, 0, OPTIONAL, "Input reference PSD which is used to compare the SFT 'quality' after gating. File consists of freq psdvalue pairs. Overwrites refFreqStart & refFreqBand option.");
+
+ XLALRegisterUvarMember ( maxIterations, UINT4, 0, OPTIONAL, "Maximum number of iterations comparing intermediate PSD to best-guess PSD. Default 3");
+ XLALRegisterUvarMember ( iterReduceHighThresh, REAL8, 0, OPTIONAL, "Reduce the high threshold by this amount each iteration");
+ XLALRegisterUvarMember ( iterReduceLowThreshFactor, REAL8, 0, OPTIONAL, "Reduce the lower threshold factor by this amount each iteration");
+ XLALRegisterUvarMember ( iterIncreaseDeadtime, REAL8, 0, OPTIONAL, "Increase the dead time by this amount each iteration");
+
+ XLALRegisterUvarMember ( outputPSD, BOOLEAN, 0, DEVELOPER, "Output PSD of SFT");
+ XLALRegisterUvarMember ( outputMagnitude, BOOLEAN, 0, DEVELOPER, "Output magnitude time series");
+ XLALRegisterUvarMember ( outputHighpassedStrain, BOOLEAN, 0, DEVELOPER, "Output highpassed input time series");
+ XLALRegisterUvarMember ( outputWelchPSD, BOOLEAN, 0, DEVELOPER, "Output Welch PSD estimate");
+ XLALRegisterUvarMember ( outputUngatedSFT, BOOLEAN, 0, DEVELOPER, "Output ungated SFT");
+ XLALRegisterUvarMember ( outputpyCBCSFT, BOOLEAN, 0, DEVELOPER, "Output SFT gated with pyCBC gating");
+ XLALRegisterUvarMember ( pyCBCFixedDurationGate, REAL8, 0, DEVELOPER, "Fixed gate duration when using pyCBC gating");
+ XLALRegisterUvarMember ( outputIntermediateSFT, BOOLEAN, 0, DEVELOPER, "Output intermediate SFTs");
+
+ XLALRegisterUvarMember ( externalGateFile, STRING, 0, OPTIONAL, "Use this file of externally provided gates instead of finding gates. Format: t_start t_end. t_start and t_end are allowed to lie outside given time window.");
+ XLALRegisterUvarMember ( additionalGateFile, STRING, 0, OPTIONAL, "Additionally provide times of gates which are used together with found gates. Format: t_start t_end. t_start and t_end are allowed to lie outside given time window.");
+
+ /* read cmdline & cfgfile */
+ BOOLEAN should_exit = 0;
+ XLAL_CHECK( XLALUserVarReadAllInput( &should_exit, argc, argv, lalAppsVCSInfoList ) == XLAL_SUCCESS, XLAL_EFUNC );
+ if ( should_exit ) {
+ exit (1);
+ }
+
+ /* Test if output directory is valid */
+ XLAL_CHECK( is_directory(uvar->outDir), XLAL_EFUNC, "Error: %s is not a valid existing directory", uvar->outDir);
+
+ return XLAL_SUCCESS;
+} /* XLALInitUserVars() */
+
+
+/**
+ * Copied from makefakedata_v5, edited: https://lscsoft.docs.ligo.org/lalsuite/lalapps/makefakedata__v5_8c_source.html
+ * Handle user-input and set up shop accordingly, and do all
+ * consistency-checks on user-input.
+ */
+int
+XLALInitgateStrain ( ConfigVars_t *cfg, UserVariables_t *uvar ){
+ XLAL_CHECK ( cfg != NULL, XLAL_EINVAL, "Invalid NULL input 'cfg'\n" );
+ XLAL_CHECK ( uvar != NULL, XLAL_EINVAL, "Invalid NULL input 'uvar'\n");
+
+ cfg->VCSInfoString = XLALVCSInfoString(lalAppsVCSInfoList, 0, "%% ");
+ XLAL_CHECK ( cfg->VCSInfoString != NULL, XLAL_EFUNC, "XLALVCSInfoString failed.\n" );
+ // ---------- check user-input consistency ----------
+
+ // ----- check if frames + frame channels given
+ BOOLEAN have_frames = (uvar->inFrames != NULL);
+ BOOLEAN have_channels= (uvar->inFrChannels != NULL);
+ XLAL_CHECK ( have_frames && have_channels, XLAL_EINVAL, "Need both --inFrames and --inFrChannels\n");
+
+ // ----- Time: From --startTime+duration
+ BOOLEAN have_startTime = XLALUserVarWasSet ( &uvar->startTime );
+ BOOLEAN have_duration = XLALUserVarWasSet ( &uvar->duration );
+ uvar->Tsft = 1800;
+ uvar->SFToverlap = 0;
+ // need BOTH startTime+duration
+ XLAL_CHECK ( have_duration && have_startTime, XLAL_EINVAL, "Need BOTH {--startTime,--duration} \n");
+
+ /* High freq cutoff has to be a power of 2 */
+ if ( XLALUserVarWasSet ( &uvar->highFreqCutoff ) ) {
+ XLAL_CHECK ( ((UINT2) uvar->highFreqCutoff & ( (UINT2) uvar->highFreqCutoff - 1)) == 0, XLAL_EINVAL, "--highFreqCutoff has to be a power of 2\n");
+ }
+ /* Handle reference PSD options */
+ if ( XLALUserVarWasSet ( &uvar->PSDReferenceFile ) ) {
+ cfg->have_PSDRefFile = (uvar->PSDReferenceFile != NULL);
+ } else {
+ cfg->have_PSDRefFile = (uvar->PSDReferenceFile != NULL);
+ }
+ XLAL_CHECK(( XLALUserVarWasSet ( &uvar->refFreqStart ) == XLALUserVarWasSet( &uvar->refFreqBand )), XLAL_EINVAL, "Need BOTH or NONE (--refFreqStart, --refFreqBand)\n");
+
+ /* Handle iteration parameters */
+ XLAL_CHECK(( cfg->threshHigh - uvar->maxIterations * uvar->iterReduceHighThresh < 0.0 ), XLAL_EINVAL, "Combination of high threshold, maximum iterations and threshold reduction leads to negative values. \n");
+ XLAL_CHECK(( cfg->threshHigh - uvar->maxIterations * uvar->iterReduceHighThresh < uvar->sigmaFactor ), XLAL_EINVAL, "Combination of high threshold, maximum iterations, threshold reduction and lower threshold factor likely leads to high threshold values smaller than the lower threshold factor at max iterations. \n");
+
+ char descWin[16];
+ snprintf(descWin, 16, "WIN%g", uvar->tukeyWindowRadius);
+ memmove(descWin+3, descWin+3+2, strlen(descWin)-3);
+// cfg->SFTTukeyWindowNaming = malloc(strlen(descWin) + 1);
+ cfg->SFTTukeyWindowNaming = strdup(descWin);
+ printf("SFT tukey window naming: %s\n", cfg->SFTTukeyWindowNaming);
+
+#ifdef HAVE_LIBLALFRAME
+ UINT4 numDetectors = uvar->inFrChannels->length;
+ XLAL_CHECK ( uvar->inFrames->length == numDetectors, XLAL_EINVAL, "Need equal number of channel names (%d) as frame specifications (%d)\n", uvar->inFrChannels->length, numDetectors );
+
+ //XLAL_CHECK ( (cfg->inputMultiTS = XLALCalloc ( 1, sizeof(*cfg->inputMultiTS))) != NULL, XLAL_ENOMEM );
+ /// cfg->inputMultiTS->length = numDetectors;
+ // XLAL_CHECK ( (cfg->inputMultiTS->data = XLALCalloc ( numDetectors, sizeof(cfg->inputMultiTS->data[0]) )) != NULL, XLAL_ENOMEM );
+ //if ( cfg->multiTimestamps == NULL ){
+ // XLAL_CHECK ( (cfg->multiTimestamps = XLALCalloc ( 1, sizeof(*cfg->multiTimestamps) )) != NULL, XLAL_ENOMEM );
+ // XLAL_CHECK ( (cfg->multiTimestamps->data = XLALCalloc ( numDetectors, sizeof(cfg->multiTimestamps->data[0]))) != NULL, XLAL_ENOMEM );
+ // cfg->multiTimestamps->length = numDetectors;
+ //}
+ for ( UINT4 X = 0; X < numDetectors; X ++ ) {
+ LALCache *cache;
+ XLAL_CHECK ( (cache = XLALCacheImport ( uvar->inFrames->data[X] )) != NULL, XLAL_EFUNC, "Failed to import cache file '%s'\n", uvar->inFrames->data[X] );
+ // ----- open frame stream ----------
+ LALFrStream *stream;
+ XLAL_CHECK ( (stream = XLALFrStreamCacheOpen ( cache )) != NULL, XLAL_EFUNC, "Failed to open stream from cache file '%s'\n", uvar->inFrames->data[X] );
+ XLALDestroyCache ( cache );
+
+ // Create identifier for file saving identifcation
+ LALStringVector* ParsedString;
+ XLAL_CHECK( ( ParsedString = XLALParseStringVector(uvar->inFrames->data[X], "/")) != NULL, XLAL_EFUNC);
+ XLAL_CHECK( ( cfg->identifier = XLALCreateStringVector(
+ ParsedString->data[ParsedString->length-1], NULL)) != NULL, XLAL_EFUNC);
+ printf("Identifier is: %s\n", *cfg->identifier->data);
+ char buffer[100];
+ snprintf(buffer, 100, "%s%s%s%s", uvar->outDir, "/log_gating", *cfg->identifier->data, ".log");
+ XLAL_CHECK( ( cfg->fileNames = XLALCreateStringVector(buffer, NULL)) != NULL, XLAL_EFUNC);
+ cfg->logFileName = cfg->fileNames->data[0];
+ XLALDestroyStringVector(ParsedString);
+
+ // ----- determine time span of frame data from the stream ----------
+ // determine time spanned by the frames in this cache
+ LIGOTimeGPS frames_startGPS, frames_endGPS;
+ XLAL_CHECK ( XLALFrStreamSeekO ( stream, 0, SEEK_SET ) == 0, XLAL_EFUNC, "Failed to move to start of input frame-stream\n");
+ XLAL_CHECK ( XLALFrStreamTell ( &frames_startGPS, stream ) == XLAL_SUCCESS, XLAL_EFUNC );
+ XLAL_CHECK ( XLALFrStreamSeekO ( stream, 0, SEEK_END ) == 0, XLAL_EFUNC, "Failed to move to end of input frame-stream\n");
+ XLAL_CHECK ( XLALFrStreamTell ( &frames_endGPS, stream ) == XLAL_SUCCESS, XLAL_EFUNC );
+ REAL8 frames_start = XLALGPSGetREAL8 ( &frames_startGPS );
+ REAL8 frames_end = XLALGPSGetREAL8 ( &frames_endGPS );
+ // REAL8 frames_span = frames_end - frames_start;
+
+ // Check if given times are correct
+ XLAL_CHECK ( frames_start <= XLALGPSGetREAL8( &uvar->startTime ), XLAL_ETIME, "Frame start is after given time starts\n");
+ XLAL_CHECK ( frames_end >= XLALGPSGetREAL8( &uvar->startTime ) + uvar->duration, XLAL_ETIME, "Start + Duration exceeds frame end\n" );
+ /*
+ LIGOTimeGPS ts_startGPS = uvar->startTime;
+ REAL8 ts_duration = uvar->duration;
+ XLAL_CHECK ( (cfg->multiTimestamps->data[X] = XLALMakeTimestamps ( ts_startGPS, ts_duration, uvar->Tsft, uvar->SFToverlap ) ) != NULL, XLAL_EFUNC );
+ // in this case we can't allow timestamps extending beyond the end of the frame timeseries, so we drop the last timestamp if necessary
+ LIGOTimeGPSVector *timestampsX = cfg->multiTimestamps->data[X];
+ REAL8 tEnd;
+ while ( (tEnd = XLALGPSGetREAL8 ( ×tampsX->data[timestampsX->length-1] ) + timestampsX->deltaT) > frames_end ) {
+ timestampsX->length --;
+ }
+ */
+ const char *channel = uvar->inFrChannels->data[X];
+ size_t limit = 0; // unlimited read
+// REAL8TimeSeries *ts;
+ //.... Use given time instead of frame time
+ XLAL_CHECK ( (cfg->inputTS = XLALFrStreamInputREAL8TimeSeries ( stream, channel, &uvar->startTime, (REAL8) uvar->duration, limit )) != NULL,
+ XLAL_EFUNC, "Frame reading failed for stream created for '%s': ts_start = {%d,%d}, duration=%.0f\n",
+ uvar->inFrames->data[X], uvar->startTime.gpsSeconds, uvar->startTime.gpsNanoSeconds, (REAL8) uvar->duration );
+
+ XLAL_CHECK ( XLALFrStreamClose ( stream ) == XLAL_SUCCESS, XLAL_EFUNC, "Stream closing failed for cache file '%s'\n", uvar->inFrames->data[X] );
+ } // for X < numDetectors
+#else
+ printf("NEED LIBLALFRAME TO FUNCTION. RECOMPILE");
+#endif
+
+ /* Calculate user-input-dependant variables */
+ uvar->taperLength = (UINT4) (uvar->padData / (cfg->inputTS->deltaT));
+
+ uvar->corruptLength = (UINT4) uvar->padData * 1.0 / (cfg->inputTS->deltaT);
+ return XLAL_SUCCESS;
+} /* XLALInitgateStrain() */
+
+/**
+ * Copied partly from MakeSFTs.c - https://lscsoft.docs.ligo.org/lalsuite/lalapps/_make_s_f_ts_8c_source.html
+ * Highpass a timeseries with a butterworth high pass
+ */
+int HighPass(ConfigVars_t *cfg, REAL8TimeSeries *Tseries, REAL8 frequency) {
+ if (frequency <= 0) {
+ return XLAL_SUCCESS;
+ }
+ PassBandParamStruc filterpar;
+ char tmpname[] = "Butterworth High Pass";
+
+ filterpar.name = tmpname;
+ filterpar.nMax = 8;
+ filterpar.f2 = frequency;
+ filterpar.a2 = 0.9;
+ filterpar.f1 = -1.0;
+ filterpar.a1 = -1.0;
+
+ XLAL_CHECK( ( XLALButterworthREAL8TimeSeries(Tseries, &filterpar)) == XLAL_SUCCESS, XLAL_EFUNC);
+
+ #if TRACKMEMUSE
+ printf("Memory use after Highpass:\n"); printmemuse();
+ #endif
+ cfg = cfg;
+ return XLAL_SUCCESS;
+}
+
+/**
+ * Taper strain at the start and the end
+ */
+int taperStrain(UserVariables_t *uvar, REAL8TimeSeries *Tseries) {
+ UINT4 endIndex = Tseries->data->length;
+ for(UINT4 index = 0; index < uvar->taperLength; index++) {
+ Tseries->data->data[index] *= ((REAL8) index) / uvar->taperLength;
+ Tseries->data->data[endIndex - index - 1] *= ((REAL8) index) / uvar->taperLength;
+ }
+
+#if TRACKMEMUSE
+ printf("Memory use after taperStrain:\n"); printmemuse();
+#endif
+ return XLAL_SUCCESS;
+} /* taperStrain() */
+
+/**
+ * Zero-pad strain to a power-of-2 length. As of now, this transforms the multiTimeSeries into the cfg->strain (single) timeseries
+ * TODO: Think about using XLALResize (but this pads zeroes only to the start)
+
+ */
+int zeroPadStrain(ConfigVars_t *cfg, REAL8TimeSeries *inTseries, REAL8TimeSeries **strain) {
+ REAL8TimeSeries *Tseries;
+ cfg->strainPaddedLength = nextPowerOf2(inTseries->data->length);
+ UINT4 padStart = cfg->strainPaddedLength / 2 - inTseries->data->length / 2;
+ UINT4 padEnd = padStart + inTseries->data->length;
+
+ //Changing the epoch here
+ LIGOTimeGPS XLAL_INIT_DECL(newTime);
+ newTime = *XLALGPSSetREAL8(&newTime, XLALGPSGetREAL8(&inTseries->epoch));
+
+ XLALGPSAdd(&newTime, ( (REAL8) padStart) * inTseries->deltaT * -1.0);
+ XLAL_CHECK ( (Tseries = XLALCreateREAL8TimeSeries(
+ "Zero padded strain",
+ &newTime,
+ 0,
+ inTseries->deltaT,
+ &lalStrainUnit,
+ cfg->strainPaddedLength)) != NULL, XLAL_EFUNC);
+
+ /* Calloc has set the array to zero. Skip index 0 to pad start. Same for pad end to end of array.*/
+ /** Set to zero anyway, better to be safe than sorry. **/
+ for(UINT4 index = 0; index < padStart; index++) {
+ Tseries->data->data[index] = 0;
+ Tseries->data->data[cfg->strainPaddedLength-1-index] = 0;
+ }
+ for(UINT4 index = padStart; index < padEnd ; index++){
+ Tseries->data->data[index] = inTseries->data->data[index-padStart];
+ }
+
+ *strain = Tseries;
+ cfg->padStart = padStart;
+ cfg->padEnd = padEnd;
+
+ XLALDestroyREAL8TimeSeries(inTseries);
+
+#if TRACKMEMUSE
+ printf("Memory use after zeroPadStrain is:\n"); printmemuse();
+#endif
+ return XLAL_SUCCESS;
+} /* zeroPadStrain() */
+
+/**
+ * Estimating the PSD of the timeseries with Welchs Method
+ * See: https://lscsoft.docs.ligo.org/lalsuite/lal/group___time_freq_f_f_t__h.html
+ */
+int estimateWelchPSD(UserVariables_t *uvar, ConfigVars_t *cfg, REAL8TimeSeries *unpaddedTseries, REAL8TimeSeries *Tseries, REAL8FrequencySeries **psdFseries){
+ UINT4 halflength = (UINT4) (uvar->psdDuration / Tseries->deltaT / 2 + 1);
+ REAL8Window* hann;
+ REAL8FFTPlan* plan;
+ XLAL_CHECK( (hann = XLALCreateHannREAL8Window((UINT4) uvar->psdDuration / Tseries->deltaT )) != NULL, XLAL_ENOMEM);
+ XLAL_CHECK( (plan = XLALCreateForwardREAL8FFTPlan((UINT4) uvar->psdDuration / Tseries->deltaT, 0)) != NULL, XLAL_ENOMEM);
+
+ REAL8FrequencySeries *psd;
+ XLAL_CHECK( (psd = XLALCreateREAL8FrequencySeries("WELCHPSD",
+ &unpaddedTseries->epoch,
+ 0,
+ 0,
+ &lalStrainUnit,
+ halflength)) != NULL, XLAL_EFUNC);
+
+ /* For making the PSD, use the unpadded - corruptLength corrected - TimeSeries */
+ XLAL_CHECK( XLALREAL8AverageSpectrumMedian(psd, unpaddedTseries,
+ (UINT4) uvar->psdDuration / unpaddedTseries->deltaT,
+ (UINT4) uvar->psdStride / unpaddedTseries->deltaT,hann, plan) == XLAL_SUCCESS, XLAL_EFUNC);
+ /* Afterwards use the zero padded Timeseries (for interpolating length etc) */
+ XLALDestroyREAL8Window(hann);
+ XLALDestroyREAL8FFTPlan(plan);
+
+ /* Interpolate PSD */
+ REAL8Sequence *x_in, *y_in, *x_out, *y_out;
+ XLAL_CHECK( (x_in = XLALCreateREAL8Sequence(psd->data->length) ) != NULL, XLAL_ENOMEM);
+ /* Populate x_in-array with frequencies */
+ REAL8 xinf0 = psd->f0;
+ REAL8 xinDeltaF = psd->deltaF;
+ for(UINT4 index = 0; index < psd->data->length; index++) {
+ x_in->data[index] = xinf0 + index * xinDeltaF;
+ }
+ /* y_in are the psd_values */
+ UINT4 new_length = (psd->data->length-1) * psd->deltaF * Tseries->data->length * Tseries->deltaT + 1;
+ y_in = psd->data;
+ XLAL_CHECK( (x_out = XLALCreateREAL8Sequence(new_length) ) != NULL, XLAL_ENOMEM);
+ /* Populate new x_array with new frequencies ...*/
+ REAL8 xoutf0 = psd->f0;
+ REAL8 xoutDeltaF = 1.0 / (Tseries->data->length * Tseries->deltaT);
+ // printf("Xoutf0 %f and xoutDeltaF %f\n", xoutf0, xoutDeltaF);
+
+ for(UINT4 index = 0; index < new_length; index++) {
+ x_out->data[index] = xoutf0 + index * xoutDeltaF;
+ }
+
+ XLAL_CHECK( (y_out = XLALCreateREAL8Sequence(new_length) ) != NULL, XLAL_ENOMEM);
+ /* Do Interpolation.*/
+ const gsl_interp_type* gsl_interp_p = gsl_interp_linear;
+ XLAL_CHECK( XLALREAL8Interpolation(x_in, y_in, x_out, y_out, psd->data->length , gsl_interp_p) == XLAL_SUCCESS, XLAL_EFUNC);
+
+ XLALDestroyREAL8Sequence(y_in);
+ XLALDestroyREAL8Sequence(x_in);
+ XLALDestroyREAL8Sequence(x_out);
+ psd->data = y_out;
+
+ /* New delta f is 0.5xSampleRate/(LengthOfStrain/2+1) */
+ psd->data->length = new_length;
+ psd->deltaF = xoutDeltaF;
+
+ /* Spectrum Invert Truncate to smear out the peaks */
+ REAL8FFTPlan* fwdplan;
+ REAL8FFTPlan* revplan;
+ XLAL_CHECK( (fwdplan = XLALCreateForwardREAL8FFTPlan((UINT4) cfg->strainPaddedLength, 0)) != NULL, XLAL_ENOMEM);
+ XLAL_CHECK( (revplan = XLALCreateReverseREAL8FFTPlan((UINT4) cfg->strainPaddedLength, 0)) != NULL, XLAL_ENOMEM);
+
+ XLAL_CHECK( ( SpectrumInvertTruncate(psd, uvar->lowFreqCutoff,
+ (UINT4) cfg->strainPaddedLength,
+ (UINT4) uvar->psdDuration / Tseries->deltaT,
+ fwdplan,
+ revplan )) == XLAL_SUCCESS, XLAL_EFUNC);
+
+ /* Invert back again */
+ for(UINT4 index = 0; index < psd->data->length; index++) {
+ psd->data->data[index] = 1.0 / psd->data->data[index];
+ }
+ /* Low Frequency Cutoff */
+ UINT4 kmin = uvar->lowFreqCutoff / psd->deltaF;
+ for(UINT4 index = 0; index < kmin; index++) {
+ psd->data->data[index] = INFINITY;
+ }
+ /* High Frequency Cutoff */
+ if(XLALUserVarWasSet(&uvar->highFreqCutoff)) {
+ UINT4 kmax = uvar->highFreqCutoff / psd->deltaF;
+ for(UINT4 index = kmax; index < psd->data->length; index++) {
+ psd->data->data[index] = INFINITY;
+ }
+ }
+
+ XLALDestroyREAL8FFTPlan(fwdplan);
+ XLALDestroyREAL8FFTPlan(revplan);
+ *psdFseries = psd;
+
+#if TRACKMEMUSE
+ printf("Memory use after estimateWelchPSD:\n"); printmemuse();
+#endif
+ return XLAL_SUCCESS;
+}/* estimateWelchPSD() */
+
+/**
+**
+**/
+int copyReferencePSDFromWelchsPSD(UserVariables_t *uvar, ConfigVars_t *cfg, REAL8FrequencySeries **referencePSD ) {
+ REAL8 refEnd_freq = uvar->refFreqStart + uvar->refFreqBand;
+ UINT4 welchStart_index = (UINT4) ((uvar->refFreqStart - cfg->welchPSD->f0) / cfg->welchPSD->deltaF);
+ UINT4 welchEnd_index = (UINT4) ((refEnd_freq - cfg->welchPSD->f0) / cfg->welchPSD->deltaF);
+ UINT4 welchLength = welchEnd_index - welchStart_index;
+
+ //Chosen frequencies may not be actual sampling frequencies. Calculate correct frequencies.
+ uvar->refFreqStart = ((REAL8) welchStart_index) * cfg->welchPSD->deltaF + cfg->welchPSD->f0;
+ refEnd_freq = ((REAL8) welchEnd_index) * cfg->welchPSD->deltaF + cfg->welchPSD->f0;
+
+ REAL8 refDeltaF = 1.0 / (uvar->duration - 2 * uvar->padData); // Later used frequency sampling
+ UINT4 refLength = welchLength * cfg->welchPSD->deltaF / refDeltaF; // Length of refPSD with different sampling
+
+ REAL8FrequencySeries *refPSD = NULL;
+ XLAL_CHECK(( refPSD = XLALCreateREAL8FrequencySeries(
+ "Harmonic Mean PSD of whole obs run",
+ &cfg->strainPadded->epoch,
+ uvar->refFreqStart,
+ refDeltaF,
+ &lalDimensionlessUnit,
+ refLength )) != NULL, XLAL_ENOMEM);
+ /* Interpolate PSD */
+ REAL8Sequence *x_in, *y_in, *x_out, *y_out;
+
+ /* Populate x_in-array with frequencies freq-start to freq-end*/
+ XLAL_CHECK( (x_in = XLALCreateREAL8Sequence(welchLength) ) != NULL, XLAL_ENOMEM);
+ for(UINT4 index = 0; index < welchLength; index++) {
+ x_in->data[index] = uvar->refFreqStart + index * cfg->welchPSD->deltaF;
+ }
+
+ /* Populate y_in array with values for given frequencies */
+ XLAL_CHECK( (y_in = XLALCreateREAL8Sequence(welchLength) ) != NULL, XLAL_ENOMEM);
+ for(UINT4 index = 0; index < welchLength; index++) {
+ y_in->data[index] = cfg->welchPSD->data->data[index + welchStart_index];
+ }
+
+ /* Populate x_out with new frequencies ...*/
+ XLAL_CHECK( (x_out = XLALCreateREAL8Sequence(refLength) ) != NULL, XLAL_ENOMEM);
+ for(UINT4 index = 0; index < refLength; index++) {
+ x_out->data[index] = uvar->refFreqStart + index * refDeltaF;
+ }
+
+ y_out = refPSD->data;
+ /* Do Interpolation.*/
+ const gsl_interp_type* gsl_interp_p = gsl_interp_linear;
+ XLAL_CHECK( XLALREAL8Interpolation(x_in, y_in, x_out, y_out, welchLength , gsl_interp_p) == XLAL_SUCCESS, XLAL_EFUNC);
+
+ XLALDestroyREAL8Sequence(x_in);
+ XLALDestroyREAL8Sequence(y_in);
+ XLALDestroyREAL8Sequence(x_out);
+
+ *referencePSD = refPSD;
+
+ return XLAL_SUCCESS;
+} /* copyReferencePSDFromWelchsPSD() */
+
+/**
+ * Copied from, changed to a different norm: https://lscsoft.docs.ligo.org/lalsuite/lal/_average_spectrum_8c_source.html
+ */
+COMPLEX16FrequencySeries*
+WhitenCOMPLEX16FrequencySeries(UserVariables_t *uvar, REAL8 strainSampleRate, COMPLEX16FrequencySeries *fseries, const REAL8FrequencySeries *psd)
+{
+ COMPLEX16 *fdata = fseries->data->data;
+ REAL8 *pdata = psd->data->data;
+ double norm;
+ if(XLALUserVarWasSet ( &uvar->highFreqCutoff )) {
+ norm = uvar->highFreqCutoff - uvar->lowFreqCutoff;
+ } else {
+ norm = strainSampleRate/2 - uvar->lowFreqCutoff;
+ }
+ LALUnit unit; /* scratch space */
+ unsigned i; /* fseries index */
+ unsigned j; /* psd index */
+
+ if((psd->deltaF != fseries->deltaF) ||
+ (fseries->f0 < psd->f0)) {
+ /* resolution mismatch, or PSD does not span fseries at low end */
+ printf("psd->deltaF: %.16f vs fseries->deltaF: %.16f \n psd->f0: %.16f vs fseries->f0: %.16f \n",
+ psd->deltaF, fseries->deltaF, fseries->f0, psd->f0);
+ XLAL_ERROR_NULL(XLAL_EINVAL);
+ }
+ j = (fseries->f0 - psd->f0) / psd->deltaF;
+ if(j * psd->deltaF + psd->f0 != fseries->f0)
+ /* fseries does not start on an integer PSD sample */
+ XLAL_ERROR_NULL(XLAL_EINVAL);
+ if(j + fseries->data->length > psd->data->length)
+ /* PSD does not span fseries at high end */
+ XLAL_ERROR_NULL(XLAL_EINVAL);
+
+ for(i = 0; i < fseries->data->length; i++, j++)
+ {
+ if(pdata[j])
+ fdata[i] *= 1 / sqrt(norm * pdata[j]);
+ else
+ /* PSD has a 0 in it, treat as a zero in the filter */
+ fdata[i] = 0.0;
+ }
+
+ /* zero the DC and Nyquist components for safety */
+ if(fseries->data->length)
+ {
+ if(fseries->f0 == 0)
+ fdata[0] = 0.0;
+ fdata[fseries->data->length - 1] = 0.0;
+ }
+
+ /* update the units of fseries. norm has units of Hz */
+ XLALUnitDivide(&unit, &psd->sampleUnits, &lalHertzUnit);
+ XLALUnitSqrt(&unit, &unit);
+ XLALUnitDivide(&fseries->sampleUnits, &fseries->sampleUnits, &unit);
+
+ return fseries;
+} /* WhitenComplex16FrequencySeries() */
+
+/**
+ * Strain --FFT--> Strain Tilde --Whitening--> Whitened Strain Tilde
+ */
+int fftAndWhitenStrain(UserVariables_t *uvar, REAL8TimeSeries *Tseries, COMPLEX16FrequencySeries **CFseries, REAL8FrequencySeries *psdFseries) {
+ XLAL_CHECK( ( (*CFseries) = XLALCreateCOMPLEX16FrequencySeries("StrainTilde",
+ &Tseries->epoch, 0, (REAL8) 1.0 / (Tseries->data->length * Tseries->deltaT),
+ &lalStrainUnit, (UINT4) (Tseries->data->length / 2 + 1) ) ) != NULL, XLAL_ENOMEM);
+ REAL8FFTPlan* plan;
+ XLAL_CHECK( (plan = XLALCreateForwardREAL8FFTPlan(Tseries->data->length, 0)) != NULL, XLAL_ENOMEM);
+ XLAL_CHECK( XLALREAL8ForwardFFT((*CFseries)->data, Tseries->data, plan) == XLAL_SUCCESS, XLAL_EFUNC);
+
+ XLALDestroyREAL8FFTPlan(plan);
+
+ //Norm with 1/sample rate
+ for(UINT4 indexx = 0; indexx < (*CFseries)->data->length; indexx++) {
+ (*CFseries)->data->data[indexx] *= Tseries->deltaT;
+ }
+
+ //Whitening: XLAL Whitening function is using the wrong norm, take this whitening...
+ XLAL_CHECK( ((*CFseries) = WhitenCOMPLEX16FrequencySeries(uvar, 1.0 / Tseries->deltaT,
+ (*CFseries), psdFseries)) != NULL, XLAL_EFUNC);
+
+ XLALDestroyREAL8FrequencySeries(psdFseries);
+
+#if TRACKMEMUSE
+ printf("Memory use after fftAndWhiten is:\n"); printmemuse();
+#endif
+ return XLAL_SUCCESS;
+} /* fftAndWhitenStrain() */
+
+/**
+ * Inverse Fourier Transform and take the absolute: |Ifft(StrainTilde)|=Magnitude
+ */
+int ifft(REAL8TimeSeries **Tseries, COMPLEX16FrequencySeries *CFseries) {
+ REAL8FFTPlan *plan;
+ XLAL_CHECK( (plan = XLALCreateReverseREAL8FFTPlan((*Tseries)->data->length, 0)) != NULL, XLAL_ENOMEM);
+ XLAL_CHECK( (XLALREAL8FreqTimeFFT((*Tseries), CFseries, plan)) == XLAL_SUCCESS, XLAL_EFUNC);
+
+ XLALDestroyREAL8FFTPlan(plan);
+ XLALDestroyCOMPLEX16FrequencySeries(CFseries);
+
+#if TRACKMEMUSE
+ printf("Memory use after iffAndAbs is:\n"); printmemuse();
+#endif
+ return XLAL_SUCCESS;
+} /* ifft() */
+
+/**
+ * Given |IFFT(Whitened Strain Tilde)|, calculate the corresponding gate parameters
+ */
+int findGlitchesAndDuration(UserVariables_t *uvar, ConfigVars_t *cfg, REAL8TimeSeries *magnitude, REAL8TimeSeries *Tseries, GateParams_t *gtParams) {
+ printf("\n__Func: findGlitchesAndDuration: \n");
+
+ /* Prepare: Harmonic Mean and standard deviation of magnitude, then abs magnitude */
+ REAL8 mean = -1.1;
+ REAL8 std = -1.1;
+ REAL8 max = -1.1;
+ XLAL_CHECK( compMagnitudeStatistics(cfg, magnitude, &mean, &std, &max, 128) == XLAL_SUCCESS, XLAL_EFUNC);
+ for(UINT4 indexx = 0; indexx < magnitude->data->length; indexx++) {
+ magnitude->data->data[indexx] = fabs(magnitude->data->data[indexx]);
+ }
+ mean = 0;
+
+ /* Make a copy of the TimeSeries */
+ REAL8TimeSeries *workingCopy = NULL;
+ SFTtype *sft = NULL;
+ REAL8FrequencySeries *sftPSD = NULL;
+
+ char description[10];
+ strcpy(description, "Reference");
+
+ /* Iterate until SFT ASD is deemed good enough */
+ REAL8 threshFactor = -1;
+ REAL8 secDead = -1;
+ REAL8 sftRatio = DBL_MAX;
+ for(UINT2 iteration = 0; iteration <= uvar->maxIterations; iteration++) {
+ /* Destroy and make a new copy (Destroying NULL is allowed for [T/F]-Series') */
+ XLALDestroyREAL8TimeSeries(workingCopy);
+ /* Working copy: Remove zero padding but leave additional padding 2048s -> 1816s e.g.
+ Start after zeropad and add. padd, length is then duration - 2x add padd */
+ XLAL_CHECK( (workingCopy = XLALCutREAL8TimeSeries(Tseries,
+ cfg->padStart + uvar->padData / Tseries->deltaT,
+ uvar->duration / Tseries->deltaT - 2 * uvar->padData / Tseries->deltaT ) ) != NULL, XLAL_ENOMEM);
+ /* Set up gating constants for this iteration */
+ printf(" Iteration %d:\n", iteration);
+ cfg->threshHigh = uvar->upThreshold - uvar->iterReduceHighThresh * iteration;
+ threshFactor = uvar->sigmaFactor - uvar->iterReduceLowThreshFactor * iteration;
+ cfg->threshLow = mean + threshFactor * std;
+ secDead = uvar->secondsDeadTime + uvar->iterIncreaseDeadtime * iteration;
+ cfg->nDead = secDead / magnitude->deltaT;
+
+ /* Save Paramters to log */
+ printf(" New Parameters: High: %f Low: %f Dead: %f MaxMagnitude: %f\n", cfg->threshHigh, cfg->threshLow, secDead, max);
+ FILE *fp = fopen(cfg->logFileName, "a");
+ XLAL_CHECK( (fp = fopen(cfg->logFileName,"a")) != NULL, XLAL_EFUNC, "Could not open file %s", cfg->logFileName);
+ fprintf(fp, "\nIteration %d\nNew Parameters: HighThreshold: %f \nLowThreshold: %f \nSecondsDead: %f\n", iteration, cfg->threshHigh, cfg->threshLow, secDead);
+
+ if (cfg->threshHigh < cfg->threshLow) {
+ printf("Lower threshold is higher than higher threshold. Segment is not well behaving (and should not be used)? Setting lower threshold = high threshold!\n");
+ cfg->threshLow = cfg->threshHigh;
+ }
+ /* There are 4 cases:
+ 1. - sftRatio is good:
+ a. - maximum is below highThreshold: Assume result is good, no/good gates are found. Stop iteration.
+ b. - maximum is above highThreshold: As above, but note that the maximum still exceeds the higher threshold. (If this is triggered often, the uvar->refPSDRatio may be too high.
+ 2. - sftRatio is bad:
+ a. - maximum is below highThreshold: Omit iteration. Go into next iteration with lower threshold(s).
+ b. - maximum is above highThreshold: GetGlitches, gateData, calculate sftRatio.
+ */
+
+ if( sftRatio < uvar->refPSDRatio ) {
+ printf("Break at iteration %d: Maximum is %f and HighThreshold is %f while sftRatio@%f < %f\n", iteration, max, cfg->threshHigh, sftRatio, uvar->refPSDRatio);
+ if( max < cfg->threshHigh ) {
+ fprintf(fp, "Break at iteration %d: Maximum@%f is lower than HighThreshold@%f while sftRatio@%f < %f\n", iteration, max, cfg->threshHigh, sftRatio, uvar->refPSDRatio);
+ } else {
+ fprintf(fp, "Break at iteration %d: Maximum@%f is HIGHER than HighThreshold@%f while sftRatio@%f < %f\n", iteration, max, cfg->threshHigh, sftRatio, uvar->refPSDRatio);
+ }
+ fclose(fp);
+ break;
+ } else {
+ if ( iteration == 0 || max > cfg->threshHigh) {
+ fprintf(fp, "Iteration %d: Maximum %f (HighThreshold %f) & sftRatio %f >= %f\n", iteration, max, cfg->threshHigh, (sftRatio == DBL_MAX) ? -0.0: sftRatio, uvar->refPSDRatio);
+ fclose(fp);
+ XLAL_CHECK( getGlitchDuration(cfg, magnitude, gtParams) == XLAL_SUCCESS, XLAL_EFUNC);
+ REAL8 refStart_freq = cfg->referencePSD->f0;
+ REAL8 refEnd_freq = ((REAL8) cfg->referencePSD->data->length) * cfg->referencePSD->deltaF + cfg->referencePSD->f0;
+ XLAL_CHECK(( sftPSD = XLALCreateREAL8FrequencySeries("PSD of SFT", &workingCopy->epoch, refStart_freq,
+ 1.0 / (workingCopy->data->length * workingCopy->deltaT), &lalDimensionlessUnit, cfg->referencePSD->data->length )) != NULL, XLAL_ENOMEM);
+
+ XLAL_CHECK( gateData(uvar, workingCopy, gtParams) == XLAL_SUCCESS, XLAL_EFUNC);
+ XLAL_CHECK( callMakeSFT(uvar, cfg, refStart_freq, refEnd_freq - refStart_freq, workingCopy, &sft, 0, description, uvar->tukeyWindowRadius) == XLAL_SUCCESS, XLAL_EFUNC);
+ XLAL_CHECK( callComputePSD(refStart_freq, refEnd_freq - refStart_freq, 101, sft, &sftPSD) == XLAL_SUCCESS, XLAL_EFUNC);
+ XLALDestroyCOMPLEX8FrequencySeries(sft);
+
+ if (uvar->outputIntermediateSFT && iteration < 10) {
+ SFTtype *intermediatesft = NULL;
+ REAL8FrequencySeries *intermediatesftPSD = NULL;
+ char descrIntermediate[100];
+ snprintf(descrIntermediate, sizeof descrIntermediate, "%s%d", "Iteration", iteration);
+
+ XLAL_CHECK(( intermediatesftPSD = XLALCreateREAL8FrequencySeries("PSD of intermediate SFT", &workingCopy->epoch, uvar->SFTFreqStart,
+ 1.0 / (workingCopy->data->length * workingCopy->deltaT), &lalDimensionlessUnit, uvar->SFTFreqBand * workingCopy->data->length * workingCopy->deltaT )) != NULL, XLAL_ENOMEM);
+
+ XLAL_CHECK( callMakeSFT(uvar, cfg, uvar->SFTFreqStart, uvar->SFTFreqBand, workingCopy, &intermediatesft, 1, descrIntermediate, uvar->tukeyWindowRadius) == XLAL_SUCCESS, XLAL_EFUNC);
+ XLAL_CHECK( callComputePSD(uvar->SFTFreqStart, uvar->SFTFreqBand, 0, intermediatesft, &intermediatesftPSD) == XLAL_SUCCESS, XLAL_EFUNC);
+
+ snprintf(descrIntermediate, sizeof descrIntermediate, "%s%s%s%d%s", uvar->outDir, "/C_Intermediate_SFT_PSD", *cfg->identifier->data, iteration, ".txt");
+ fp = fopen(descrIntermediate,"wb");
+ XLAL_CHECK( (fp = fopen(descrIntermediate,"wb")) != NULL, XLAL_EFUNC, "Could not open file %s", descrIntermediate);
+ XLALWriteREAL8FrequencySeries2fp(fp, intermediatesftPSD);
+ fclose(fp);
+ XLALDestroyREAL8FrequencySeries(intermediatesftPSD);
+ }
+ XLAL_CHECK( compareHMeanToThis(cfg, sftPSD, &sftRatio) == XLAL_SUCCESS, XLAL_EFUNC);
+ XLALDestroyREAL8FrequencySeries(sftPSD);
+ } else {
+ printf("Omitting iteration %d: Maximum@%f is lower than HighThreshold@%f while sftRatio@%f is >= %f\n", iteration, max, cfg->threshHigh, sftRatio, uvar->refPSDRatio);
+ fprintf(fp, "Omitting iteration %d: Maximum@%f is lower than HighThreshold@%f while sftRatio@%f is >= %f\n", iteration, max, cfg->threshHigh, sftRatio, uvar->refPSDRatio);
+ fclose(fp);
+ }
+ }
+ char logsftRatio[100];
+ snprintf(logsftRatio, sizeof logsftRatio, "%s%s%s%s", uvar->outDir, "/sftRatio.txt", *cfg->identifier->data, ".txt");
+
+ XLAL_CHECK( (fp = fopen(logsftRatio,"a")) != NULL, XLAL_EFUNC, "Could not open file %s", logsftRatio);
+ fprintf(fp, "%1.16f\n", sftRatio);
+ fclose(fp);
+ }
+ //These functions free all memory associated with a LAL time series. It is safe to pass NULL to these functions.
+ XLALDestroyREAL8TimeSeries(workingCopy);
+ XLALDestroyREAL8TimeSeries( cfg->magnitude);
+ XLALDestroyREAL8FrequencySeries( cfg->referencePSD);
+// XLALDestroyREAL8FrequencySeries(sftPSD);
+#if TRACKMEMUSE
+ printf("Memory use after findGlitchesAndDuration is:\n"); printmemuse();
+#endif
+ return XLAL_SUCCESS;
+} /* findGlitchesAndDuration() */
+
+/**
+ * Calculate harmonic mean of mean and harmonic mean of standard deviation over x segments of Timeseries
+ */
+int compMagnitudeStatistics(ConfigVars_t *cfg, REAL8TimeSeries *Tseries, REAL8 *mean, REAL8 *std, REAL8 *max, UINT2 segments) {
+ *mean = 0.0;
+ *std = 0.0;
+ *max = -1.0;
+ REAL8 x;
+ // Check that (not-padded) timeseries is divideable by segments.
+ UINT4 TsLength = Tseries->data->length;
+ XLAL_CHECK( TsLength % segments == 0, XLAL_ESIZE );
+ UINT4 segLength = TsLength / segments;
+
+ REAL8Vector *means;
+ REAL8Vector *stds;
+ XLAL_CHECK( ( means = XLALCreateREAL8Vector(segments)) != NULL, XLAL_ENOMEM);
+ XLAL_CHECK( ( stds = XLALCreateREAL8Vector(segments)) != NULL, XLAL_ENOMEM);
+ //Calculate mean and std for each segment
+ for(UINT4 segIndex = 0; segIndex < segments; segIndex++) {
+ means->data[segIndex] = 0;
+ stds->data[segIndex] = 0;
+ for(UINT4 i = 0; i < segLength; i++) {
+ x = Tseries->data->data[segIndex * segLength + i];
+ means->data[segIndex] += x;
+ if(fabs(x) > *max) { *max = fabs(x); } //Adjust maximum if neccessary
+ }
+ means->data[segIndex] /= (REAL8) segLength;
+
+ for(UINT4 i = 0; i < segLength; i++) {
+ x = Tseries->data->data[segIndex * segLength + i] - means->data[segIndex];
+ stds->data[segIndex] += x*x;
+ }
+ stds->data[segIndex] /= (REAL8) segLength;
+ stds->data[segIndex] = sqrt(stds->data[segIndex]);
+ }
+
+ //Calculate the harmonic mean of the mean and std
+ for(UINT4 segIndex = 0; segIndex < segments; segIndex++) {
+ *mean += 1.0 / means->data[segIndex];
+ *std += 1.0 / stds->data[segIndex];
+ }
+ *mean = segments / *mean;
+ *std = segments / *std;
+
+ /* Log */
+ FILE *fp = NULL;
+ XLAL_CHECK( (fp = fopen(cfg->logFileName,"a")) != NULL, XLAL_EFUNC, "Could not open file %s", cfg->logFileName);
+ printf("__Mean: %.6f Std: %.6f Max: %.6f\n", *mean, *std, *max);
+ fprintf(fp, "Mean: %.6f Std: %.6f Max: %.6f\n", *mean, *std, *max);
+ fclose(fp);
+
+ XLALDestroyREAL8Vector(means);
+ XLALDestroyREAL8Vector(stds);
+
+ return XLAL_SUCCESS;
+} /* calculateHarmonicMeanOfMeanAndStd() */
+
+/**
+ * Load the harmonic mean PSD computed from multiple other time segments
+ */
+int loadHMeanPSD(UserVariables_t *uvar, ConfigVars_t *cfg) {
+ FILE *fp = NULL;
+ XLAL_CHECK( (fp = fopen(uvar->PSDReferenceFile,"r")) != NULL, XLAL_EFUNC, "Could not open file %s", uvar->PSDReferenceFile);
+ //Get length of Reference, as well as start and end frequency
+ REAL8 throwaway, maxFrequency, minFrequency=0;
+ UINT4 PSDLength = 0;
+ while(EOF != fscanf(fp, "%lf %lf\n", &maxFrequency, &throwaway)) {
+ if(PSDLength==0) {
+ minFrequency = maxFrequency;
+ }
+ PSDLength++;
+ }
+ XLAL_CHECK(( cfg->referencePSD = XLALCreateREAL8FrequencySeries(
+ "Harmonic Mean PSD of whole obs run",
+ &cfg->strainPadded->epoch,
+ minFrequency,
+ 1.0 / (uvar->duration - 2 * uvar->padData),
+ &lalStrainUnit,
+ PSDLength
+ )) != NULL, XLAL_ENOMEM);
+ PSDLength=0;
+ rewind(fp); //Reset filepointer to start of file
+ while(EOF != fscanf(fp, "%lf %lf\n", &throwaway, &cfg->referencePSD->data->data[PSDLength++])) {
+//Do not delete // printf("%lf\n", cfg->referencePSD->data->data[counter-1]);
+ }
+ fclose(fp);
+ return XLAL_SUCCESS;
+} /* loadHMeanPSD() */
+
+/**
+ * Load file with start and end of chosen gates
+ */
+int loadAndMergeAdditionalGatesFile(UserVariables_t *uvar, GateParams_t *gtParams) {
+ printf("Found %d gates\n", gtParams->gateAmount);
+
+ GateParams_t XLAL_INIT_DECL(addGates);
+ addGates.gateAmount = 0;
+ XLAL_CHECK( (addGates.gateStart = XLALCreateREAL8Vector(0) ) != NULL, XLAL_ENOMEM);
+ XLAL_CHECK( (addGates.gateEnd = XLALCreateREAL8Vector(0) ) != NULL, XLAL_ENOMEM);
+ XLAL_CHECK( (addGates.gateDuration = XLALCreateREAL8Vector(0) ) != NULL, XLAL_ENOMEM);
+
+ FILE *fp = NULL;
+ XLAL_CHECK( (fp = fopen(uvar->additionalGateFile,"r")) != NULL, XLAL_EFUNC, "Could not open file %s", uvar->additionalGateFile);
+
+ REAL8 gtStart, gtEnd, gtPrevEnd = 0;
+ while(EOF != fscanf(fp, "%lf %lf\n", >Start, >End)) {
+ /* Gate end or start within this time */
+ if ( gtEnd > XLALGPSGetREAL8(&uvar->startTime) && gtStart < XLALGPSGetREAL8(&uvar->startTime) + uvar->duration ) {
+ XLAL_CHECK( gtEnd - gtStart >= 0.0 , XLAL_EINVAL, "Provided end before start, format for additional gates file is 'gateStart gateEnd', given: %lf %lf", gtStart, gtEnd );
+ XLAL_CHECK( TRUE || gtPrevEnd < gtStart, XLAL_EINVAL, "Last gates end ends after next gates start, they may not overlap and must be sorted in ascending order, given:\nX %lf\n%lf %lf", gtPrevEnd, gtStart, gtEnd);
+ gtPrevEnd = gtEnd;
+
+ XLAL_CHECK( (addGates.gateStart = XLALResizeREAL8Vector(addGates.gateStart, addGates.gateAmount+1) ) != NULL, XLAL_ENOMEM);
+ XLAL_CHECK( (addGates.gateEnd = XLALResizeREAL8Vector(addGates.gateEnd, addGates.gateAmount+1) ) != NULL, XLAL_ENOMEM);
+ XLAL_CHECK( (addGates.gateDuration = XLALResizeREAL8Vector(addGates.gateDuration, addGates.gateAmount+1) ) != NULL, XLAL_ENOMEM);
+ addGates.gateStart->data[addGates.gateAmount] = gtStart;
+ addGates.gateEnd->data[addGates.gateAmount] = gtEnd;
+ addGates.gateDuration->data[addGates.gateAmount] = gtEnd - gtStart;
+ addGates.gateAmount += 1;
+ }
+ }
+ printf("Additional %d gates\n", addGates.gateAmount);
+
+ fclose(fp);
+
+ /* Merge both lists, otherwise we may apply an unwanted double-taper */
+ /* There are no additional gates, nothing to do but return, gtParams still valid: */
+ if (addGates.gateAmount == 0) {
+ XLALDestroyREAL8Vector(addGates.gateStart);
+ XLALDestroyREAL8Vector(addGates.gateEnd);
+ XLALDestroyREAL8Vector(addGates.gateDuration);
+ return XLAL_SUCCESS;
+ } else if(gtParams->gateAmount == 0) {
+ //Gatestrain found no gates, just copy additional into gtParams
+ XLAL_CHECK( (gtParams->gateStart = XLALResizeREAL8Vector(gtParams->gateStart, addGates.gateAmount ) ) != NULL, XLAL_ENOMEM);
+ XLAL_CHECK( (gtParams->gateEnd = XLALResizeREAL8Vector(gtParams->gateEnd, addGates.gateAmount ) ) != NULL, XLAL_ENOMEM);
+ XLAL_CHECK( (gtParams->gateDuration = XLALResizeREAL8Vector(gtParams->gateDuration, addGates.gateAmount) ) != NULL, XLAL_ENOMEM);
+ memcpy(gtParams->gateStart->data, addGates.gateStart->data, sizeof(addGates.gateStart->data[0]) * addGates.gateAmount);
+ memcpy(gtParams->gateEnd->data, addGates.gateEnd->data, sizeof(addGates.gateEnd->data[0]) * addGates.gateAmount);
+ memcpy(gtParams->gateDuration->data, addGates.gateDuration->data, sizeof(addGates.gateDuration->data[0]) * addGates.gateAmount);
+ gtParams->gateAmount = addGates.gateAmount;
+ XLALDestroyREAL8Vector(addGates.gateStart);
+ XLALDestroyREAL8Vector(addGates.gateEnd);
+ XLALDestroyREAL8Vector(addGates.gateDuration);
+ return XLAL_SUCCESS;
+ } else {
+ /* There are at least 1 gates in each */
+
+ /* Create boolean vector over time with 1/time-baseline spacing.
+ * Iterate over gates and set 1 everywhere a gate is.
+ * Afterwards iterate over vector and open gate everywhere 0->1 and close when 1->0
+ * Bonus? Taper overlap?
+ */
+ REAL8 tStart = XLALGPSGetREAL8(&uvar->startTime);
+ REAL8 tEnd = XLALGPSGetREAL8(&uvar->startTime) + uvar->duration;
+ REAL8 dt = 1.0 / uvar->Tsft;
+ UINT4 NSamples = (tEnd - tStart) / dt;
+ UINT4 NTaper = uvar->autogatingTaper / dt;
+ BOOLEAN * timeZeroedOut = XLALCalloc( NSamples, sizeof(BOOLEAN));
+
+ UINT4 NStart, NEnd;
+
+ for (UINT4 i = 0; i < addGates.gateAmount; i++) {
+ //Get Start and End in units of samples
+ NStart = fmax(0, floor(( addGates.gateStart->data[i] - tStart) / dt));
+ NEnd = fmin(NSamples, ceil(( addGates.gateEnd->data[i] - tStart)/ dt));
+ for (UINT4 j = NStart; j < NEnd; j++) {
+ timeZeroedOut[j] = TRUE;
+ }
+ }
+ for (UINT4 i = 0; i < gtParams->gateAmount; i++) {
+ //Get Start and End in units of samples
+ NStart = fmax(0, floor(( gtParams->gateStart->data[i] - tStart) / dt));
+ NEnd = fmin(NSamples, ceil(( gtParams->gateEnd->data[i] - tStart)/ dt));
+ for (UINT4 j = NStart; j < NEnd; j++) {
+ timeZeroedOut[j] = TRUE;
+ }
+ }
+
+ //Empty gtParams - set to zero is enough, gets resized next
+ gtParams->gateAmount = 0;
+ BOOLEAN previous = FALSE;
+ BOOLEAN next;
+ for (UINT4 j = 0; j < NSamples; j++) {
+ next = timeZeroedOut[j];
+ if (previous == next) {
+ //empty
+ } else if(previous == FALSE && next == TRUE) {
+ //Open gate
+ XLAL_CHECK( (gtParams->gateStart = XLALResizeREAL8Vector(gtParams->gateStart, gtParams->gateAmount+1) ) != NULL, XLAL_ENOMEM);
+ XLAL_CHECK( (gtParams->gateEnd = XLALResizeREAL8Vector(gtParams->gateEnd, gtParams->gateAmount+1) ) != NULL, XLAL_ENOMEM);
+ XLAL_CHECK( (gtParams->gateDuration = XLALResizeREAL8Vector(gtParams->gateDuration, gtParams->gateAmount+1) ) != NULL, XLAL_ENOMEM);
+ gtParams->gateStart->data[gtParams->gateAmount] = j * dt + tStart;
+
+ } else if (previous == TRUE && next == FALSE) {
+ BOOLEAN breakflag = FALSE;
+ UINT4 i = j; //Save current index
+ for (UINT4 k = 0; k < NTaper; k++) { //Spool forward and see if there is a next-true value inside the taper distance
+ if(i+k+1 >= NSamples) {break;} //Don't exceed bounds
+ if (timeZeroedOut[i+k+1] == TRUE) { //Check next value(s) for true within taper
+ j=i+k; //Advance j to 1 step before it's true again.
+ breakflag = TRUE; //Set breakflag, we don't want to close the gate
+ }
+ }
+ if (breakflag == FALSE) {
+ //Close gate
+ gtParams->gateEnd->data[gtParams->gateAmount] = j * dt + tStart;
+ gtParams->gateDuration->data[gtParams->gateAmount] = gtParams->gateEnd->data[gtParams->gateAmount] - gtParams->gateStart->data[gtParams->gateAmount];
+ gtParams->gateAmount+=1;
+ } else {
+ continue; //Continue, i.e. don't hit previous = next since next is FALSE now but we still have the gate open. Next j is true again.
+ }
+ }
+ previous = next;
+ }
+ //It may be that the last sample should be zeroed out, so we have to close the last gate manually.
+ if (previous == TRUE) {
+ gtParams->gateEnd->data[gtParams->gateAmount] = tEnd;
+ gtParams->gateDuration->data[gtParams->gateAmount] = gtParams->gateEnd->data[gtParams->gateAmount] - gtParams->gateStart->data[gtParams->gateAmount];
+ gtParams->gateAmount+=1;
+ }
+
+ XLALFree(timeZeroedOut);
+ XLALDestroyREAL8Vector(addGates.gateStart);
+ XLALDestroyREAL8Vector(addGates.gateEnd);
+ XLALDestroyREAL8Vector(addGates.gateDuration);
+ REAL8 totalDuration = 0;
+ for (UINT4 i = 0; i < gtParams->gateAmount; i++) {
+ totalDuration += gtParams->gateDuration->data[i];
+ }
+ printf("Merged %d gates with total duration %.2f\n", gtParams->gateAmount, totalDuration);
+ return XLAL_SUCCESS;
+ }
+} /* loadAndMergeAdditionalGatesFile */
+
+/**
+ * Load file with start and end of chosen gates
+ */
+int gateWithExternalGatesFile(UserVariables_t *uvar, ConfigVars_t *cfg, REAL8TimeSeries *Tseries) {
+ /* Reset gtParams */
+ GateParams_t XLAL_INIT_DECL(externalGates);
+ externalGates.gateAmount = 0;
+ XLALDestroyREAL8Vector(externalGates.gateStart);
+ XLALDestroyREAL8Vector(externalGates.gateEnd);
+ XLALDestroyREAL8Vector(externalGates.gateDuration);
+ XLAL_CHECK( (externalGates.gateStart = XLALCreateREAL8Vector(0) ) != NULL, XLAL_ENOMEM);
+ XLAL_CHECK( (externalGates.gateEnd = XLALCreateREAL8Vector(0) ) != NULL, XLAL_ENOMEM);
+ XLAL_CHECK( (externalGates.gateDuration = XLALCreateREAL8Vector(0) ) != NULL, XLAL_ENOMEM);
+
+ FILE *fp = NULL;
+ XLAL_CHECK( (fp = fopen(uvar->externalGateFile,"r")) != NULL, XLAL_EFUNC, "Could not open file %s", uvar->externalGateFile);
+
+ REAL8 gtStart, gtEnd;
+ while(EOF != fscanf(fp, "%lf %lf\n", >Start, >End)) {
+ /* Gate end or start within this time */
+ if ( gtEnd > XLALGPSGetREAL8(&uvar->startTime) && gtStart < XLALGPSGetREAL8(&uvar->startTime) + uvar->duration ) {
+ XLAL_CHECK( gtEnd - gtStart >= 0.0 , XLAL_EINVAL);
+
+ XLAL_CHECK( (externalGates.gateStart = XLALResizeREAL8Vector(externalGates.gateStart, externalGates.gateAmount+1) ) != NULL, XLAL_ENOMEM);
+ XLAL_CHECK( (externalGates.gateEnd = XLALResizeREAL8Vector(externalGates.gateEnd, externalGates.gateAmount+1) ) != NULL, XLAL_ENOMEM);
+ XLAL_CHECK( (externalGates.gateDuration = XLALResizeREAL8Vector(externalGates.gateDuration, externalGates.gateAmount+1) ) != NULL, XLAL_ENOMEM);
+ externalGates.gateStart->data[externalGates.gateAmount] = gtStart;
+ externalGates.gateEnd->data[externalGates.gateAmount] = gtEnd;
+ externalGates.gateDuration->data[externalGates.gateAmount] = gtEnd - gtStart;
+ externalGates.gateAmount += 1;
+ }
+ }
+ fclose(fp);
+
+ REAL8TimeSeries *strainCopy = NULL;
+
+ //Make a copy of the input strain
+ XLAL_CHECK( ( strainCopy = XLALCutREAL8TimeSeries(Tseries, 0, uvar->duration / Tseries->deltaT )) != NULL, XLAL_ENOMEM);
+
+ /* Gate data */
+ XLAL_CHECK( gateData(uvar, strainCopy, &externalGates) == XLAL_SUCCESS, XLAL_EFUNC);
+
+ if (uvar->outputHighpassedStrain) {
+ char buffer[100];
+ snprintf(buffer, sizeof buffer, "%s%s%s%s", uvar->outDir, "/C_B4SecHighpassExternalStrain", *cfg->identifier->data, ".txt");
+ XLAL_CHECK( (fp = fopen(buffer,"w")) != NULL, XLAL_EFUNC, "Could not open file %s", buffer);
+ XLALWriteREAL8TimeSeries2fp(fp, strainCopy);
+ fclose(fp);
+ }
+
+ /* Highpass again */
+ XLAL_CHECK(HighPass(cfg, strainCopy, uvar->Highpass) == XLAL_SUCCESS, XLAL_EFUNC);
+
+ /* Remove Padding */
+ XLAL_CHECK( (strainCopy = XLALResizeREAL8TimeSeries(strainCopy,
+ uvar->padData / strainCopy->deltaT, (uvar->duration - 2*uvar->padData) / strainCopy->deltaT) ) != NULL, XLAL_ENOMEM);
+
+ SFTtype *sft = NULL;
+ char description[32];
+ snprintf(description, 32, "ExternalGateFile_%s", cfg->SFTTukeyWindowNaming);
+
+ /* Make the SFT */
+ XLAL_CHECK( callMakeSFT(uvar, cfg, uvar->SFTFreqStart, uvar->SFTFreqBand, strainCopy, &sft, 1, description, uvar->tukeyWindowRadius) == XLAL_SUCCESS, XLAL_EFUNC);
+
+ XLALDestroyREAL8TimeSeries(strainCopy);
+ XLALDestroyCOMPLEX8FrequencySeries(sft);
+ XLALDestroyREAL8Vector(externalGates.gateStart);
+ XLALDestroyREAL8Vector(externalGates.gateEnd);
+ XLALDestroyREAL8Vector(externalGates.gateDuration);
+
+ return XLAL_SUCCESS;
+} /* gateWithExternalGatesFile() */
+
+
+/**
+ * Given a magnitude and thresholds, find values exceeding the high threshold and adjacent values
+ * exceeding the lower values to get a duration for a glitch.
+ * Gating constants are inside the ConfigVars_t *cfg variable
+ */
+int getGlitchDuration(ConfigVars_t *cfg, REAL8TimeSeries *magnitude, GateParams_t *gtParams) {
+ UINT4 leftCounter = 0;
+ UINT4 rightCounter = 0;
+ UINT4 curIndex = -1;
+ REAL8 timeStart = XLALGPSGetREAL8(&magnitude->epoch);
+
+ /* Reset gtParams */
+ gtParams->gateAmount = 0;
+ XLALDestroyREAL8Vector(gtParams->gateStart);
+ XLALDestroyREAL8Vector(gtParams->gateEnd);
+ XLALDestroyREAL8Vector(gtParams->gateDuration);
+ XLAL_CHECK( (gtParams->gateStart = XLALCreateREAL8Vector(0) ) != NULL, XLAL_ENOMEM);
+ XLAL_CHECK( (gtParams->gateEnd = XLALCreateREAL8Vector(0) ) != NULL, XLAL_ENOMEM);
+ XLAL_CHECK( (gtParams->gateDuration = XLALCreateREAL8Vector(0) ) != NULL, XLAL_ENOMEM);
+
+ //Iterate from 'left' to 'right' through the magnitude
+ for(UINT4 gtIndex = 0; gtIndex < magnitude->data->length; gtIndex++) {
+ if(magnitude->data->data[gtIndex] > cfg->threshHigh) {
+ leftCounter=1;
+ rightCounter=1;
+ curIndex = gtIndex;
+ while(leftCounter < cfg->nDead) {
+ if(curIndex < leftCounter) { //Same as curIndex - leftCounter < 0, but no unsigned error
+ break;
+ } else if(magnitude->data->data[curIndex - leftCounter] > cfg->threshLow) {
+ /* Going left, a value exceeding ThreshLow was detected. */
+ curIndex = curIndex - leftCounter;
+ leftCounter = 1;
+ } else {
+ /* If no value exceeding ThreshLow was detected, progress until leftCounter > nDead */
+ leftCounter++;
+ }
+ }
+ XLAL_CHECK( (gtParams->gateStart = XLALResizeREAL8Vector(gtParams->gateStart, gtParams->gateAmount+1) ) != NULL, XLAL_ENOMEM);
+ gtParams->gateStart->data[gtParams->gateAmount] = timeStart + curIndex * magnitude->deltaT;
+ printf(" Gate Start (GPS Time): %f\n", gtParams->gateStart->data[gtParams->gateAmount]);
+
+ curIndex = gtIndex;
+ while(rightCounter < cfg->nDead) {
+ if(curIndex + rightCounter >= magnitude->data->length) {
+ break;
+ } else if(magnitude->data->data[curIndex + rightCounter] > cfg->threshLow) {
+ curIndex = curIndex + rightCounter;
+ rightCounter = 1;
+ } else {
+ rightCounter++;
+ }
+ }
+ XLAL_CHECK( (gtParams->gateEnd = XLALResizeREAL8Vector(gtParams->gateEnd, gtParams->gateAmount+1) ) != NULL, XLAL_ENOMEM);
+ gtParams->gateEnd->data[gtParams->gateAmount] = timeStart + curIndex * magnitude->deltaT;
+ printf(" Gate End (GPS Time): %f\n", gtParams->gateEnd->data[gtParams->gateAmount]);
+
+ XLAL_CHECK( (gtParams->gateDuration = XLALResizeREAL8Vector(gtParams->gateDuration, gtParams->gateAmount+1) ) != NULL, XLAL_ENOMEM);
+ gtParams->gateDuration->data[gtParams->gateAmount] = gtParams->gateEnd->data[gtParams->gateAmount] - gtParams->gateStart->data[gtParams->gateAmount];
+ printf(" Gate duration: %f seconds\n\n", (gtParams->gateDuration->data[gtParams->gateAmount]));
+
+ gtParams->gateAmount += 1;
+ gtIndex = curIndex;
+ }
+ }
+ FILE *fp = NULL;
+ XLAL_CHECK( (fp = fopen(cfg->logFileName,"a")) != NULL, XLAL_EFUNC, "Could not open file %s", cfg->logFileName);
+ fprintf(fp, "Gate(s) found: %d\n", gtParams->gateAmount);
+ REAL8 gtDur=0;
+ for(UINT2 index = 0; index < gtParams->gateAmount; index++) {
+ fprintf(fp, "Start End Dur: %.16f %.16f %.16f\n", gtParams->gateStart->data[index], gtParams->gateEnd->data[index], gtParams->gateDuration->data[index]);
+ gtDur+=gtParams->gateDuration->data[index];
+ }
+ fclose(fp);
+ printf("Gate(s) found with duration: %d %.4fs\n", gtParams->gateAmount, gtDur);
+ return XLAL_SUCCESS;
+} /* getGlitchDuration() */
+
+/**
+ * Take a timeseries and gate with gateParams
+ */
+int gateData(UserVariables_t *uvar, REAL8TimeSeries *Tseries, GateParams_t *gtParams) {
+ REAL8Window *tukey;
+ REAL8 gateStartIndex = -1;
+ REAL8 gateDuration = 1;
+ REAL8 tukeyShapeParameter = -1;
+ UINT4 tseriesLength = Tseries->data->length;
+ /* The shape parameter e[0,1] sets what fraction of the window is spanned by the tapers.
+ β=0 yields the rectangle window, β=1 yields the Hann window.
+ We want padding only in taper, gate overall 0. If the tapered gate duration is t+d, the fraction which is
+ only taper is t/(t+d).For big d, the shape limits to 0, for small d the shape goes to 1.
+ */
+ for(UINT4 gtIndex = 0; gtIndex < gtParams->gateAmount; gtIndex++) {
+ tukey = NULL;
+ /* Duration of gate is initial duration + taper (in seconds) before and after */
+ gateDuration = (2 * uvar->autogatingTaper + gtParams->gateDuration->data[gtIndex]) / Tseries->deltaT;
+ tukeyShapeParameter = (2 * uvar->autogatingTaper / Tseries->deltaT) / gateDuration;
+
+ XLAL_CHECK( (tukey = XLALCreateTukeyREAL8Window(gateDuration, tukeyShapeParameter)) != NULL, XLAL_ENOMEM);
+ gateStartIndex = (gtParams->gateStart->data[gtIndex] - XLALGPSGetREAL8(&Tseries->epoch) - uvar->autogatingTaper) / Tseries->deltaT;
+
+ XLAL_CHECK(gateStartIndex > 0 || gateStartIndex < Tseries->data->length, XLAL_ERANGE);
+ //Apply (inverted) tukey to Tseries
+ for(UINT4 tukIndex = 0; tukIndex < tukey->data->length; tukIndex++) {
+ //Tukey is gate duration long: Therefore step thru tukey window from start and apply. Catch out of bounds which may happen due to taper.
+ if( (UINT4) gateStartIndex + tukIndex > 0 && (UINT4) gateStartIndex + tukIndex < tseriesLength ) {
+ Tseries->data->data[(UINT4) gateStartIndex + tukIndex] *= (1.0 - tukey->data->data[tukIndex]);
+ }
+ }
+ XLALDestroyREAL8Window(tukey);
+ }
+
+#if TRACKMEMUSE
+ printf("Memory use after gateData is:\n"); printmemuse();
+#endif
+ return XLAL_SUCCESS;
+} /* gateData() */
+
+/**
+ * Fields copied from MakeSFTs.c - https://lscsoft.docs.ligo.org/lalsuite/lalapps/_make_s_f_ts_8c_source.html
+ * Make SFT of the TimeSeries.
+ */
+int callMakeSFT(UserVariables_t *uvar, ConfigVars_t *cfg, REAL8 freqStart, REAL8 freqBand, REAL8TimeSeries *Tseries, SFTtype **sft, BOOLEAN write2File, char* description, REAL8 tukeyWindowRadius) {
+ SFTConfigVarsTag SFTVars;
+
+ /* Initialize default values */
+ SFTVars.FrCacheFile=NULL;
+ SFTVars.makeGPSDirs=0;
+ SFTVars.sftVersion=2;
+ SFTVars.commentField=NULL;
+ SFTVars.htdata = 0;
+ SFTVars.makeTmpFile = 0;
+// SFTVars.PSSCleaning = 0;
+// SFTVars.PSSCleanHPf = 100.0;
+// SFTVars.PSSCleanExt = 1;
+
+ /* Fill / Override with meaningful values */
+ SFTVars.HPf = uvar->Highpass;
+ SFTVars.T = (int) Tseries->data->length * Tseries->deltaT; //e.g. 1800
+ char stringT[12];
+ sprintf(stringT, "%d",SFTVars.T);
+ SFTVars.stringT = stringT;
+ SFTVars.GPSStart = XLALGPSGetREAL8(&(Tseries->epoch));
+ SFTVars.GPSEnd = XLALGPSGetREAL8(&(Tseries->epoch)) + Tseries->data->length * Tseries->deltaT;
+
+ SFTVars.ChannelName = uvar->inFrChannels->data[0];
+ char ifo[3];
+ strncpy(ifo, uvar->inFrChannels->data[0], 2);
+ ifo[2] = '\0';
+ SFTVars.IFO = ifo;
+ SFTVars.miscDesc = description;
+
+ //PSS Options....
+ SFTVars.windowOption = 1;
+ SFTVars.windowR = tukeyWindowRadius;
+ SFTVars.overlapFraction = 0.0;
+ SFTVars.useSingle = 0;
+ char frameStructType[] = "PROC_REAL8";
+ SFTVars.frameStructType = frameStructType;
+
+ /** Write SFT options */
+ SFTVars.SFTpath = uvar->outDir;
+ SFTVars.writeFile = write2File;
+
+ FILE *fp = NULL;
+ XLAL_CHECK( (fp = fopen(cfg->logFileName,"a")) != NULL, XLAL_EFUNC, "Could not open file %s", cfg->logFileName);
+ fprintf(fp, "\nRelevant makeSFT Parameters used: HPf %f - T %d - GPSStart %d - GPSEnd %d - sftVersion %d - ChannelName %s - miscDesc %s - IFO %s - freqStart %f - freqBand %f - windowR %f\n", SFTVars.HPf, SFTVars.T, SFTVars.GPSStart, SFTVars.GPSEnd, SFTVars.sftVersion, SFTVars.ChannelName, SFTVars.miscDesc, SFTVars.IFO, freqStart, freqBand, tukeyWindowRadius);
+ fclose(fp);
+
+ char commentField[4096];
+ snprintf(commentField, 4096, "\nRelevant makeSFT Parameters used: HPf %f - T %d - GPSStart %d - GPSEnd %d - sftVersion %d - ChannelName %s - miscDesc %s - IFO %s - freqStart %f - freqBand %f - windowR %f\n", SFTVars.HPf, SFTVars.T, SFTVars.GPSStart, SFTVars.GPSEnd, SFTVars.sftVersion, SFTVars.ChannelName, SFTVars.miscDesc, SFTVars.IFO, freqStart, freqBand, tukeyWindowRadius);
+ SFTVars.commentField = commentField;
+
+ XLAL_CHECK( makeSFT(SFTVars, Tseries, freqStart, freqBand, sft) == XLAL_SUCCESS, XLAL_EFUNC);
+ return XLAL_SUCCESS;
+} /* callMakeSFT() */
+
+/**
+ * Fields copied from ComputePSD.c - https://lscsoft.docs.ligo.org/lalsuite/lalapps/_compute_p_s_d_8c_source.html
+ * Compute PSD of the given SFT
+ */
+int callComputePSD(REAL8 freqStart, REAL8 freqBand, UINT2 blocksRngMed, SFTtype *sft, REAL8FrequencySeries **Fseries) {
+ PSDUserVarsTag PSDVars;
+ /* default: read all SFT bins */
+ PSDVars.outputPSD = NULL;
+ PSDVars.outputNormSFT = FALSE;
+ PSDVars.outFreqBinEnd = FALSE;
+
+ PSDVars.dumpMultiPSDVector = FALSE;
+
+ PSDVars.binSizeHz = 0.0;
+ PSDVars.binSize = 1;
+
+ PSDVars.binStep = 0.0;
+ PSDVars.binStep = 1;
+
+ PSDVars.blocksRngMed = blocksRngMed; //This does smoothing of PSD
+ /* Meaningful options: */
+ PSDVars.fStart = freqStart;
+ PSDVars.fBand = freqBand;
+
+ /* No output */
+ //char outputPSD[] = "./tempPSD.txt";
+ //PSDVars.outputPSD = outputPSD;
+
+ XLAL_CHECK( computePSD( PSDVars, sft, Fseries) == XLAL_SUCCESS, XLAL_EFUNC);
+
+ return XLAL_SUCCESS;
+} /* callComputePSD() */
+
+/**
+ * Compare this SFT PSD to the harmonic mean over multiple segments
+ */
+int compareHMeanToThis(ConfigVars_t *cfg, REAL8FrequencySeries *sftPSD, REAL8 *sftRatio) {
+ REAL8 mean = 0;
+ UINT4 vecLength = sftPSD->data->length;
+ for(UINT4 index = 0; index < vecLength -1 ; index++) {
+ mean += sftPSD->data->data[index] / cfg->referencePSD->data->data[index] / vecLength;
+ }
+ printf(" Mean Ratio: %f\n", mean);
+ mean = sqrt(mean);
+ printf(" Sqrt Mean Ratio: %f\n", mean);
+ *sftRatio = mean;
+
+ return XLAL_SUCCESS;
+} /* compareHMeanToThis() */
+
+/**
+ * Copied from: https://lscsoft.docs.ligo.org/lalsuite/lalapps/makefakedata__v5_8c_source.html
+ * Output timeseries into a new *.gwf file
+ */
+int outputTimeSeriesToFile(UserVariables_t *uvar, ConfigVars_t *cfg, REAL8TimeSeries *strain){
+#ifdef HAVE_LIBLALFRAME
+ size_t len;
+ printf("Writing Timeseries to file\n");
+ if ( uvar->outDir != NULL) {
+ XLAL_CHECK ( XLALCheckValidDescriptionField ( uvar->outLabel ) == XLAL_SUCCESS, XLAL_EFUNC );
+ len = strlen(uvar->outDir) + strlen(uvar->outLabel) + 100;
+ char *fname;
+
+ char *hist = XLALUserVarGetLog (UVAR_LOGFMT_CMDLINE);
+ /* Not used for single timeseries
+ if ( XLALUserVarWasSet ( &uvar->outFrChannels ) ) {
+ XLAL_CHECK ( uvar->outFrChannels->length == mTseries->length, XLAL_EINVAL, "--outFrChannels: number of channel names (%d) must agree with number of IFOs (%d)\n",
+ uvar->outFrChannels->length, mTseries->length );
+ }
+ */
+ for ( UINT4 X=0; X < 1 /*mTseries->length*/; X ++ ) {
+ // REAL8TimeSeries *Tseries = mTseries->data[X];
+ REAL8TimeSeries *Tseries = strain;
+
+ /* use standard frame output filename format */
+ char IFO[2] = { Tseries->name[0], Tseries->name[1] };
+ LIGOTimeGPS startTimeGPS = Tseries->epoch;
+ REAL8 duration = Tseries->data->length * Tseries->deltaT;
+ XLAL_CHECK ( (fname = LALCalloc (1, len )) != NULL, XLAL_ENOMEM );
+ size_t written = snprintf ( fname, len, "%s/%c-%c%c_%s-%d-%d.gwf",
+ uvar->outDir, IFO[0], IFO[0], IFO[1], uvar->outLabel, startTimeGPS.gpsSeconds, (int)duration );
+ XLAL_CHECK ( written < len, XLAL_ESIZE, "Frame-filename exceeds expected maximal length (%zu): '%s'\n", len, fname );
+
+ /* define the output frame */
+ LALFrameH *outFrame;
+ XLAL_CHECK ( (outFrame = XLALFrameNew ( &startTimeGPS, duration, uvar->outLabel, 1, 0, 0 )) != NULL, XLAL_EFUNC );
+
+ /* add timeseries to the frame - make sure to change the timeseries name since this is used as the channel name */
+ char buffer[LALNameLength];
+ // if output frame channel names given, use those
+ if ( XLALUserVarWasSet ( &uvar->outFrChannels ) ) {
+ written = snprintf ( buffer, sizeof(buffer), "%s", uvar->outFrChannels->data[X] );
+ if ( buffer[2] == ':' ) { // check we got correct IFO association
+ XLAL_CHECK ( (buffer[0] == Tseries->name[0]) && (buffer[1] == Tseries->name[1]), XLAL_EINVAL,
+ "Possible IFO mismatch: outFrChannel[%d] = '%s', IFO = '%c%c': be careful about --outFrChannel ordering\n", X, buffer, Tseries->name[0], Tseries->name[1] );
+ } // if buffer[2]==':'
+ } else if ( XLALUserVarWasSet ( &uvar->inFrChannels ) ) { // otherwise: if input frame channel names given, use them for output, append "-<outLabel>"
+ written = snprintf ( buffer, sizeof(buffer), "%s-%s", uvar->inFrChannels->data[X], uvar->outLabel );
+ } else { // otherwise: fall back to <IFO>:<outLabel> channel name
+ written = snprintf ( buffer, sizeof(buffer), "%c%c:%s", Tseries->name[0], Tseries->name[1], uvar->outLabel );
+ }
+ XLAL_CHECK ( written < LALNameLength, XLAL_ESIZE, "Output frame name exceeded max length (%d): '%s'\n", LALNameLength, buffer );
+ strcpy ( Tseries->name, buffer );
+
+ XLAL_CHECK ( (XLALFrameAddREAL8TimeSeriesProcData ( outFrame, Tseries ) == XLAL_SUCCESS ) , XLAL_EFUNC );
+
+ /* Here's where we add extra information into the frame - first we add the command line args used to generate it */
+ XLALFrameAddFrHistory ( outFrame, __FILE__, hist );
+
+ /* then we add the version string */
+ XLALFrameAddFrHistory ( outFrame, __FILE__, cfg->VCSInfoString );
+
+ /* output the frame to file - compression level 1 (higher values make no difference) */
+ XLAL_CHECK ( XLALFrameWrite ( outFrame, fname ) == XLAL_SUCCESS , XLAL_EFUNC );
+
+ /* free the frame, frame file name and history memory */
+ XLALFrameFree ( outFrame );
+ XLALFree ( fname );
+
+ } // for X < numDetectors
+ XLALFree ( hist );
+
+ }
+
+#if TRACKMEMUSE
+ printf("Memory use after outputTimeSeriesToFile is:\n"); printmemuse();
+#endif
+
+ return XLAL_SUCCESS;
+#else
+ printf("NEED LIBLALFRAME TO FUNCTION. RECOMPILE\n");
+ (void)(uvar);
+ (void)(cfg);
+ (void)(strain); /* suppress unused var warning/error */
+ return XLAL_SUCCESS;
+#endif
+
+} /* outputTimeSeriesToFile() */
+
+
+/**
+ * This routine frees up all the memory.
+ */
+int
+XLALFreeMem ( UserVariables_t *uvar, ConfigVars_t *cfg, GateParams_t *gtParams )
+{
+ XLAL_CHECK ( cfg != NULL, XLAL_EINVAL );
+ /* Free config-Variables and userInput stuff */
+ XLALFree(uvar->outLabel);
+ XLALDestroyUserVars();
+ XLALDestroyStringVector( cfg->identifier);
+ XLALDestroyStringVector( cfg->fileNames);
+ XLALFree ( cfg->VCSInfoString );
+ XLALDestroyREAL8TimeSeries( cfg->strainPadded);
+ XLALDestroyREAL8TimeSeries( cfg->originalStrain);
+ XLALDestroyREAL8Vector(gtParams->gateStart);
+ XLALDestroyREAL8Vector(gtParams->gateEnd);
+ XLALDestroyREAL8Vector(gtParams->gateDuration);
+
+ return XLAL_SUCCESS;
+} /* XLALFreeMem() */
+
+
+/**
+ * pyCBC gate data Method. Returns gate parameters which can be used in the gate data method. Writes a log file if gtAmount>0.
+ */
+int pycbcGateDataMethod(UserVariables_t *uvar, ConfigVars_t *cfg, REAL8TimeSeries *magnitude, GateParams_t *gtParams) {
+ /* Reset gtParams */
+ gtParams->gateAmount = 0;
+ XLALDestroyREAL8Vector(gtParams->gateStart);
+ XLALDestroyREAL8Vector(gtParams->gateEnd);
+ XLALDestroyREAL8Vector(gtParams->gateDuration);
+ XLAL_CHECK( (gtParams->gateStart = XLALCreateREAL8Vector(0) ) != NULL, XLAL_ENOMEM);
+
+ REAL8 timeStart = XLALGPSGetREAL8(&magnitude->epoch);
+
+ UINT4 curr_index = 0;
+ REAL8 curr_val = 0;
+ REAL8 cluster_window = 0.5 / magnitude->deltaT;
+ for(UINT4 index = 0; index < magnitude->data->length; index++) {
+ if( magnitude->data->data[index] > uvar->upThreshold) {
+ if( (index - curr_index) > cluster_window || gtParams->gateAmount == 0) {
+ XLAL_CHECK( (gtParams->gateStart = XLALResizeREAL8Vector(gtParams->gateStart, gtParams->gateAmount+1) ) != NULL, XLAL_ENOMEM);
+ gtParams->gateStart->data[gtParams->gateAmount] = timeStart + index * magnitude->deltaT;
+ printf("Gate Mid (GPS Time): %f\n", gtParams->gateStart->data[gtParams->gateAmount]);
+
+ gtParams->gateAmount += 1;
+
+ curr_index = index;
+ curr_val = magnitude->data->data[index];
+ } else if(magnitude->data->data[index] > curr_val) {
+ gtParams->gateStart->data[gtParams->gateAmount - 1] = timeStart + index * magnitude->deltaT;
+
+ curr_index = index;
+ curr_val = magnitude->data->data[index];
+ }
+ } //else do nothing
+ }
+ char fileName[100];
+ snprintf(fileName, 100, "%s%s%s%s", uvar->outDir, "/pyCBCGates", *cfg->identifier->data, ".txt");
+
+ FILE *fp = NULL;
+ XLAL_CHECK( (fp = fopen(fileName,"w")) != NULL, XLAL_EFUNC, "Could not open file %s", fileName);
+ XLAL_CHECK( (gtParams->gateEnd = XLALCreateREAL8Vector(gtParams->gateAmount) ) != NULL, XLAL_ENOMEM);
+ XLAL_CHECK( (gtParams->gateDuration = XLALCreateREAL8Vector(gtParams->gateAmount) ) != NULL, XLAL_ENOMEM);
+ for(UINT2 index = 0; index < gtParams->gateAmount; index++) {
+ fprintf(fp, "%.16f\n", gtParams->gateStart->data[index]);
+ gtParams->gateDuration->data[index] = uvar->pyCBCFixedDurationGate;
+ gtParams->gateStart->data[index] = gtParams->gateStart->data[index] - uvar->pyCBCFixedDurationGate / 2.0;
+ gtParams->gateEnd->data[index] = gtParams->gateStart->data[index] + uvar->pyCBCFixedDurationGate;
+ }
+ fclose(fp);
+
+ return XLAL_SUCCESS;
+}
+
+
+/**
+ * Return the smallest integer power of 2 larger than the argument
+ */
+UINT4 nextPowerOf2(UINT4 x) {
+ return (UINT4) pow ( 2, ceil ( log2 ( x ) ) );
+ }
+
+/**
+ * Copied from: https://lscsoft.docs.ligo.org/lalsuite/lalapps/makefakedata__v5_8c_source.html
+ * determine if the given filepath is an existing directory or not
+ */
+BOOLEAN
+is_directory ( const CHAR *fname )
+{
+ struct stat stat_buf;
+ if ( stat (fname, &stat_buf) )
+ return 0;
+ if ( ! S_ISDIR (stat_buf.st_mode) )
+ return 0;
+ else
+ return 1;
+} /* is_directory() */
+
+/**
+ * Copied and converted to REAL8 from: https://lscsoft.docs.ligo.org/lalsuite/lalapps/makefakedata__v5_8c_source.html
+ * Write a REAL8TimeSeries in a 2-column ASCII format (lines of "GPS_i data_i")
+ * into an open file 'fp'
+ */
+int
+XLALWriteREAL8TimeSeries2fp ( FILE *fp, const REAL8TimeSeries *TS )
+{
+ printf("Writing TimeSeries to file... \nSTART");
+ XLAL_CHECK ( fp != NULL, XLAL_EINVAL );
+ XLAL_CHECK ( TS != NULL, XLAL_EINVAL );
+ XLAL_CHECK ( (TS->data != NULL) && (TS->data->length >0) && (TS->data->data != NULL), XLAL_EINVAL );
+ REAL8 t0 = XLALGPSGetREAL8( &(TS->epoch) );
+ XLAL_CHECK ( xlalErrno == 0, XLAL_EFAULT, "XLALGPSGetREAL8(%d,%d) failed\n", TS->epoch.gpsSeconds, TS->epoch.gpsNanoSeconds );
+ for( UINT4 i = 0; i < TS->data->length; i++)
+ {
+ REAL8 ti = t0 + i * TS->deltaT;
+ XLAL_CHECK ( fprintf( fp, "%16.9f %1.16e\n", ti, TS->data->data[i] ) > 0, XLAL_EIO );
+ }
+ printf(" ... END\n");
+ return XLAL_SUCCESS;
+} /* XLALWriteREAL8TimeSeries2fp() */
+
+/**
+ * Copied and converted to REAL8 from: https://lscsoft.docs.ligo.org/lalsuite/lalapps/makefakedata__v5_8c_source.html
+ * Write a REAL8FrequencySeries in a 2-column ASCII format (lines of "freq_i data_i")
+ * into an open file 'fp'
+ */
+int
+XLALWriteREAL8FrequencySeries2fp ( FILE* fp, const REAL8FrequencySeries *FS )
+{
+ printf("Writing REAL8FrequencySeries to file... \nSTART");
+ XLAL_CHECK ( fp != NULL, XLAL_EINVAL );
+ XLAL_CHECK ( FS != NULL, XLAL_EINVAL );
+ XLAL_CHECK ( (FS->data != NULL) && (FS->data->length >0) && (FS->data->data != NULL), XLAL_EINVAL );
+ REAL8 f0 = FS->f0;
+ for( UINT4 i = 0; i < FS->data->length; i++)
+ {
+ REAL8 fi = f0 + i * FS->deltaF;
+ XLAL_CHECK ( fprintf( fp, "%16.9f %1.16e\n", fi, FS->data->data[i] ) > 0, XLAL_EIO );
+ }
+ printf(" ... END\n");
+ return XLAL_SUCCESS;
+} /* XLALWriteREAL8TimeSeries2fp() */
+
+/**
+ * Copied and changed to other norm from: https://lscsoft.docs.ligo.org/lalsuite/lal/_average_spectrum_8c_source.html
+ * This version norms by dividing thruough the length of the (padded) timeseries (in seconds) and applying a hann window
+*/
+int SpectrumInvertTruncate(
+ REAL8FrequencySeries *spectrum,
+ REAL8 lowfreq,
+ UINT4 seglen,
+ UINT4 trunclen,
+ REAL8FFTPlan *fwdplan,
+ REAL8FFTPlan *revplan
+ )
+{
+ UINT4 cut;
+ UINT4 k;
+
+ if ( ! spectrum )
+ XLAL_ERROR( XLAL_EFAULT );
+ if ( ! spectrum->data )
+ XLAL_ERROR( XLAL_EINVAL );
+ if ( spectrum->deltaF <= 0.0 )
+ XLAL_ERROR( XLAL_EINVAL );
+ if ( lowfreq < 0 )
+ XLAL_ERROR( XLAL_EINVAL );
+ if ( ! seglen || spectrum->data->length != seglen/2 + 1 )
+ XLAL_ERROR( XLAL_EBADLEN );
+ if ( trunclen && ! revplan )
+ XLAL_ERROR( XLAL_EFAULT );
+
+ cut = lowfreq / spectrum->deltaF;
+ if ( cut < 1 ) /* need to get rid of DC at least */
+ cut = 1;
+
+ if ( trunclen ) /* truncate while inverting */
+ {
+ COMPLEX16Vector *vtilde; /* frequency-domain vector workspace */
+ REAL8Vector *vector; /* time-domain vector workspace */
+ REAL8 normfac;
+
+ vtilde = XLALCreateCOMPLEX16Vector( seglen / 2 + 1 );
+ vector = XLALCreateREAL8Vector( seglen );
+
+ /* clear the low frequency components */
+ memset( vtilde->data, 0, cut * sizeof( *vtilde->data ) );
+
+ /* stick the root-inv-spectrum into the complex frequency-domain vector */
+ for ( k = cut; k < vtilde->length - 1; ++k )
+ {
+ if ( spectrum->data->data[k] )
+ vtilde->data[k] = 1.0 / sqrt( spectrum->data->data[k] ); /* purely real */
+ else
+ vtilde->data[k] = 0.0;
+ }
+ /* no Nyquist */
+ vtilde->data[vtilde->length - 1] = 0.0;
+
+ /* construct time-domain version of root-inv-spectrum */
+ XLALREAL8ReverseFFT( vector, vtilde, revplan );
+
+ /* now truncate it: zero time from trunclen/2 to length - trunclen/2 */
+ memset( vector->data + trunclen/2, 0,
+ (vector->length - trunclen) * sizeof( *vector->data ) );
+
+ /* Start changes from lalsuite function */
+ /* Norm by the length of the TimeSeries in seconds */
+ REAL8Window* hann;
+ XLAL_CHECK( (hann = XLALCreateHannREAL8Window((UINT4) trunclen) ) != NULL, XLAL_ENOMEM);
+
+ for( UINT4 i = 0; i < (UINT4) trunclen/2; i++) {
+ vector->data[i] *= hann->data->data[(UINT4) trunclen/2 + i];
+ vector->data[vector->length - i - 1] *= hann->data->data[(UINT4) trunclen/2 + i];
+ }
+ XLALDestroyREAL8Window(hann);
+ /* End changes from lalsuite function */
+
+ /* reconstruct spectrum */
+ XLAL_CHECK(XLALREAL8PowerSpectrum(spectrum->data, vector, fwdplan) == XLAL_SUCCESS, XLAL_EFUNC);
+
+ memset( spectrum->data->data, 0, cut * sizeof( *spectrum->data->data ) );
+ spectrum->data->data[spectrum->data->length - 1] = 0.0;
+
+ /* rescale spectrum: 0.5 to undo factor of two from REAL8PowerSpectrum */
+ /* seglen^2 for the reverse fft (squared because power spectrum) */
+ /* Note: cast seglen to real so that seglen*seglen is a real rather */
+ /* than a four-byte integer which might overflow... */
+ normfac = 0.5 / ( ((REAL8)seglen) * ((REAL8)seglen) );
+ for ( k = cut; k < spectrum->data->length - 1; ++k )
+ spectrum->data->data[k] *= normfac;
+
+ /* clear the low frequency components... again, and Nyquist too */
+ /* cleanup workspace */
+ XLALDestroyCOMPLEX16Vector( vtilde );
+ XLALDestroyREAL8Vector( vector );
+ }
+ else /* otherwise just invert the spectrum */
+ {
+ /* clear the low frequency components */
+ memset( spectrum->data->data, 0, cut * sizeof( *spectrum->data->data ) );
+
+ /* invert the high frequency (non-Nyquist) components */
+ for ( k = cut; k < spectrum->data->length - 1; ++k )
+ {
+ if ( spectrum->data->data[k] )
+ spectrum->data->data[k] = 1.0 / spectrum->data->data[k];
+ else
+ spectrum->data->data[k] = 0.0;
+ }
+
+ /* zero Nyquist */
+ spectrum->data->data[spectrum->data->length - 1] = 0.0;
+ }
+
+ /* now correct the units */
+ XLALUnitRaiseINT2( &spectrum->sampleUnits, &spectrum->sampleUnits, -1 );
+
+ return 0;
+}
diff --git a/lalapps/src/pulsar/MakeData/testMakeSFTs.sh b/lalapps/src/pulsar/MakeData/testMakeSFTs.sh
index 516039c38bea33d15b3183cb8d0ab22cfcdc3c2f..b5a4fda1d4220417139852f49c81bc5bfa1dcf91 100644
--- a/lalapps/src/pulsar/MakeData/testMakeSFTs.sh
+++ b/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