ComputeFStatistic_resamp.c

Go to the documentation of this file.

/*
 * Copyright (C) 2008 Pinkesh Patel, Xavier Siemens
 * Copyright (C) 2007 Reinhard Prix, Iraj Gholami , Pinkesh Patel, Xavier Siemens
 * Copyright (C) 2005, 2006 Reinhard Prix, Iraj Gholami
 * Copyright (C) 2004 Reinhard Prix
 * Copyright (C) 2002, 2003, 2004 M.A. Papa, X. Siemens, Y. Itoh
 *
 *  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
 */

/*********************************************************************************/
/** \author P.Patel, X.Siemens, R. Prix, I. Gholami, Y. Ioth,M. Papa
 * \file
 * \brief
 * Calculate the F-statistic for a given parameter-space of pulsar GW signals.
 * Implements the so-called "F-statistic" as introduced in \ref JKS98.
 *
 *********************************************************************************/
#include "config.h"

/* System includes */
#include <math.h>
#include <stdio.h>
#include <time.h>
#ifdef HAVE_UNISTD_H
#include <unistd.h>
#endif


int finite(double);

/* LAL-includes */
#include <lal/AVFactories.h>
#include <lal/LALInitBarycenter.h>
#include <lal/UserInput.h>
#include <lal/SFTfileIO.h>
#include <lal/ExtrapolatePulsarSpins.h>
#include <lal/FrequencySeries.h>

#include <lal/NormalizeSFTRngMed.h>
#include <lal/ComputeFstat.h>
#include <lal/LALHough.h>

#include <lal/LogPrintf.h>
#include <lal/DopplerFullScan.h>
#include <lal/ComplexFFT.h>
#include <lal/LALBarycenter.h>

#include <lalapps.h>
#include <lal/Window.h>
#include <fftw3.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_interp.h>

/* local includes */

#include "../HeapToplist.h"

RCSID("$Id: ComputeFStatistic_resamp.c,v 1.30 2008/08/21 07:44:03 ppatel Exp $");
NRCSID(TEMPORARY,"$Blah$");

/*---------- DEFINES ----------*/

#define MAXFILENAMELENGTH 256   /* Maximum # of characters of a SFT filename */

#define EPHEM_YEARS  "00-04"    /**< default range: override with --ephemYear */

#define TRUE (1==1)
#define FALSE (1==0)

/*----- SWITCHES -----*/
#define NUM_SPINS 4             /* number of spin-values to consider: {f, fdot, f2dot, ... } */

/*----- Error-codes -----*/
#define COMPUTEFSTATISTIC_ENULL         1
#define COMPUTEFSTATISTIC_ESYS          2
#define COMPUTEFSTATISTIC_EINPUT        3
#define COMPUTEFSTATISTIC_EMEM          4
#define COMPUTEFSTATISTIC_ENONULL       5
#define COMPUTEFSTATISTIC_EXLAL         6

#define COMPUTEFSTATISTIC_MSGENULL      "Arguments contained an unexpected null pointer"
#define COMPUTEFSTATISTIC_MSGESYS       "System call failed (probably file IO)"
#define COMPUTEFSTATISTIC_MSGEINPUT     "Invalid input"
#define COMPUTEFSTATISTIC_MSGEMEM       "Out of memory. Bad."
#define COMPUTEFSTATISTIC_MSGENONULL    "Output pointer is non-NULL"
#define COMPUTEFSTATISTIC_MSGEXLAL      "XLALFunction-call failed"

/*----- Macros -----*/

/* convert GPS-time to REAL8 */
#define GPS2REAL8(gps) (1.0 * (gps).gpsSeconds + 1.e-9 * (gps).gpsNanoSeconds )
#define SQ(x) ( (x) * (x) )

#define MYMAX(x,y) ( (x) > (y) ? (x) : (y) )
#define MYMIN(x,y) ( (x) < (y) ? (x) : (y) )

#define LAL_INT4_MAX 2147483647

/*---------- internal types ----------*/

