libFStatNomad.cpp

/*
 *  Copyright (C) 2012 Miroslav Shaltev.
 *
 *  May include code snippets by:
 *      Holger Pletsch, Karl Wette, Reinhard Prix,
 *      Badri Krishnan, Alicia Sintes, Bernd Machenschalk.
 *
 *  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
 *
 */


#define __STDC_CONSTANT_MACROS

#include "libFStatNomad.h"

static FstatCandidate empty_FstatCandidate;
static ConfigVariables empty_ConfigVariables;

using namespace std;

REAL8 FStatNomad::DoPredictFStat(INT4 fstatsearchtype, LALSegList *segList) {

    AntennaPatternMatrix Mmunu;

    init_SP(fstatsearchtype,segList);
    SkyPosition skypos;
    MultiAMCoeffs *multiAMcoef = NULL;
    UINT4 X, i, numSFTs;
    LALStatus lstatus = blank_status;
    /* normalize skyposition: correctly map into [0,2pi]x[-pi/2,pi/2] */
    skypos.longitude = SignalPoint.alpha;
    skypos.latitude = SignalPoint.delta;
    skypos.system = COORDINATESYSTEM_EQUATORIAL;
    LALNormalizeSkyPosition ( &lstatus, &skypos, &skypos);

    REAL8 rho2 = 0;

    for (i = 0; i < Fstat_in_vec->length; i++ ) {
        MultiSFTVector*  multiSFTs = XLALLoadMultiSFTs( catalogSeq.data + i, freqmin, freqmax );
        MultiNoiseWeights *multiNoiseWeights = NULL;
        MultiPSDVector *psd = NULL;
        LALNormalizeMultiSFTVect ( &lstatus, &psd, multiSFTs, blocksRngMed() );
        LALComputeMultiNoiseWeights ( &lstatus, &multiNoiseWeights, psd, blocksRngMed(), 0 );
        LALDestroyMultiPSDVector ( &lstatus, &psd );
        MultiDetectorStateSeries *multiDetStates = NULL;
        LALGetMultiDetectorStates( &lstatus, &multiDetStates, multiSFTs, edat);


        LALGetMultiAMCoeffs ( &lstatus, &multiAMcoef, multiDetStates, skypos );
        /* noise-weighting of Antenna-patterns and compute A,B,C */
        if ( XLALWeightMultiAMCoeffs ( multiAMcoef, multiNoiseWeights ) != XLAL_SUCCESS ) {
            LogPrintf (LOG_NORMAL, "XLALWeightMultiAMCoeffs() failed with error = %d\n\n", xlalErrno );
            ABORT ( &lstatus, FSTATFCNOMAD_EXLAL, FSTATFCNOMAD_MSGEXLAL );
        }

        /* OK: we only need the antenna-pattern matrix M_mu_nu */