Commit f72097f6 authored by Miroslav Shaltev's avatar Miroslav Shaltev
Browse files

fix FStatNomad cache writer, AfterMath read data and compute metric

parent 85986dd6
......@@ -55,7 +55,22 @@ int main(int argc, char *argv[]) {
MAFTERMath = new AFTERMath(argc,argv);
MAFTERMath->XLALInitUserVars(argc,argv);
// read input data
MAFTERMath->ReadInputData();
// compute metric
MAFTERMath->SegList(LocalXLALReadSegmentsFromFile(MAFTERMath->uvar->segmentList));
MAFTERMath->metric = gsl_matrix_calloc(MAFTERMath->DimMetric(MAFTERMath->uvar->sType,USETYPE_SEARCH),MAFTERMath->DimMetric(MAFTERMath->uvar->sType,USETYPE_SEARCH));
MAFTERMath->ComputeMetric(MAFTERMath->Data(0));
// compute distance for each point to each point
// compare 2F for each point to each point
return EXIT_SUCCESS;
}
......@@ -76,17 +91,27 @@ INT4 AFTERMath::XLALInitUserVars ( int argc, char *argv[] )
uvar->fnameout = (char*)LALCalloc( strlen(FNAMEOUT) + 1, sizeof(CHAR) );
strcpy(uvar->fnameout, FNAMEOUT);
uvar->sType = STYPE_NONE;
uvar->refTime = 0;
uvar->metricType = MTAVFSTAT;
if ( (uvar->MetricIFO = XLALCreateStringVector ( "H1,L1", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
}
/* register user input variables */
XLALRegisterUvarMember( version,BOOLEAN, 'V', SPECIAL, "Output version information");
XLALRegisterUvarMember( fnamein, STRING, 0, OPTIONAL, "Input filename");
XLALRegisterUvarMember( fnameout, STRING, 0, OPTIONAL, "Output filename");
XLALRegisterUvarMember( ephemEarth, STRING, 0, OPTIONAL, "Location of Earth ephemeris file");
XLALRegisterUvarMember( ephemSun, STRING, 0, OPTIONAL, "Location of Sun ephemeris file");
XLALRegisterUvarMember( fnameout, STRING, 0, OPTIONAL, "Output filename");
XLALRegisterUvarMember( ephemEarth, STRING, 0, OPTIONAL, "Location of Earth ephemeris file");
XLALRegisterUvarMember( ephemSun, STRING, 0, OPTIONAL, "Location of Sun ephemeris file");
XLALRegisterUvarMember( segmentList, STRING, 0, OPTIONAL, "File containing a segment list used to compute Fisher ellipsoid: lines of form <startGPS endGPS duration[h] NumSFTs>");
XLALRegisterUvarMember( sType, INT4, 0, OPTIONAL, "Search / metric / cache type: Isolated ( 0 - ALPHA_DELTA, 1 - ALPHA_DELTA_FREQ_F1DOT, 2 - ALPHA_DELTA_FREQ_F1DOT_F2DOT, 3 - ALPHA_DELTA_FREQ_F1DOT_F2DOT_F3DOT, 4 - FREQ_F1DOT, 5 - FREQ_F1DOT_F2DOT, 6 - FREQ_F1DOT_F2DOT_F3DOT ), Binary ( 7 - ALPHA_DELTA_FREQ_ASINI_TASC_OMEGA_KAPPA_ETA, 8 - ALPHA_DELTA_FREQ_F1DOT_ASINI_TASC_OMEGA_KAPPA_ETA, 9 - ALPHA_DELTA_FREQ_F1DOT_F2DOT_ASINI_TASC_OMEGA_KAPPA_ETA, 10 - ALPHA_DELTA_FREQ_F1DOT_F2DOT_F3DOT_ASINI_TASC_OMEGA_KAPPA_ETA, 11 - FREQ_ASINI_TASC_OMEGA_KAPPA_ETA, 12 - FREQ_F1DOT_ASINI_TASC_OMEGA_KAPPA_ETA, 13 - FREQ_F1DOT_F2DOT_ASINI_TASC_OMEGA_KAPPA_ETA, 14 - FREQ_F1DOT_F2DOT_F3DOT_ASINI_TASC_OMEGA_KAPPA_ETA )");
XLALRegisterUvarMember( refTime, REAL8, 0, OPTIONAL, "Ref. time for pulsar pars [Default: 0]");
XLALRegisterUvarMember( metricType, INT4, 0, OPTIONAL, "Type of metric to compute: 0=phase-metric, 1=average F-metric");
XLALRegisterUvarMember( MetricIFO, STRINGVector, 0, OPTIONAL, "Use this IFOs to compute metric");
/* read cmdline & cfgfile */
BOOLEAN should_exit = 0;
XLAL_CHECK( XLALUserVarReadAllInput( &should_exit, argc, argv ) == XLAL_SUCCESS, XLAL_EFUNC );
......@@ -134,10 +159,9 @@ AFTERMath::AFTERMath(int argc, char *argv[]) {
UserVariables = empty_UserInput; // initializes this struct to {0}
/* user variables */
edat = NULL;
}
AFTERMath::~AFTERMath() {}
......
......@@ -36,6 +36,9 @@
/* lalapps includes */
#include <lalapps.h>
#include <vector>
#include <gsl/gsl_vector.h>
#define AFTERMATH_ENORM 0
#define AFTERMATH_ESUB 1
......@@ -87,14 +90,30 @@ typedef enum {
STYPE_ISO_FREQ_F1DOT_F2DOT_F3DOT,
STYPE_BNS_ALPHA_DELTA_FREQ_ASINI_TASC_OMEGA_KAPPA_ETA,
STYPE_BNS_ALPHA_DELTA_FREQ_F1DOT_ASINI_TASC_OMEGA_KAPPA_ETA,
STYPE_BNS_ALPHA_DELTA_FREQ_F1DOT_F2DOT_ASINI_TASC_OMEGA_KAPPA_ETA,
STYPE_BNS_ALPHA_DELTA_FREQ_F1DOT_F2DOT_F3DOT_ASINI_TASC_OMEGA_KAPPA_ETA,
STYPE_BNS_ALPHA_DELTA_FREQ_F1DOT_F2DOT_ASINI_TASC_OMEGA_KAPPA_ETA,
STYPE_BNS_ALPHA_DELTA_FREQ_F1DOT_F2DOT_F3DOT_ASINI_TASC_OMEGA_KAPPA_ETA,
STYPE_BNS_FREQ_ASINI_TASC_OMEGA_KAPPA_ETA,
STYPE_BNS_FREQ_F1DOT_ASINI_TASC_OMEGA_KAPPA_ETA,
STYPE_BNS_FREQ_F1DOT_F2DOT_ASINI_TASC_OMEGA_KAPPA_ETA,
STYPE_BNS_FREQ_F1DOT_F2DOT_F3DOT_ASINI_TASC_OMEGA_KAPPA_ETA
STYPE_BNS_FREQ_F1DOT_F2DOT_ASINI_TASC_OMEGA_KAPPA_ETA,
STYPE_BNS_FREQ_F1DOT_F2DOT_F3DOT_ASINI_TASC_OMEGA_KAPPA_ETA
} STypeID;
typedef enum {
DP_NONE = -1,
DP_ALPHA,
DP_DELTA,
DP_FREQ,
DP_F1DOT,
DP_F2DOT,
DP_F3DOT,
DP_TPSSB,
DP_ARGP,
DP_ASINI,
DP_ECC,
DP_PERIOD,
DP_STAT
} DirectionPositionID;
/******************************************************
* Protection against C++ name mangling
......@@ -110,7 +129,11 @@ extern "C" {
CHAR *fnameout;
CHAR *ephemEarth;
CHAR *ephemSun;
CHAR *segmentList;
CHAR *segmentList;
INT4 sType;
REAL8 refTime;
INT4 metricType;
LALStringVector* MetricIFO;
} UserInput_t;
#ifdef __cplusplus
} /* Close C++ protection */
......@@ -129,10 +152,6 @@ typedef struct {
REAL8 ecc;
REAL8 period;
REAL8 stat;
REAL8 h0;
REAL8 cosi;
REAL8 psi;
REAL8 phi0;
} ExPoint;
typedef struct
......@@ -149,6 +168,7 @@ void PrintPoint(INT4 nseg,const ExPoint *p,const CHAR* msg, BOOLEAN bns);
LALSegList *LocalXLALReadSegmentsFromFile ( const char *fname );
class AFTERMath {
public:
AFTERMath(int argc, char *argv[]);
......@@ -157,28 +177,113 @@ public:
UserInput_t UserVariables;
UserInput_t *uvar;
void BinarySearch(BOOLEAN v) {
var_BinarySearch = v;
};
BOOLEAN BinarySearch() {
return var_BinarySearch;
};
void MetricIFO(LALStringVector* v) {
var_MetricIFO = v;
INT4 ReadInputData();
gsl_matrix *metric;
void SegList(LALSegList *v) { segList = v; };
INT4 DimMetric(INT4 SType,INT4 UseType) {
switch (SType) {
case STYPE_ISO_ALPHA_DELTA:
if (UseType == USETYPE_EXTENDS) {
return 2;
}
else {
return 3;
}
break;
case STYPE_ISO_ALPHA_DELTA_FREQ_F1DOT:
if (UseType == USETYPE_EXTENDS) {
return 4;
}
else {
return 5;
}
break;
case STYPE_ISO_ALPHA_DELTA_FREQ_F1DOT_F2DOT:
if (UseType == USETYPE_EXTENDS) {
return 5;
}
else {
return 6;
}
break;
case STYPE_ISO_ALPHA_DELTA_FREQ_F1DOT_F2DOT_F3DOT:
if (UseType == USETYPE_EXTENDS) {
return 6;
}
else {
return 7;
}
break;
case STYPE_ISO_FREQ_F1DOT:
return 2;
break;
case STYPE_ISO_FREQ_F1DOT_F2DOT:
return 3;
break;
case STYPE_ISO_FREQ_F1DOT_F2DOT_F3DOT:
return 4;
break;
case STYPE_BNS_ALPHA_DELTA_FREQ_ASINI_TASC_OMEGA_KAPPA_ETA:
if (UseType == USETYPE_EXTENDS) {
return 8;
}
else {
return 9;
}
break;
case STYPE_BNS_ALPHA_DELTA_FREQ_F1DOT_ASINI_TASC_OMEGA_KAPPA_ETA:
if (UseType == USETYPE_EXTENDS) {
return 9;
}
else {
return 10;
}
break;
case STYPE_BNS_FREQ_ASINI_TASC_OMEGA_KAPPA_ETA:
return 6;
break;
case STYPE_BNS_FREQ_F1DOT_ASINI_TASC_OMEGA_KAPPA_ETA:
return 7;
break;
default:
return -1;
}
};
INT4 ComputeMetric(gsl_vector *v);
gsl_vector* Data(UINT8 v){
return indata.at(v);
};
private:
CHAR *VCSInfoString;
BOOLEAN var_BinarySearch;
LALStringVector* var_MetricIFO;
EphemerisData *edat;
LALSegList *segList;
INT4 compute_phase_metric(DopplerPhaseMetric *phase_metric, gsl_matrix *g, LALSegList *seglist, LIGOTimeGPS refTimeGPS, INT4 MType, ExPoint *p, INT4 EType);
INT4 compute_fstat_metric(DopplerFstatMetric *fstat_metric, gsl_matrix *g, LALSegList *seglist, LIGOTimeGPS refTimeGPS, INT4 SType, ExPoint *p, INT4 EType);
DopplerPhaseMetric *phase_metric;
DopplerFstatMetric *fstat_metric;
BOOLEAN check_matrix(gsl_matrix *g);
void DoDisplayMetric(gsl_matrix *m);
std::vector<gsl_vector*> indata;
INT4 compute_metric(gsl_matrix *g, LALSegList *seglist, INT4 SType, ExPoint *p, INT4 EType, INT4 MAType);
bool BinarySearch(){
if (uvar->sType > STYPE_BNS_ALPHA_DELTA_FREQ_ASINI_TASC_OMEGA_KAPPA_ETA) {
return true;
}
else {
return false;
}
}
};
#endif /* Double-include protection. */
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include "AfterMath.h"
/** Function to read a segment list from given filename, returns a *sorted* SegmentList
......@@ -82,3 +87,34 @@ LocalXLALReadSegmentsFromFile ( const char *fname /**< name of file containing s
return segList;
} /* XLALReadSegmentsFromFile() */
INT4 AFTERMath::ReadInputData(){
std::string line;
std::ifstream dfile (uvar->fnamein);
if (dfile.is_open())
{
while ( getline (dfile,line) )
{
gsl_vector *c;
std::vector<double> v;
std::stringstream in(line);
double temp;
while(in >> temp) {
v.push_back(temp);
}
c = gsl_vector_calloc(v.size());
gsl_vector_set_zero(c);
for (int i=0; i < v.size(); i++){
gsl_vector_set(c,i,v.at(i));
}
indata.push_back(c);
// std::cout << line << '\n';
}
dfile.close();
printf("Found %d records in %s.\n",indata.size(),uvar->fnamein);
}
else std::cout << "Unable to open file";
return 0;
}
......@@ -34,11 +34,21 @@ BOOLEAN AFTERMath::check_matrix(gsl_matrix *g) {
return TRUE;
}
INT4 AFTERMath::compute_metric(gsl_matrix *g, LALSegList *seglist, INT4 SType, ExPoint *p, INT4 EType, INT4 MAType) {
INT4 AFTERMath::compute_phase_metric(DopplerPhaseMetric *phase_metric, gsl_matrix *g, LALSegList *seglist, LIGOTimeGPS refTimeGPS, INT4 SType, ExPoint *p, INT4 EType) {
LIGOTimeGPS refTimeGPS;
XLAL_INIT_DECL(refTimeGPS);
XLALGPSSetREAL8(&refTimeGPS, uvar->refTime);
if ( MAType == MTPHASE ) {
DopplerPhaseMetric *phase_metric;
}
else if ( MAType == MTAVFSTAT ) {
DopplerFstatMetric *fstat_metric;
}
if (seglist != NULL) {
// XLALGPSSetREAL8( &refTimeGPS, usefulParams.refTime );
int Nseg = seglist->length;
// #ifdef DEBUG
PrintPoint(Nseg,p,"compute metric in this point",BinarySearch());
......@@ -212,7 +222,7 @@ INT4 AFTERMath::compute_phase_metric(DopplerPhaseMetric *phase_metric, gsl_matri
dop->refTime = refTimeGPS;
if ( XLALParseMultiLALDetector ( &config.multiIFO, var_MetricIFO ) != XLAL_SUCCESS ) {
if ( XLALParseMultiLALDetector ( &config.multiIFO, uvar->MetricIFO ) != XLAL_SUCCESS ) {
LogPrintf (LOG_CRITICAL, "XLALParseMultiDetectorInfo() failed to parse detector names and/or weights. errno = %d.\n\n", xlalErrno);
return xlalErrno;
}
......@@ -229,6 +239,7 @@ INT4 AFTERMath::compute_phase_metric(DopplerPhaseMetric *phase_metric, gsl_matri
/* ----- compute metric full metric + Fisher matrix ---------- */
if ( MAType == MTPHASE ) {
if ( (phase_metric = XLALComputeDopplerPhaseMetric ( &metricParams, edat )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Something failed in XLALDopplerPhaseMetric(). xlalErrno = %d\n\n", xlalErrno);
return xlalErrno;
......@@ -248,217 +259,8 @@ INT4 AFTERMath::compute_phase_metric(DopplerPhaseMetric *phase_metric, gsl_matri
exit(-1);
}
}
if (check_matrix(g)) {
#ifdef DEBUG
DoDisplayMetric(g);
#endif
}
else {
LogPrintf(LOG_CRITICAL,"check_matrix failed for the rescaled metric. Abort!\n");
exit(-1);
}
}
else {
LogPrintf(LOG_NORMAL,"Cannot compute metric from NULL segment list!\n");
return -1;
}
}
INT4 AFTERMath::compute_fstat_metric(DopplerFstatMetric *fstat_metric, gsl_matrix *g, LALSegList *seglist, LIGOTimeGPS refTimeGPS, INT4 SType, ExPoint *p, INT4 EType) {
if (seglist != NULL) {
// XLALGPSSetREAL8( &refTimeGPS, usefulParams.refTime );
int Nseg = seglist->length;
// #ifdef DEBUG
PrintPoint(Nseg,p,"compute metric in this point",BinarySearch());
// #endif
LogPrintf(LOG_NORMAL,"for refTime: %d\n",p->f1dot,refTimeGPS.gpsSeconds);
int s;
CMConfigVariables config = empty_ConfigVariables;
DopplerMetricParams XLAL_INIT_DECL(metricParams);
LALStringVector* coords;
switch (SType) {
case STYPE_ISO_ALPHA_DELTA:
if (EType == USETYPE_EXTENDS) {
if ( (coords = XLALCreateStringVector ( "alpha", "delta", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
}
else {
if ( (coords = XLALCreateStringVector ( "n3x_equ", "n3y_equ", "n3z_equ", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
}
break;
case STYPE_ISO_ALPHA_DELTA_FREQ_F1DOT:
if (EType == USETYPE_EXTENDS) {
if ( (coords = XLALCreateStringVector ( "alpha", "delta", "freq", "f1dot", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
}
else {
if ( (coords = XLALCreateStringVector ( "n3x_equ", "n3y_equ", "n3z_equ", "freq", "f1dot", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
}
break;
case STYPE_ISO_ALPHA_DELTA_FREQ_F1DOT_F2DOT:
if (EType == USETYPE_EXTENDS) {
if ( (coords = XLALCreateStringVector ( "alpha", "delta", "freq", "f1dot", "f2dot", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
}
else {
if ( (coords = XLALCreateStringVector ( "n3x_equ", "n3y_equ", "n3z_equ", "freq", "f1dot", "f2dot", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
}
break;
case STYPE_ISO_ALPHA_DELTA_FREQ_F1DOT_F2DOT_F3DOT:
if (EType == USETYPE_EXTENDS) {
if ( (coords = XLALCreateStringVector ( "alpha", "delta", "freq", "f1dot", "f2dot", "f3dot", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
}
else {
if ( (coords = XLALCreateStringVector ( "n3x_equ", "n3y_equ", "n3z_equ", "freq", "f1dot", "f2dot", "f3dot", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
}
break;
case STYPE_ISO_FREQ_F1DOT:
if ( (coords = XLALCreateStringVector ( "freq", "f1dot", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
break;
case STYPE_ISO_FREQ_F1DOT_F2DOT:
if ( (coords = XLALCreateStringVector ( "freq", "f1dot", "f2dot", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
break;
case STYPE_ISO_FREQ_F1DOT_F2DOT_F3DOT:
if ( (coords = XLALCreateStringVector ( "freq", "f1dot", "f2dot", "f3dot", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
break;
case STYPE_BNS_ALPHA_DELTA_FREQ_ASINI_TASC_OMEGA_KAPPA_ETA:
if ( (coords = XLALCreateStringVector ( "alpha", "delta", "freq", "asini", "tasc", "porb", "kappa", "eta", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
break;
case STYPE_BNS_ALPHA_DELTA_FREQ_F1DOT_ASINI_TASC_OMEGA_KAPPA_ETA:
if ( (coords = XLALCreateStringVector ( "alpha", "delta", "freq", "f1dot", "asini", "tasc", "porb", "kappa", "eta", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
break;
case STYPE_BNS_ALPHA_DELTA_FREQ_F1DOT_F2DOT_ASINI_TASC_OMEGA_KAPPA_ETA:
if ( (coords = XLALCreateStringVector ( "alpha", "delta", "freq", "f1dot", "f2dot", "asini", "tasc", "porb", "kappa", "eta", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
break;
case STYPE_BNS_ALPHA_DELTA_FREQ_F1DOT_F2DOT_F3DOT_ASINI_TASC_OMEGA_KAPPA_ETA:
if ( (coords = XLALCreateStringVector ( "alpha", "delta", "freq", "f1dot", "f2dot", "f3dot", "asini", "tasc", "porb", "kappa", "eta", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
break;
case STYPE_BNS_FREQ_ASINI_TASC_OMEGA_KAPPA_ETA:
if ( (coords = XLALCreateStringVector ( "freq", "asini", "tasc", "porb", "kappa", "eta", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
break;
case STYPE_BNS_FREQ_F1DOT_ASINI_TASC_OMEGA_KAPPA_ETA:
if ( (coords = XLALCreateStringVector ( "freq", "f1dot", "asini", "tasc", "porb", "kappa", "eta", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
break;
case STYPE_BNS_FREQ_F1DOT_F2DOT_ASINI_TASC_OMEGA_KAPPA_ETA:
if ( (coords = XLALCreateStringVector ( "freq", "f1dot", "f2dot", "asini", "tasc", "porb", "kappa", "eta", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
break;
case STYPE_BNS_FREQ_F1DOT_F2DOT_F3DOT_ASINI_TASC_OMEGA_KAPPA_ETA:
if ( (coords = XLALCreateStringVector ( "freq", "f1dot", "f2dot", "f3dot", "asini", "tasc", "porb", "kappa", "eta", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
break;
default:
LogPrintf (LOG_CRITICAL, "Uknown SType: %d\n", SType );
return -1;
}
if ( XLALDopplerCoordinateNames2System ( &config.coordSys, coords ) ) {
LogPrintf (LOG_CRITICAL, "Call to XLALDopplerCoordinateNames2System() failed. errno = %d\n\n", xlalErrno );
return xlalErrno;
}
config.signalParams.Amp.h0 = 1.;
config.signalParams.Amp.cosi = 0.;
config.signalParams.Amp.psi = 0.;
config.signalParams.Amp.phi0 = 0.;
PulsarDopplerParams *dop = &(config.signalParams.Doppler);
XLAL_INIT_DECL(*dop);
dop->Alpha = p->alpha;
dop->Delta = p->delta;
dop->fkdot[0] = p->freq;
dop->fkdot[1] = p->f1dot;
dop->fkdot[2] = p->f2dot;
dop->fkdot[3] = p->f3dot;
if (BinarySearch()) {
dop->tp.gpsSeconds = p->TpSSB.gpsSeconds;
dop->tp.gpsNanoSeconds = p->TpSSB.gpsNanoSeconds;
dop->argp = p->argp;
dop->asini = p->asini;
dop->ecc = p->ecc;
dop->period = p->period;
}
dop->refTime = refTimeGPS;
if ( XLALParseMultiLALDetector ( &config.multiIFO, var_MetricIFO ) != XLAL_SUCCESS ) {
LogPrintf (LOG_CRITICAL, "XLALParseMultiDetectorInfo() failed to parse detector names and/or weights. errno = %d.\n\n", xlalErrno);
return xlalErrno;
}
metricParams.coordSys = config.coordSys;
metricParams.detMotionType = DETMOTION_ORBIT; /* FIXME spin + motion */
// FIXME metricParams.metricType = (MetricType_t)MAType;
metricParams.multiIFO = config.multiIFO;
metricParams.multiNoiseFloor = config.multiNoiseFloor;
metricParams.signalParams = config.signalParams;
metricParams.projectCoord = - 1; /* no projection */
metricParams.segmentList = *seglist;
/* ----- compute metric full metric + Fisher matrix ---------- */
}
else if (MAType == MTAVFSTAT) {
if ( (fstat_metric = XLALComputeDopplerFstatMetric ( &metricParams, edat )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Something failed in XLALDopplerFstatMetric(). xlalErrno = %d\n\n", xlalErrno);
return xlalErrno;
......@@ -480,12 +282,13 @@ INT4 AFTERMath::compute_fstat_metric(DopplerFstatMetric *fstat_metric, gsl_matri
}
}
if (check_matrix(g)) {
#ifdef DEBUG
// #ifdef DEBUG
DoDisplayMetric(g);
#endif
// #endif
}
else {
LogPrintf(LOG_CRITICAL,"check_matrix failed for the rescaled metric. Abort!\n");
......@@ -498,3 +301,21 @@ INT4 AFTERMath::compute_fstat_metric(DopplerFstatMetric *fstat_metric, gsl_matri
}
}
INT4 AFTERMath::ComputeMetric(gsl_vector *v){
ExPoint p;
p.alpha = gsl_vector_get(v,DP_ALPHA);
p.delta = gsl_vector_get(v,DP_DELTA);
p.freq = gsl_vector_get(v,DP_FREQ);
p.f1dot = gsl_vector_get(v,DP_F1DOT);