/** What info do we want to store in our toplist? */
typedef struct {
  PulsarDopplerParams doppler;          /**< Doppler params of this 'candidate' */
  Fcomponents  Fstat;                   /**< the Fstat-value (plus Fa,Fb) for this candidate */
  CmplxAntennaPatternMatrix Mmunu;              /**< antenna-pattern matrix Mmunu = 0.5* Sinv*Tsft * [ Ad, Cd; Cd; Bd ] */
} FstatCandidate;


/** moving 'Scanline window' of candidates on the scan-line,
 * which is used to find local 1D maxima.
 */
typedef struct
{
  UINT4 length;
  FstatCandidate *window;               /**< array holding candidates */
  FstatCandidate *center;               /**< pointer to middle candidate in window */
} scanlineWindow_t;

/** Configuration settings required for and defining a coherent pulsar search.
 * These are 'pre-processed' settings, which have been derived from the user-input.
 */
typedef struct {
  REAL8 Tsft;                               /**< length of one SFT in seconds */
  LIGOTimeGPS refTime;                      /**< reference-time for pulsar-parameters in SBB frame */
  /* -------------------- Resampling -------------------- */
  REAL8 FFTFreqBand;                        /**< treated outside of DopplerScan for resampling-technique */
  /* ---------------------------------------------------- */
  DopplerRegion searchRegion;               /**< parameter-space region (at *internalRefTime*) to search over */
  DopplerFullScanState *scanState;          /**< current state of the Doppler-scan */
  PulsarDopplerParams stepSizes;            /**< user-preferences on Doppler-param step-sizes */
  EphemerisData *ephemeris;                 /**< ephemeris data (from LALInitBarycenter()) */
  MultiSFTVector *multiSFTs;                /**< multi-IFO SFT-vectors */
  MultiDetectorStateSeries *multiDetStates; /**< pos, vel and LMSTs for detector at times t_i */
  MultiNoiseWeights *multiNoiseWeights;     /**< normalized noise-weights of those SFTs */
  ComputeFParams CFparams;                  /**< parameters for Fstat (e.g Dterms, SSB-prec,...) */
  CHAR *logstring;                          /**< log containing max-info on this search setup */
  toplist_t* FstatToplist;                  /**< sorted 'toplist' of the NumCandidatesToKeep loudest candidates */
  scanlineWindow_t *scanlineWindow;         /**< moving window of candidates on scanline to find local maxima */
} ConfigVariables;

LALUnit empty_Unit;

/*---------- Global variables ----------*/
extern UINT4 vrbflg;            /**< defined in lalapps.c */

/* ----- Resampling Variables ----------*/
REAL8 uvar_F;
/*REAL8 uvar_Het;*/


/* ----- User-variables: can be set from config-file or command-line */
INT4 uvar_Dterms;
CHAR *uvar_IFO;
BOOLEAN uvar_SignalOnly;
BOOLEAN uvar_UseNoiseWeights;
REAL8 uvar_Freq;
REAL8 uvar_FreqBand;
/* REAL8 uvar_dFreq;   deactivated for resampling version */
REAL8 uvar_Alpha;
REAL8 uvar_dAlpha;
REAL8 uvar_AlphaBand;
REAL8 uvar_Delta;
REAL8 uvar_dDelta;
REAL8 uvar_DeltaBand;
/* 1st spindown */
REAL8 uvar_f1dot;
REAL8 uvar_df1dot;
REAL8 uvar_f1dotBand;
/* 2nd spindown */
REAL8 uvar_f2dot;
REAL8 uvar_df2dot;
REAL8 uvar_f2dotBand;
/* 3rd spindown */
REAL8 uvar_f3dot;
REAL8 uvar_df3dot;
REAL8 uvar_f3dotBand;
/* --- */
REAL8 uvar_TwoFthreshold;
CHAR *uvar_ephemDir;
CHAR *uvar_ephemYear;
INT4  uvar_gridType;
INT4  uvar_metricType;
BOOLEAN uvar_projectMetric;
REAL8 uvar_metricMismatch;
CHAR *uvar_skyRegion;
CHAR *uvar_DataFiles;
BOOLEAN uvar_help;
CHAR *uvar_outputLogfile;
CHAR *uvar_outputFstat;
CHAR *uvar_outputLoudest;

