Commit 53194e09 authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

adding glitch and raw data to output frame


git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@787 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent 5d9f74b3
......@@ -63,7 +63,7 @@ void parse_command_line(int argc, char **argv, struct Data *data);
static REAL8TimeSeries *readTseries(CHAR *cachefile, CHAR *channel, LIGOTimeGPS start, REAL8 length);
static void output_frame(REAL8TimeSeries *timeData, CHAR *frameType, CHAR *channel, CHAR *ifo);
static void output_frame(REAL8TimeSeries *timeData, REAL8TimeSeries *timeRes, REAL8TimeSeries *timeGlitch, CHAR *frameType, CHAR *ifo);
void dfour1(double data[], unsigned long nn, int isign);
void drealft(double data[], unsigned long n, int isign);
......@@ -117,7 +117,9 @@ int main(int argc, char *argv[])
// ------------------------------------------------- //
parse_command_line(argc, argv, data);
REAL8TimeSeries* timeRes=NULL;
REAL8TimeSeries* timeData=NULL;
REAL8TimeSeries* timeGlitch=NULL;
LIGOTimeGPS epoch;
// Set time & retrieve data by reading frame cache
......@@ -126,16 +128,26 @@ int main(int argc, char *argv[])
fprintf(stdout," Type: %s\n",data->fr_type);
fprintf(stdout," Channel: %s\n",data->fr_chanl);
fprintf(stdout,"\n");
timeData=readTseries(data->fr_cache,data->fr_chanl,epoch,data->fr_seglen);
timeRes = readTseries(data->fr_cache,data->fr_chanl,epoch,data->fr_seglen);
timeData = readTseries(data->fr_cache,data->fr_chanl,epoch,data->fr_seglen);
timeGlitch = readTseries(data->fr_cache,data->fr_chanl,epoch,data->fr_seglen);
char * version="T1700406_v3";
char outframeType[128];
char outframeChannel[128];
char outframeGlitchChannel[128];
sprintf(outframeType, "%s_%s", data->fr_type, version);
sprintf(outframeChannel,"%s_%s", data->fr_chanl, version);
sprintf(outframeGlitchChannel,"%s_glitch", data->fr_chanl);
/* set the channel name */
strncpy(timeData->name, data->fr_chanl, LALNameLength);
strncpy(timeRes->name, outframeChannel, LALNameLength);
strncpy(timeGlitch->name, outframeGlitchChannel, LALNameLength);
fprintf(stdout,"Creating glitch-subtracted frame:\n");
fprintf(stdout," Type: %s\n",outframeType);
fprintf(stdout," Channel: %s\n",outframeChannel);
......@@ -289,12 +301,15 @@ int main(int argc, char *argv[])
clean_data[i] = timeData->data->data[i] - norm*glitch_model_padded[i];
//fill LAL data structure with cleaned data for frame-making
timeData->data->data[i] = clean_data[i];
timeRes->data->data[i] = clean_data[i];
//fill another structure with glitch model
timeGlitch->data->data[i] = norm*glitch_model_padded[i];
}
// Output cleaned data!
output_frame(timeData, outframeType, outframeChannel, data->ifo);
output_frame(timeData, timeRes, timeGlitch, outframeType, data->ifo);
XLALDestroyREAL8TimeSeries(timeData);
......@@ -540,7 +555,11 @@ static REAL8TimeSeries *readTseries(CHAR *cachefile, CHAR *channel, LIGOTimeGPS
}
static void output_frame(REAL8TimeSeries *timeData, CHAR *frameType, CHAR *channel, CHAR *ifo)
static void output_frame(REAL8TimeSeries *timeData,
REAL8TimeSeries *timeRes,
REAL8TimeSeries *timeGlitch,
CHAR *frameType,
CHAR *ifo)
{
CHAR fname[2048];
INT4 duration;
......@@ -574,14 +593,12 @@ static void output_frame(REAL8TimeSeries *timeData, CHAR *frameType, CHAR *chann
duration = gpsEnd - gpsStart;
snprintf( fname, FILENAME_MAX, "%c-%s-%d-%d.gwf", ifo[0], frameType, gpsStart, duration );
/* set the channel name */
strncpy(timeData->name, channel, LALNameLength);
/* define frame */
frame = XLALFrameNew( &timeData->epoch, duration, "LIGO", 0, 1,
detectorFlags );
frame = XLALFrameNew( &timeData->epoch, duration, "LIGO", 0, 1, detectorFlags );
/* add channel to frame */
XLALFrameAddREAL8TimeSeriesSimData( frame, timeRes );
XLALFrameAddREAL8TimeSeriesSimData( frame, timeGlitch );
XLALFrameAddREAL8TimeSeriesSimData( frame, timeData );
fprintf( stdout, "Writing data to frame: '%s'\n", fname );
......
......@@ -800,7 +800,7 @@ void print_bayesline_spectrum(char filename[], double *f, double *power, double
int i;
FILE *waveout = fopen(filename,"w");
for(i=0; i<N; i++) fprintf(waveout,"%e %e %e %e\n", f[i], power[i], Sbase[i], Sline[i]);
for(i=0; i<N; i++) fprintf(waveout,"%.16g %.16g %.16g %.16g\n", f[i], power[i], Sbase[i], Sline[i]);
fclose(waveout);
}
......@@ -1150,7 +1150,7 @@ void print_time_domain_waveforms(char filename[], double *h, int N, double *Snf,
for(i=0; i<N; i++)
{
t = (double)(i)/(double)(N)*Tobs;
if(t>=tmin && t<tmax)fprintf(waveout,"%e %e\n", t-t0, ht[i]/norm);
if(t>=tmin && t<tmax)fprintf(waveout,"%.16g %.16g\n", t-t0, ht[i]/norm);
}
free(ht);
......@@ -1158,6 +1158,54 @@ void print_time_domain_waveforms(char filename[], double *h, int N, double *Snf,
fclose(waveout);
}
void print_colored_time_domain_waveforms(char filename[], double *h, int N, double *Snf, double eta, double Tobs, int imin, int imax, double tmin, double tmax)
{
/*
TODO: invFFT comes out time-reversed.
LAL uses exp(-i2pift) for their Fourier transforms
Why doesn't drealft(ty-1,N,1) work for inverse FFT?
In the meantime, hacked invFFT by changing sign of imaginary part
*/
int i;
double x,t;
FILE *waveout = fopen(filename,"w");
double *ht = malloc(N*sizeof(double));
for(i=0; i<N; i++)ht[i]=h[i];
ht[0] = 0.0;
ht[1] = 0.0;
for(i=1; i< N/2; i++)
{
if(i>imin && i<imax)
{
x = 1.0;//sqrt(Snf[i]*eta);
ht[2*i] /= x;
ht[2*i+1] /= -x;
}
else
{
ht[2*i]=ht[2*i+1]=0.0;
}
}
drealft(ht-1,N,-1);
double norm = 0.5*sqrt((double)N);
double t0 = Tobs-2.0;
for(i=0; i<N; i++)
{
t = (double)(i)/(double)(N)*Tobs;
if(t>=tmin && t<tmax)fprintf(waveout,"%.16g %.16g\n", t-t0, ht[i]/norm);
}
free(ht);
fclose(waveout);
}
void print_time_domain_hdot(char filename[], double *h, int N, double *Snf, double eta, double Tobs, int imin, int imax, double tmin, double tmax)
{
/*
......
......@@ -46,6 +46,7 @@ void print_time_domain_waveforms(char filename[], double *h, int N, double *Snf,
void print_time_domain_hdot(char filename[], double *h, int N, double *Snf, double eta, double Tobs, int imin, int imax, double tmin, double tmax);
void print_frequency_domain_waveforms(char filename[], double *h, int N, double *Snf, double Tobs, double eta, int imin, int imax);
void print_frequency_domain_data(char filename[], double *h, int N, double Tobs, int imin, int imax);
void print_colored_time_domain_waveforms(char filename[], double *h, int N, double *Snf, double eta, double Tobs, int imin, int imax, double tmin, double tmax);
void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *prior, ProcessParamsTable *commandLine);
void parse_glitch_parameters(struct Data *data, struct Model *model, FILE **paramsfile, double **grec);
......
......@@ -1226,12 +1226,18 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
{
sprintf(filename,"waveforms/%s_glitch_%d.dat.%i", modelname, frame,ifo);
print_time_domain_waveforms(filename, model[chain->index[0]]->glitch[ifo]->templates, N, model[chain->index[0]]->Snf[ifo],model[chain->index[0]]->Sn[ifo],data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
sprintf(filename,"waveforms/%s_colored_glitch_%d.dat.%i", modelname, frame,ifo);
print_colored_time_domain_waveforms(filename, model[chain->index[0]]->glitch[ifo]->templates, N, model[chain->index[0]]->Snf[ifo],model[chain->index[0]]->Sn[ifo],data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
}
if(data->signalFlag)
{
sprintf(filename,"waveforms/%s_wave_%d.dat.%i", modelname, frame, ifo);
print_time_domain_waveforms(filename, model[chain->index[0]]->response[ifo], N, model[chain->index[0]]->Snf[ifo],model[chain->index[0]]->Sn[ifo],data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
sprintf(filename,"waveforms/%s_colored_wave_%d.dat.%i", modelname, frame, ifo);
print_colored_time_domain_waveforms(filename, model[chain->index[0]]->response[ifo], N, model[chain->index[0]]->Snf[ifo],model[chain->index[0]]->Sn[ifo],data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
}
if(data->bayesLineFlag)
......
......@@ -60,7 +60,7 @@ void set_bayesline_priors(char *channel, struct BayesLineParams *bayesline, doub
else
{
bayesline->priors->SAmin = 5.0e-48*Tobs/4.0;
bayesline->priors->SAmax = 5.0e-42*Tobs/4.0;
bayesline->priors->SAmax = 5.0e-36*Tobs/4.0;
bayesline->priors->LQmin = 1.0e1;
bayesline->priors->LQmax = 1.0e7;
bayesline->priors->LAmin = 1.0e-45;
......@@ -68,7 +68,13 @@ void set_bayesline_priors(char *channel, struct BayesLineParams *bayesline, doub
}
}
int checktime(double *p, double **r,double Tobs)
{
//p[0]=t0
//r[0][0]=tmin
//r[0][1]=tmax
}
int checkrange(double *p, double **r)
{
int k;
......
......@@ -503,11 +503,11 @@ void intrinsic_fisher_update(double *params, double *dparams, double *Snf, doubl
double sqrt3=1.73205080756888;
double invSNR = 1./SNR;
dparams[0] = 1.0/(LAL_TWOPI*params[1]*SNR);
dparams[1] = 2.0*params[1]/(params[2]*SNR);
dparams[2] = 2.0*params[2]/(sqrt3*SNR);
dparams[3] = params[3]*invSNR;
dparams[4] = 1.0*invSNR;
dparams[0] = 1.0/(LAL_TWOPI*params[1]*SNR); //sqrt(1/t0t0)
dparams[1] = 2.0*params[1]/(params[2]*SNR); //sqrt(1/f0f0)
dparams[2] = 2.0*params[2]/(sqrt3*SNR); //sqrt(1/QQ)
dparams[3] = params[3]*invSNR; //sqrt(1/AA)
dparams[4] = 1.0*invSNR; //sqrt(1/phiphi)
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment