Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALApps 10.1.0.1-b246709
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
getresp.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2007 Jolien Creighton
3*
4* This program is free software; you can redistribute it and/or modify
5* it under the terms of the GNU General Public License as published by
6* the Free Software Foundation; either version 2 of the License, or
7* (at your option) any later version.
8*
9* This program is distributed in the hope that it will be useful,
10* but WITHOUT ANY WARRANTY; without even the implied warranty of
11* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12* GNU General Public License for more details.
13*
14* You should have received a copy of the GNU General Public License
15* along with with program; see the file COPYING. If not, write to the
16* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17* MA 02110-1301 USA
18*/
19
20#include "config.h"
21
22#include <math.h>
23#include <string.h>
24
25#include <lal/LALStdlib.h>
26#include <lal/LALStdio.h>
27#include <lal/AVFactories.h>
28#include <lal/LALCache.h>
29#include <lal/FrequencySeries.h>
30#include <lal/LALFrStream.h>
31#include <lal/Calibration.h>
32#include <lal/Units.h>
33#include <lal/LIGOMetadataRingdownUtils.h>
34
35#include <LALAppsVCSInfo.h>
36#include "getresp.h"
37#include "errutil.h"
38#include "gpstime.h"
39
40/* a few useful constants */
41static LALUnit strainPerCount = {0,{0,0,0,0,0,1,-1},{0,0,0,0,0,0,0}};
42
43
44/* routine read a response function or construct an impulse response function */
46 const char *cacheName,
47 const char *ifoName,
48 LIGOTimeGPS *epoch,
49 REAL8 dataDuration,
50 REAL8 dataSampleRate,
51 REAL4 responseScale,
52 int strainData,
53 const char *channel_name
54 )
55{
57
58 if ( strainData )
59 response = get_impulse_response( ifoName, epoch, dataDuration,
60 dataSampleRate, responseScale );
61 else
62 response = get_frame_response( cacheName, ifoName, epoch, dataDuration,
63 dataSampleRate, responseScale, channel_name );
64
65 return response;
66}
67
68
69/* routine to construct an impulse (uniform in frequency) response function */
71 const char *ifoName,
72 LIGOTimeGPS *epoch,
73 REAL8 dataDuration,
74 REAL8 dataSampleRate,
75 REAL4 responseScale
76 )
77{
79 UINT4 npoints;
80 UINT4 k;
81
82 /* create frequency series */
83 response = LALCalloc( 1, sizeof( *response ) );
84 if ( ! response )
85 return NULL;
86
87 /* allocate data memory and set metadata */
88 npoints = floor( dataDuration * dataSampleRate + 0.5 ); /* round */
89 snprintf( response->name, sizeof( response->name ),
90 "%s:CAL-RESPONSE", ifoName );
91 response->deltaF = 1.0/dataDuration;
92 response->epoch = *epoch;
93 response->sampleUnits = strainPerCount;
94 response->data = XLALCreateCOMPLEX8Vector( npoints/2 + 1 );
95
96 /* uniform response function */
97 for ( k = 0; k < response->data->length; ++k )
98 {
99 response->data->data[k] = crectf( responseScale, 0.0 );
100 }
101
102 return response;
103}
104
105
106/* routine to construct a response function from calibration frame files */
108 const char *cacheName,
109 const char *ifoName,
110 LIGOTimeGPS *epoch,
111 REAL8 dataDuration,
112 REAL8 dataSampleRate,
113 REAL4 responseScale,
114 const char *channel_name
115 )
116{
118 LALCache *cache = NULL;
119 COMPLEX8FrequencySeries *response;
120 COMPLEX8FrequencySeries *refResponse;
121 COMPLEX8FrequencySeries *refSensing;
123 COMPLEX8TimeSeries *alphabeta;
124 CalibrationFunctions calfuncs;
126 char ifo[3];
127 UINT4 npoints;
128 UINT4 k;
129 memset(&calfacts, 0, sizeof(calfacts));
130
131 /* create frequency series */
132 npoints = floor( dataDuration * dataSampleRate + 0.5 ); /* round */
133 response = XLALCreateCOMPLEX8FrequencySeries(channel_name, epoch, 0.0, 1.0 / dataDuration, &strainPerCount, npoints / 2 + 1);
134
135 if(!response)
136 {
138 return NULL;
139 }
140
141 verbose("obtaining calibration information from cache file %s\n", cacheName);
142
143 /* create cache from cachefile name */
144 cache = XLALCacheImport( cacheName );
145
146 /* get reference functions and factors */
147 refResponse = get_reference_response_function( cache, ifoName );
148 refSensing = get_reference_sensing_function( cache, ifoName );
150 alphabeta = get_open_loop_gain_factor( cache, epoch, ifoName );
151
152 /* done with cache... get rid of it */
154
155 /* the ifo is two characters long */
156 strncpy( ifo, ifoName, 2 );
157 ifo[2] = 0;
158
159 /* setup calfuncs structure needed by LALUpdateCalibration */
160 calfuncs.responseFunction = refResponse;
161 calfuncs.sensingFunction = refSensing;
162
163 /* setup calfacts structure needed by LALUpdateCalibration */
164 calfacts.openLoopFactor = alphabeta;
165 calfacts.sensingFactor = alpha;
166 calfacts.epoch = alpha->epoch;
167 calfacts.ifo = ifo;
168 ns_to_epoch( &calfacts.duration, sec_to_ns( alpha->deltaT ) );
169
170 /* update reference calibration to current calibration */
171 verbose( "updating reference calibration to epoch %d.%09d\n",
172 calfacts.epoch.gpsSeconds, calfacts.epoch.gpsNanoSeconds );
173 LAL_CALL( LALUpdateCalibration( &status, &calfuncs, &calfuncs, &calfacts ),
174 &status );
175 /* these are the values of alpha and alphabeta that were used */
176 verbose( "calibrating with alpha=%g and alphabeta=%g\n", crealf(calfacts.alpha),
177 crealf(calfacts.alphabeta) );
178
179 /* convert response so that it has the correct resolution */
180 LAL_CALL( LALResponseConvert( &status, response, refResponse ), &status );
181
182 /* scale response function */
183 for ( k = 0; k < response->data->length; ++k )
184 {
185 response->data->data[k] *= ((REAL4) responseScale);
186 }
187
188 /* cleanup memory in reference functions */
189 XLALDestroyCOMPLEX8Vector( alphabeta->data );
190 LALFree( alphabeta );
192 LALFree( alpha );
193 XLALDestroyCOMPLEX8Vector( refSensing->data );
194 LALFree( refSensing );
195 XLALDestroyCOMPLEX8Vector( refResponse->data );
196 LALFree( refResponse );
197
198 return response;
199}
200
201
202/* routine to read the reference response function from a frame file cache */
204 const char *ifoName )
205{
206 COMPLEX8FrequencySeries *refResponse;
207 LALCache *refCache = NULL;
208 LALFrStream *stream = NULL;
209
210 /* create frequency series */
211 refResponse = LALCalloc( 1, sizeof( *refResponse ) );
212 if ( ! refResponse )
213 return NULL;
214
215 snprintf( refResponse->name, sizeof( refResponse->name ),
216 "%s:CAL-RESPONSE", ifoName );
217
218 /* make a new cache with only the CAL_REF frames and open a frame stream */
219 refCache = XLALCacheDuplicate( calCache );
220 XLALCacheSieve( refCache, 0, 0, NULL, "CAL_REF", NULL );
221 stream = XLALFrStreamCacheOpen( refCache );
222
223 /* get the frequency series from the frame stream and correct units */
224 XLALFrStreamGetCOMPLEX8FrequencySeries( refResponse, stream );
225 refResponse->sampleUnits = strainPerCount;
226
227 /* close stream and destroy cache */
228 XLALFrStreamClose( stream );
229 XLALDestroyCache( refCache );
230
231 return refResponse;
232}
233
234
235/* routine to read the reference sensing function from a frame file cache */
237 const char *ifoName )
238{
239 COMPLEX8FrequencySeries *refSensing;
240 LALCache *refCache = NULL;
241 LALFrStream *stream = NULL;
242
243 /* create frequency series */
244 refSensing = LALCalloc( 1, sizeof( *refSensing ) );
245 if ( ! refSensing )
246 return NULL;
247
248 snprintf( refSensing->name, sizeof( refSensing->name ),
249 "%s:CAL-CAV_GAIN", ifoName );
250
251 /* make a new cache with only the CAL_REF frames and open a frame stream */
252 refCache = XLALCacheDuplicate( calCache );
253 XLALCacheSieve( refCache, 0, 0, NULL, "CAL_REF", NULL );
254 stream = XLALFrStreamCacheOpen( refCache );
255
256 /* get the frequency series from the frame stream and correct units */
257 XLALFrStreamGetCOMPLEX8FrequencySeries( refSensing, stream );
258 refSensing->sampleUnits = lalDimensionlessUnit;
259
260 /* close stream and destroy cache */
261 XLALFrStreamClose( stream );
262 XLALDestroyCache( refCache );
263
264 return refSensing;
265}
266
267
268/* routine to read one cavity gain factor from a frame file cache */
270 LIGOTimeGPS *epoch, const char *ifoName )
271{
273 LALCache *facCache = NULL;
274 LALFrStream *stream = NULL;
275
276 /* create time series with one data point */
277 alpha = LALCalloc( 1, sizeof( *alpha ) );
278 alpha->data = XLALCreateCOMPLEX8Vector( 1 ); /* only 1 point */
279
280 snprintf( alpha->name, sizeof( alpha->name ), "%s:CAL-CAV_FAC", ifoName );
281
282 /* make a new cache with only the CAL_FAC frames and open a frame stream */
283 facCache = XLALCacheDuplicate( calCache );
284 XLALCacheSieve( facCache, 0, 0, NULL, "CAL_FAC", NULL );
285 stream = XLALFrStreamCacheOpen( facCache );
286
287 /* read the first point *after* the specified epoch from the frame stream
288 * and correct the units */
289 /* FIXME: should the point be *after* the specified epoch??? */
290 XLALFrStreamSeek( stream, epoch );
292 alpha->sampleUnits = lalDimensionlessUnit;
293
294 /* close stream and destroy cache */
295 XLALFrStreamClose( stream );
296 XLALDestroyCache( facCache );
297
298 return alpha;
299}
300
301
302/* routine to read one open loop gain factor from a frame file cache */
304 LIGOTimeGPS *epoch, const char *ifoName )
305{
306 COMPLEX8TimeSeries *alphabeta;
307 LALCache *facCache = NULL;
308 LALFrStream *stream = NULL;
309
310 /* create time series with one data point */
311 alphabeta = LALCalloc( 1, sizeof( *alphabeta ) );
312 alphabeta->data = XLALCreateCOMPLEX8Vector( 1 ); /* only 1 point */
313
314 snprintf( alphabeta->name, sizeof( alphabeta->name ),
315 "%s:CAL-OLOOP_FAC", ifoName );
316
317 /* make a new cache with only the CAL_FAC frames and open a frame stream */
318 facCache = XLALCacheDuplicate( calCache );
319 XLALCacheSieve( facCache, 0, 0, NULL, "CAL_FAC", NULL );
320 stream = XLALFrStreamCacheOpen( facCache );
321
322 /* read the first point *after* the specified epoch from the frame stream
323 * and correct the units */
324 /* FIXME: should the point be *after* the specified epoch??? */
325 XLALFrStreamSeek( stream, epoch );
326 XLALFrStreamGetCOMPLEX8TimeSeries( alphabeta, stream );
328
329 /* close stream and destroy cache */
330 XLALFrStreamClose( stream );
331 XLALDestroyCache( facCache );
332
333 return alphabeta;
334}
#define LAL_CALL(function, statusptr)
int k
#define LALCalloc(m, n)
#define LALFree(p)
int verbose
Definition: chirplen.c:77
COMPLEX8FrequencySeries * get_response(const char *cacheName, const char *ifoName, LIGOTimeGPS *epoch, REAL8 dataDuration, REAL8 dataSampleRate, REAL4 responseScale, int strainData, const char *channel_name)
Definition: getresp.c:45
COMPLEX8TimeSeries * get_open_loop_gain_factor(LALCache *calCache, LIGOTimeGPS *epoch, const char *ifoName)
Definition: getresp.c:303
COMPLEX8FrequencySeries * get_frame_response(const char *cacheName, const char *ifoName, LIGOTimeGPS *epoch, REAL8 dataDuration, REAL8 dataSampleRate, REAL4 responseScale, const char *channel_name)
Definition: getresp.c:107
COMPLEX8FrequencySeries * get_impulse_response(const char *ifoName, LIGOTimeGPS *epoch, REAL8 dataDuration, REAL8 dataSampleRate, REAL4 responseScale)
Definition: getresp.c:70
COMPLEX8TimeSeries * get_cavity_gain_factor(LALCache *calCache, LIGOTimeGPS *epoch, const char *ifoName)
Definition: getresp.c:269
static LALUnit strainPerCount
Definition: getresp.c:41
COMPLEX8FrequencySeries * get_reference_response_function(LALCache *calCache, const char *ifoName)
Definition: getresp.c:203
COMPLEX8FrequencySeries * get_reference_sensing_function(LALCache *calCache, const char *ifoName)
Definition: getresp.c:236
INT8 sec_to_ns(REAL8 sec)
Definition: gpstime.c:57
LIGOTimeGPS * ns_to_epoch(LIGOTimeGPS *epoch, INT8 ns)
Definition: gpstime.c:48
void LALUpdateCalibration(LALStatus *status, CalibrationFunctions *output, CalibrationFunctions *input, CalibrationUpdateParams *params)
void LALResponseConvert(LALStatus *status, COMPLEX8FrequencySeries *output, COMPLEX8FrequencySeries *input)
COMPLEX8FrequencySeries * XLALCreateCOMPLEX8FrequencySeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaF, const LALUnit *sampleUnits, size_t length)
void XLALDestroyCache(LALCache *cache)
LALCache * XLALCacheImport(const char *fname)
LALCache * XLALCacheDuplicate(const LALCache *cache)
int XLALCacheSieve(LALCache *cache, INT4 t0, INT4 t1, const char *srcregex, const char *dscregex, const char *urlregex)
double REAL8
#define XLAL_INIT_DECL(var,...)
#define crectf(re, im)
uint32_t UINT4
float REAL4
int XLALFrStreamSeek(LALFrStream *stream, const LIGOTimeGPS *epoch)
int XLALFrStreamClose(LALFrStream *stream)
LALFrStream * XLALFrStreamCacheOpen(LALCache *cache)
int XLALFrStreamGetCOMPLEX8TimeSeries(COMPLEX8TimeSeries *series, LALFrStream *stream)
int XLALFrStreamGetCOMPLEX8FrequencySeries(COMPLEX8FrequencySeries *series, LALFrStream *stream)
const LALUnit lalDimensionlessUnit
COMPLEX8Vector * XLALCreateCOMPLEX8Vector(UINT4 length)
void XLALDestroyCOMPLEX8Vector(COMPLEX8Vector *vector)
#define XLAL_ERROR_NULL(...)
XLAL_ENOMEM
char * ifo
Definition: gwf2xml.c:57
static LALStatus status
Definition: inspinj.c:552
CHAR name[LALNameLength]
COMPLEX8Sequence * data
CHAR name[LALNameLength]
COMPLEX8Sequence * data
COMPLEX8 * data
COMPLEX8FrequencySeries * responseFunction
COMPLEX8FrequencySeries * sensingFunction
COMPLEX8TimeSeries * sensingFactor
COMPLEX8TimeSeries * openLoopFactor
INT4 gpsNanoSeconds
enum @1 epoch
REAL4 alpha
Definition: tmpltbank.c:432