INT4 uvar_NumCandidatesToKeep;
INT4 uvar_clusterOnScanline;

CHAR *uvar_gridFile;
REAL8 uvar_dopplermax;
INT4 uvar_RngMedWindow;
REAL8 uvar_refTime;
REAL8 uvar_internalRefTime;
INT4 uvar_SSBprecision;

INT4 uvar_minStartTime;
INT4 uvar_maxEndTime;
CHAR *uvar_workingDir;
REAL8 uvar_timerCount;
INT4 uvar_upsampleSFTs;

REAL8 BaryRefTime;


/* Defining a multi-IFO complex time series. Keeping the Real and Imaginary parts as seperate vectors, since it is necessary to interpolate them seperately */

typedef struct
{
  UINT4 length;                    /* Number of IFOs */
  REAL8Sequence** Real;           /* Real part of the time series */
  REAL8Sequence** Imag;           /* Imaginary part of the time series */
  REAL8 f_het;
  REAL8 deltaT;
  LIGOTimeGPS epoch; 
}MultiCOMPLEX8TimeSeries;

/* The Downsampled Time Series*/
MultiCOMPLEX8TimeSeries* TSeries;
REAL8 Fmin,Fmax;

/* A container for fftw_complex vector data */
typedef struct
{
  UINT4 length;
  fftw_complex* data;
}FFTWCOMPLEXSeries;

typedef struct
{
  UINT4 length;
  FFTWCOMPLEXSeries** data;
}MultiFFTWCOMPLEXSeries;

/* A buffer for resampling */
typedef struct
{
  const MultiDetectorStateSeries *multiDetStates;/**< buffer for each detStates (store pointer) and skypos */
  REAL8 Alpha, Delta;                           /**< skyposition of candidate */
  MultiSSBtimes *multiSSB;
  MultiSSBtimes *multiBinary;
  MultiAMCoeffs *multiAMcoef;
  MultiCmplxAMCoeffs *multiCmplxAMcoef;
  MultiFFTWCOMPLEXSeries *Saved_a,*Saved_b;      
}ReSampBuffer;


/* A contiguity structure required by the preprocessing function in order to store the information pertaining to the contiguity of SFTs and the gaps between them */

typedef struct
{
  UINT4  length;                    /* Number of Contiguous blocks */
  UINT4* NumContinuous;             /* Number of Contiguous SFTs in each block */
  REAL8* Gap;                       /* Gap between two Contiguous blocks in seconds */
}Contiguity;


/* ---------- local prototypes ---------- */
/* Resampling prototypes (start) */

LIGOTimeGPS REAL82GPS(REAL8 Time);
void ComputeFStat_resamp (LALStatus *,REAL8FrequencySeries *fstatVector, const PulsarDopplerParams *doppler, const MultiSFTVector *multiSFTs, const MultiNoiseWeights *multiWeights, const MultiDetectorStateSeries *multiDetStates,const ComputeFParams *params, ReSampBuffer *Buffer,ConfigVariables* GV);
void CalcTimeSeries(MultiSFTVector *multiSFTs);
MultiCOMPLEX8TimeSeries* XLALCreateMultiCOMPLEX8TimeSeries(UINT4 i);
void XLALDestroyMultiCOMPLEX8TimeSeries(MultiCOMPLEX8TimeSeries* T);
void XLALDestroyReSampBuffer ( ReSampBuffer *cfb);
void XLALDestroyMultiFFTWCOMPLEXSeries(MultiFFTWCOMPLEXSeries *X);
void XLALDestroyMultiCOMPLEX8TimeSeries(MultiCOMPLEX8TimeSeries *T);
void XLALDestroyFFTWCOMPLEXSeries(FFTWCOMPLEXSeries *X);
void XLALDestroyREAL8Sequence(REAL8Sequence *X);
void XLALDestroyMultiCmplxAMCoeffs(MultiCmplxAMCoeffs *X);
REAL8Sequence* XLALCreateREAL8Sequence(UINT4 length);
MultiFFTWCOMPLEXSeries *XLALCreateMultiFFTWCOMPLEXSeries(UINT4 length);
FFTWCOMPLEXSeries *XLALCreateFFTWCOMPLEXSeries(UINT4 length);
UINT4 factorial(UINT4 x);
REAL8 magsquare(fftw_complex f);
UINT4 Round(REAL8 X);
INT4 CombineSFTs(FFTWCOMPLEXSeries *L,SFTVector *sft_vect,REAL8 FMIN,REAL8 FMAX,INT4 number,INT4 startindex);
void ApplyWindow(REAL8Window *Win, FFTWCOMPLEXSeries *X);
void Reshuffle(FFTWCOMPLEXSeries *X);
void PrintL(FFTWCOMPLEXSeries* L,REAL8 min,REAL8 step);
void ApplyHetCorrection(REAL8Sequence *BaryTimes, REAL8Sequence *DetectorTimes, REAL8 dt, REAL8Sequence *Real, REAL8Sequence *Imag, REAL8 TSFT);
void ApplySpinDowns(REAL8 *SpinDowns, REAL8 dt, FFTWCOMPLEXSeries *FaIn, FFTWCOMPLEXSeries *FbIn,REAL8 StartTime, REAL8 BaryStartTime);
void ApplyAandB(REAL8Sequence *FineBaryTimes,REAL8Sequence *BaryTimes,REAL8Sequence *a,REAL8Sequence *b,REAL8Sequence *Real,REAL8Sequence *Imag,FFTWCOMPLEXSeries *FaIn, FFTWCOMPLEXSeries *FbIn);
double sinc(double t);
void retband(REAL8 t0, REAL8 dt, REAL8* t,REAL8* x, REAL8* y,UINT4 n,UINT4 size, UINT4 terms);
REAL8 strob(REAL8* Xdata, REAL8* Ydata, REAL8 X, UINT4 N1);
REAL8Sequence* ResampleSeries(REAL8Sequence *X_Real,REAL8Sequence *X_Imag,REAL8Sequence *Y_Real,REAL8Sequence *Y_Imag,REAL8 dt,REAL8Vector *BaryTimes, REAL8Sequence *DetectorTimes,REAL8 TSFT);


/* Resampling prototypes (end) */

int main(int argc,char *argv[]);
void initUserVars (LALStatus *);
void InitFStat ( LALStatus *, ConfigVariables *cfg );
void Freemem(LALStatus *,  ConfigVariables *cfg);

void WriteFStatLog (LALStatus *, CHAR *argv[], const CHAR *log_fname);
void checkUserInputConsistency (LALStatus *);
int outputBeamTS( const CHAR *fname, const AMCoeffs *amcoe, const DetectorStateSeries *detStates );
void InitEphemeris (LALStatus *, EphemerisData *edat, const CHAR *ephemDir, const CHAR *ephemYear, LIGOTimeGPS epoch, BOOLEAN isLISA);
void getUnitWeights ( LALStatus *, MultiNoiseWeights **multiWeights, const MultiSFTVector *multiSFTs );

int write_FstatCandidate_to_fp ( FILE *fp, const FstatCandidate *thisFCand );
int write_PulsarCandidate_to_fp ( FILE *fp,  const PulsarCandidate *pulsarParams, const FstatCandidate *Fcand );

int compareFstatCandidates ( const void *candA, const void *candB );
void getLogString ( LALStatus *status, CHAR **logstr, const ConfigVariables *cfg );

const char *va(const char *format, ...);        /* little var-arg string helper function */

/* ---------- scanline window functions ---------- */
scanlineWindow_t *XLALCreateScanlineWindow ( UINT4 windowWings );
void XLALDestroyScanlineWindow ( scanlineWindow_t *scanlineWindow );
int XLALAdvanceScanlineWindow ( const FstatCandidate *nextCand, scanlineWindow_t *scanWindow );
BOOLEAN XLALCenterIsLocalMax ( const scanlineWindow_t *scanWindow );


/*---------- empty initializers ---------- */
static const ConfigVariables empty_ConfigVariables;
static const FstatCandidate empty_FstatCandidate;
static const ReSampBuffer empty_ReSampBuffer;

/*----------------------------------------------------------------------*/
/* Function definitions start here */
/*----------------------------------------------------------------------*/

/**
 * MAIN function of ComputeFStatistic code.
 * Calculate the F-statistic over a given portion of the parameter-space
 * and write a list of 'candidates' into a file(default: 'Fstats').
 */
int main(int argc,char *argv[])
{
  LALStatus status = blank_status;      /* initialize status */
  /*MultiCOMPLEX8TimeSeries* TSeries;*/

  FILE *fpFstat = NULL;
  ReSampBuffer Buffer = empty_ReSampBuffer;
  REAL8 numTemplates, templateCounter;
  REAL8 tickCounter;
  time_t clock0;
  REAL8FrequencySeries *fstatVector = NULL;
  PulsarDopplerParams dopplerpos = empty_PulsarDopplerParams;           /* current search-parameters */
  FstatCandidate loudestFCand = empty_FstatCandidate, thisFCand = empty_FstatCandidate;
  UINT4 k;
  ConfigVariables GV = empty_ConfigVariables;           /**< global container for various derived configuration settings */

  lalDebugLevel = 0;
  vrbflg = 1;   /* verbose error-messages */

  /* set LAL error-handler */
  lal_errhandler = LAL_ERR_EXIT;

  /* register all user-variable */
  LAL_CALL (LALGetDebugLevel(&status, argc, argv, 'v'), &status);
  LAL_CALL (initUserVars(&status), &status);

  /* do ALL cmdline and cfgfile handling */
  LAL_CALL (LALUserVarReadAllInput(&status, argc,argv), &status);

  if (uvar_help)        /* if help was requested, we're done here */
    exit (0);

  /* set log-level */
  LogSetLevel ( lalDebugLevel );

  /* keep a log-file recording all relevant parameters of this search-run */
  if ( uvar_outputLogfile ) {
    LAL_CALL (WriteFStatLog ( &status, argv, uvar_outputLogfile ), &status );
  }

  /* do some sanity checks on the user-input before we proceed */
  LAL_CALL ( checkUserInputConsistency(&status), &status);

  /* Initialization the common variables of the code, */
  /* like ephemeries data and template grids: */
  LAL_CALL ( InitFStat(&status, &GV), &status);

  /* if a complete output of the F-statistic file was requested,
   * we open and prepare the output-file here */
  if (uvar_outputFstat)
    {
      if ( (fpFstat = fopen (uvar_outputFstat, "wb")) == NULL)
        {
          LALPrintError ("\nError opening file '%s' for writing..\n\n", uvar_outputFstat);
          return (COMPUTEFSTATISTIC_ESYS);
        }

      fprintf (fpFstat, "%s", GV.logstring );
    } /* if outputFstat */

  /* count number of templates */
  numTemplates = XLALNumDopplerTemplates ( GV.scanState );

  /* prepare Fstat-vector over frequencies to hold output results */
  {
    REAL8 dFreq = 1.0 / GV.multiDetStates->Tspan;
    /*UINT4 numFreqBins = floor( GV.FFTFreqBand / dFreq +1e-6) + 1;
      REAL8 Freq0 = GV.searchRegion.fkdot[0];*/
    UINT4 numFreqBins = floor(uvar_FreqBand / dFreq + 1e-6);
    REAL8 Freq0 = uvar_Freq;
    fstatVector = XLALCreateREAL8FrequencySeries ("Fstat vector", &GV.searchRegion.refTime, Freq0, dFreq, &empty_Unit, numFreqBins );
    if ( fstatVector == NULL ) {
      fprintf ( stderr, "Oops, out of memory!\n");
      return  COMPUTEFSTATISTIC_EMEM;
    }
  } /* setup Fstat-vector */

  /*----------------------------------------------------------------------
   * main loop: demodulate data for each point in the sky-position grid
   * and for each value of the frequency-spindown
   */
  templateCounter = 0.0;
  tickCounter = 0;
  clock0 = time(NULL);

  /*Call the CalcTimeSeries Function Here*/
  LogPrintf (LOG_DEBUG, "Calculating Time Series.\n");
  CalcTimeSeries(GV.multiSFTs);
  fprintf(stderr, "\n WARNING!!! Only the middle half of the band you asked for or is usable. Rest of it is destroyed by Interpolation. Please ask for a larger band. In the future, this will be done automatically. \n"); 
  LogPrintf (LOG_DEBUG, "Done Calculating Time Series.\n");

  LogPrintf (LOG_DEBUG, "Starting Main Resampling Loop.\n");
  while ( XLALNextDopplerPos( &dopplerpos, GV.scanState ) == 0 )
    {
      /* main function call: compute F-statistic over frequency-band  */ 
      LAL_CALL( ComputeFStat_resamp ( &status, fstatVector, &dopplerpos, GV.multiSFTs, GV.multiNoiseWeights,GV.multiDetStates, &GV.CFparams, &Buffer, &GV), &status );

      /* Progress meter */
      templateCounter += 1.0;
      if ( lalDebugLevel && ( ++tickCounter > uvar_timerCount) )
        {
          REAL8 diffSec = time(NULL) - clock0 ;  /* seconds since start of loop*/
          REAL8 taup = diffSec / templateCounter ;
          REAL8 timeLeft = (numTemplates - templateCounter) *  taup;
          tickCounter = 0.0;
          LogPrintf (LOG_DEBUG, "Progress: %g/%g = %.2f %% done, Estimated time left: %.0f s\n",
                     templateCounter, numTemplates, templateCounter/numTemplates * 100.0, timeLeft);
        }

      for ( k=0; k < fstatVector->data->length; k++)
        {
      
          REAL8 thisF = fstatVector->data->data[k];
          /*REAL8 thisFreq = fstatVector->f0 + k * fstatVector->deltaF;*/
          REAL8 thisFreq = uvar_Freq + k* fstatVector->deltaF;
          /*fprintf(stderr,"f0 = %f , deltaF = %f\n",fstatVector->f0,fstatVector->deltaF);*/
          /* sanity check on the result */
          if ( !finite ( thisF ) )
            {
              LogPrintf(LOG_CRITICAL, "non-finite F = %.16g\n", thisF );
              LogPrintf (LOG_CRITICAL, "[Alpha,Delta] = [%.16g,%.16g],\nfkdot=[%.16g,%.16g,%.16g,%16.g]\n",
                         dopplerpos.Alpha, dopplerpos.Delta,
                         thisFreq, dopplerpos.fkdot[1], dopplerpos.fkdot[2], dopplerpos.fkdot[3] );
          return -1;
            }

          /* propagate fkdot from internalRefTime back to refTime for outputting results */
          /* FIXE: only do this for candidates we're going to write out */
          dopplerpos.fkdot[0] = thisFreq;
          /*LAL_CALL ( LALExtrapolatePulsarSpins ( &status, dopplerpos.fkdot, GV.refTime, dopplerpos.fkdot, GV.searchRegion.refTime ), &status );*/
          dopplerpos.refTime = GV.refTime;

          /* correct normalization in --SignalOnly case:
           * we didn't normalize data by 1/sqrt(Tsft * 0.5 * Sh) in terms of
           * the single-sided PSD Sh: the SignalOnly case is characterized by
           * setting Sh->1, so we need to divide Fa,Fb by sqrt(0.5*Tsft) and F by (0.5*Tsft)
           */
          if ( uvar_SignalOnly )
            {
              REAL8 norm = 1.0 / sqrt( 0.5 * GV.Tsft );
              thisF *= norm * norm;
              thisF += 2;               /* compute E[2F]:= 4 + SNR^2 */
            }
          thisFCand.Fstat.F = thisF;
          thisFCand.doppler = dopplerpos;

          /* push new value onto scan-line buffer */
          XLALAdvanceScanlineWindow ( &thisFCand, GV.scanlineWindow );

          /* two types of threshold: fixed (TwoFThreshold) and dynamic (NumCandidatesToKeep) */
          if ( XLALCenterIsLocalMax ( GV.scanlineWindow )                                       /* must be 1D local maximum */
               && (2.0 * GV.scanlineWindow->center->Fstat.F >= uvar_TwoFthreshold) )    /* fixed threshold */
            {
              FstatCandidate *writeCand = GV.scanlineWindow->center;

              /* insert this into toplist if requested */
              if ( GV.FstatToplist  )                   /* dynamic threshold */
                {
                  if ( insert_into_toplist(GV.FstatToplist, (void*)writeCand ) )
                    LogPrintf ( LOG_DETAIL, "Added new candidate into toplist: 2F = %f\n", 2.0 * writeCand->Fstat.F );
                  else
                    LogPrintf ( LOG_DETAIL, "NOT added the candidate into toplist: 2F = %f\n", 2 * writeCand->Fstat.F );
                }
              else if ( fpFstat )                               /* no toplist :write out immediately */
                {
                  if ( write_FstatCandidate_to_fp ( fpFstat, writeCand ) != 0 )
                    {
                      LogPrintf (LOG_CRITICAL, "Failed to write candidate to file.\n");
                      return -1;
                    }
                } /* if outputFstat */

            } /* if 2F > threshold */

          /* separately keep track of loudest candidate (for --outputLoudest) */
          if ( thisFCand.Fstat.F > loudestFCand.Fstat.F )
            loudestFCand = thisFCand;

        } /* inner loop about frequency-bins from resampling frequ-band */

    } /* while more Doppler positions to scan */

  /* ----- if using toplist: sort and write it out to file now ----- */
  if ( fpFstat && GV.FstatToplist )
    {
      UINT4 el;

      /* sort toplist */
      LogPrintf ( LOG_DEBUG, "Sorting toplist ... ");
      qsort_toplist ( GV.FstatToplist, compareFstatCandidates );
      LogPrintfVerbatim ( LOG_DEBUG, "done.\n");

      for ( el=0; el < GV.FstatToplist->elems; el ++ )
        {
          const FstatCandidate *candi;
          if ( ( candi = (const FstatCandidate *) toplist_elem ( GV.FstatToplist, el )) == NULL ) {
            LogPrintf ( LOG_CRITICAL, "Internal consistency problems with toplist: contains fewer elements than expected!\n");
            return -1;
          }
          if ( write_FstatCandidate_to_fp ( fpFstat, candi ) != 0 )
            {
              LogPrintf (LOG_CRITICAL, "Failed to write candidate to file.\n");
              return -1;
            }
        } /* for el < elems in toplist */

    } /* if fpFstat && toplist */

  XLALDestroyMultiCOMPLEX8TimeSeries(TSeries);

  if ( fpFstat )
    {
      fprintf (fpFstat, "%%DONE\n");
      fclose (fpFstat);
      fpFstat = NULL;
    }

  /* ----- estimate amplitude-parameters for the loudest canidate and output into separate file ----- */
  if ( uvar_outputLoudest )
    {
      FILE *fpLoudest;
      PulsarCandidate pulsarParams = empty_PulsarCandidate;
      pulsarParams.Doppler = loudestFCand.doppler;

      LAL_CALL(LALEstimatePulsarAmplitudeParams (&status, &pulsarParams, &loudestFCand.Fstat, &GV.searchRegion.refTime, &loudestFCand.Mmunu ),
               &status );

      if ( (fpLoudest = fopen (uvar_outputLoudest, "wb")) == NULL)
        {
          LALPrintError ("\nError opening file '%s' for writing..\n\n", uvar_outputLoudest);
          return COMPUTEFSTATISTIC_ESYS;
        }

      /* write header with run-info */
      fprintf (fpLoudest, "%s", GV.logstring );

      /* write this 'candidate' to disc */
      if ( write_PulsarCandidate_to_fp ( fpLoudest,  &pulsarParams, &loudestFCand) != XLAL_SUCCESS )
        {
          LogPrintf(LOG_CRITICAL, "call to write_PulsarCandidate_to_fp() failed!\n");
          return COMPUTEFSTATISTIC_ESYS;
        }

      fclose (fpLoudest);

      gsl_matrix_free ( pulsarParams.AmpFisherMatrix );

    } /* write loudest candidate to file */

  LogPrintf (LOG_DEBUG, "Search finished.\n");

  /* Free memory */
  LogPrintf (