Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALPulsar 7.1.1.1-8a6b96f
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
SFTfunctions.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2014 Evan Goetz
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 <lal/UserInput.h>
21#include <lal/CWMakeFakeData.h>
22#include <lal/SinCosLUT.h>
23#include <lal/GenerateSpinOrbitCW.h>
24#include <lal/TimeSeries.h>
25#include <lal/Units.h>
26#include "SFTfunctions.h"
27#include "TwoSpectSpecFunc.h"
28#include "statistics.h"
29
30/**
31 * Find the SFT data specified by user input
32 * \param [in] params Pointer to the UserInput_t
33 * \return Pointer to SFTCatalog
34 */
36{
37
39
40 fprintf( LOG, "Finding SFTs... " );
41 fprintf( stderr, "Finding SFTs... " );
42
43 //Set the start and end times in the LIGO GPS format
47 XLALGPSSetREAL8( &end, params->t0 + params->Tobs - params->SFToverlap );
49
50 //Setup the constraints
51 SFTConstraints XLAL_INIT_DECL( constraints );
52 if ( params->IFO->length == 1 ) {
53 constraints.detector = params->IFO->data[0];
54 }
55 constraints.minStartTime = &start;
56 constraints.maxStartTime = &end;
57
58 //Find SFT files
59 SFTCatalog *catalog = NULL;
60 XLAL_CHECK_NULL( ( catalog = XLALSFTdataFind( params->inputSFTs, &constraints ) ) != NULL, XLAL_EFUNC );
61
62 fprintf( LOG, "done\n" );
63 fprintf( stderr, "done\n" );
64
65 return catalog;
66
67} // findSFTdata()
68
69/**
70 * Extract the SFT coefficients from the band of interest
71 * \param [in] catalog Pointer to the SFTCatalog
72 * \param [in] minfbin Frequency value of the minimum frequency bin
73 * \param [in] maxfbin Frequency value of the maximum frequency bin
74 * \return Pointer to MultiSFTVector
75 */
76MultiSFTVector *extractSFTband( const SFTCatalog *catalog, const REAL8 minfbin, const REAL8 maxfbin )
77{
78
79 XLAL_CHECK_NULL( catalog != NULL, XLAL_EINVAL );
80
81 fprintf( LOG, "Extracting band from SFTs... " );
82 fprintf( stderr, "Extracting band from SFTs... " );
83
84 //Extract the data
85 MultiSFTVector *sftvector = NULL;
86 XLAL_CHECK_NULL( ( sftvector = XLALLoadMultiSFTs( catalog, minfbin, maxfbin ) ) != NULL, XLAL_EFUNC );
87
88 fprintf( LOG, "done\n" );
89 fprintf( stderr, "done\n" );
90
91 return sftvector;
92
93} // extractSFTband()
94
95/**
96 * Get a MultiSFTVector from the user-input values
97 * \param [in] params Pointer to UserInput_t
98 * \param [in] minfbin Frequency value of the minimum frequency bin
99 * \param [in] maxfbin Frequency value of the maximum frequency bin
100 * \return Pointer to MultiSFTVector
101 */
103{
104 //Get catalog of SFTs
105 SFTCatalog *catalog = NULL;
106 XLAL_CHECK_NULL( ( catalog = findSFTdata( params ) ) != NULL, XLAL_EFUNC );
107
108 MultiSFTVector *sftvector = NULL;
109 if ( params->markBadSFTs && !params->signalOnly ) {
110 //Extract band to get a MultiSFTVector
111 MultiSFTVector *tmpsftvector = NULL;
112 XLAL_CHECK_NULL( ( tmpsftvector = extractSFTband( catalog, minfbin, maxfbin ) ) != NULL, XLAL_EFUNC );
114
115 //Get the timestamps of the SFTs applying the KS/Kuipers tests if desired
118
119 //Get the SFT subset
121
122 XLALDestroyMultiSFTVector( tmpsftvector );
124 } else {
125 XLAL_CHECK_NULL( ( sftvector = extractSFTband( catalog, minfbin, maxfbin ) ) != NULL, XLAL_EFUNC );
127 }
128
129 XLALDestroySFTCatalog( catalog );
130
131 return sftvector;
132} // getMultiSFTVector()
133
134/**
135 * Add SFTs together from a MultiSFTVector
136 * \param [in] multiSFTvector Pointer to a MultiSFTVector containing the SFT data
137 * \param [in] multissb Pointer to a MultiSSBtimes structure
138 * \param [in] multiAMcoeffs Pointer to a MultiAMCoeffs structure
139 * \param [in] jointTimestamps Pointer to a LIGOTimeGPSVector of joint SFT times from all detectors
140 * \param [in] backgroundRatio Pointer to a REAL4VectorAlignedArray of the running means of each IFO divided by the running mean of IFO_0
141 * \param [in] cosiSign Value of 1, 0, -1 where 1 means average cosi over [0,1], 0 means average over [-1,1], and -1 means average over [-1,0]
142 * \param [in] NSparams Values for a set of assumed NS spin/binary motion parameters
143 * \param [in] params Pointer to UserInput_t
144 * \param [out] backgroundScaling Pointer to REAL4VectorAligned of background scaling values
145 * \return REAL4VectorAligned of powers of the coherently combined SFTs
146 */
147REAL4VectorAligned *coherentlyAddSFTs( const MultiSFTVector *multiSFTvector, const MultiSSBtimes *multissb, const MultiAMCoeffs *multiAMcoeffs, const LIGOTimeGPSVector *jointTimestamps, const REAL4VectorAlignedArray *backgroundRatio, const INT4 cosiSign, const assumeNSparams *NSparams, const UserInput_t *params, REAL4VectorAligned *backgroundScaling )
148{
149
150 //Sanity check the values
151 XLAL_CHECK_NULL( multiSFTvector != NULL && multissb != NULL && multiAMcoeffs != NULL && jointTimestamps != NULL && params != NULL, XLAL_EINVAL );
152
153 fprintf( stderr, "Coherently adding SFT data... " );
154
155 //parse noise list
156 MultiNoiseFloor multiNoiseFloor;
157 XLAL_CHECK_NULL( XLALParseMultiNoiseFloor( &multiNoiseFloor, params->avesqrtSh, params->IFO->length ) == XLAL_SUCCESS, XLAL_EFUNC );
158
159 //initialize backgroundScaling and normalizationScaling to values of 0
160 memset( backgroundScaling->data, 0, sizeof( REAL4 )*backgroundScaling->length );
161
162 //get total number of SFTs possible, number of frequency bins in each FFT of the background, and number of original SFT bins to skip because of background calculation
163 UINT4 numFbinsInBackground = multiSFTvector->data[0]->data[0].data->length - 20 - ( params->blksize - 1 );
164 UINT4 numSFTbins2skip = ( multiSFTvector->data[0]->data[0].data->length - numFbinsInBackground ) / 2;
165
166 //Allocate the combined SFTs SFTVector
167 SFTVector *combinedSFTs = NULL;
168 XLAL_CHECK_NULL( ( combinedSFTs = XLALCreateSFTVector( jointTimestamps->length, 0 ) ) != NULL, XLAL_EFUNC );
169
170 //Create an INT4Vector to determine at which vector we are using in the multiSFTvector
171 INT4Vector *whichSFTinMultiSFTvector = NULL, *whichIFOsToBeUsed = NULL;
172 XLAL_CHECK_NULL( ( whichSFTinMultiSFTvector = XLALCreateINT4Vector( multiSFTvector->length ) ) != NULL, XLAL_EFUNC );
173 XLAL_CHECK_NULL( ( whichIFOsToBeUsed = XLALCreateINT4Vector( multiSFTvector->length ) ) != NULL, XLAL_EFUNC );
174 memset( whichSFTinMultiSFTvector->data, 0, whichSFTinMultiSFTvector->length * sizeof( INT4 ) ); //All index values are set to zero at the beginning
175
176 //Pre-compute the different sin(2*psi) and cos(2*psi) values
177 REAL4VectorAligned *twoPsiVec = NULL, *sin2psiVec = NULL, *cos2psiVec = NULL, *Fplus0s = NULL, *Fcross0s = NULL, *FplusXs = NULL, *FcrossXs = NULL, *aVals = NULL, *bVals = NULL;
178 XLAL_CHECK_NULL( ( twoPsiVec = XLALCreateREAL4VectorAligned( 16, 32 ) ) != NULL, XLAL_EFUNC );
179 XLAL_CHECK_NULL( ( sin2psiVec = XLALCreateREAL4VectorAligned( twoPsiVec->length, 32 ) ) != NULL, XLAL_EFUNC );
180 XLAL_CHECK_NULL( ( cos2psiVec = XLALCreateREAL4VectorAligned( twoPsiVec->length, 32 ) ) != NULL, XLAL_EFUNC );
181 for ( UINT4 ii = 0; ii < twoPsiVec->length; ii++ ) {
182 twoPsiVec->data[ii] = 0.03125 * ii;
183 }
184 XLAL_CHECK_NULL( XLALVectorSinCos2PiREAL4( sin2psiVec->data, cos2psiVec->data, twoPsiVec->data, twoPsiVec->length ) == XLAL_SUCCESS, XLAL_EFUNC );
185 XLAL_CHECK_NULL( ( Fplus0s = XLALCreateREAL4VectorAligned( twoPsiVec->length, 32 ) ) != NULL, XLAL_EFUNC );
186 XLAL_CHECK_NULL( ( Fcross0s = XLALCreateREAL4VectorAligned( twoPsiVec->length, 32 ) ) != NULL, XLAL_EFUNC );
187 XLAL_CHECK_NULL( ( FplusXs = XLALCreateREAL4VectorAligned( twoPsiVec->length, 32 ) ) != NULL, XLAL_EFUNC );
188 XLAL_CHECK_NULL( ( FcrossXs = XLALCreateREAL4VectorAligned( twoPsiVec->length, 32 ) ) != NULL, XLAL_EFUNC );
189 XLAL_CHECK_NULL( ( aVals = XLALCreateREAL4VectorAligned( twoPsiVec->length, 32 ) ) != NULL, XLAL_EFUNC );
190 XLAL_CHECK_NULL( ( bVals = XLALCreateREAL4VectorAligned( twoPsiVec->length, 32 ) ) != NULL, XLAL_EFUNC );
191
192 //Pre-compute the GW frequencies and the bin number of the GW frequencies and frequencies of the SFT bins
193 alignedREAL8Vector *GWfrequencies = NULL, *GWfrequencyBins = NULL, *SFTfrequencies = NULL, *scaledSFTfrequencies = NULL;
194 XLAL_CHECK_NULL( ( GWfrequencies = createAlignedREAL8Vector( multiSFTvector->data[0]->data[0].data->length, 32 ) ) != NULL, XLAL_EFUNC );
195 XLAL_CHECK_NULL( ( GWfrequencyBins = createAlignedREAL8Vector( GWfrequencies->length, 32 ) ) != NULL, XLAL_EFUNC );
196 XLAL_CHECK_NULL( ( SFTfrequencies = createAlignedREAL8Vector( GWfrequencies->length, 32 ) ) != NULL, XLAL_EFUNC );
197 SFTfrequencies->data[0] = multiSFTvector->data[0]->data[0].f0;
198 for ( UINT4 ii = 1; ii < SFTfrequencies->length; ii++ ) {
199 SFTfrequencies->data[ii] = SFTfrequencies->data[ii - 1] + multiSFTvector->data[0]->data[0].deltaF;
200 }
201 if ( NSparams->assumeNSGWfreq == NULL ) {
202 memcpy( GWfrequencies->data, SFTfrequencies->data, sizeof( REAL8 )*SFTfrequencies->length );
203 } else {
204 for ( UINT4 ii = 0; ii < GWfrequencies->length; ii++ ) {
205 GWfrequencies->data[ii] = *( NSparams->assumeNSGWfreq );
206 }
207 XLAL_CHECK_NULL( ( scaledSFTfrequencies = createAlignedREAL8Vector( SFTfrequencies->length, 32 ) ) != NULL, XLAL_EFUNC );
208 }
209 XLAL_CHECK_NULL( VectorScaleREAL8( GWfrequencyBins, GWfrequencies, params->Tsft, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
210
211 //Pre-allocate delta values for IFO_0 and IFO_X, also 2*pi*f_k*tau
212 alignedREAL8Vector *delta0vals = NULL, *deltaXvals = NULL, *delta0valsSubset = NULL, *deltaXvalsSubset = NULL, *floorDelta0valsSubset = NULL, *floorDeltaXvalsSubset = NULL, *diffFloorDeltaValsSubset = NULL, *roundDelta0valsSubset = NULL, *roundDeltaXvalsSubset = NULL, *diffRoundDeltaValsSubset = NULL, *TwoPiGWfrequenciesTau = NULL, *absDiffFloorDeltaValsSubset = NULL;
213 REAL4VectorAligned *TwoPiGWfrequenciesTau_f = NULL, *sinTwoPiGWfrequenciesTau = NULL, *cosTwoPiGWfrequenciesTau = NULL;
214 XLAL_CHECK_NULL( ( delta0vals = createAlignedREAL8Vector( GWfrequencyBins->length, 32 ) ) != NULL, XLAL_EFUNC );
215 XLAL_CHECK_NULL( ( deltaXvals = createAlignedREAL8Vector( GWfrequencyBins->length, 32 ) ) != NULL, XLAL_EFUNC );
216 XLAL_CHECK_NULL( ( delta0valsSubset = createAlignedREAL8Vector( GWfrequencyBins->length - 20, 32 ) ) != NULL, XLAL_EFUNC );
217 XLAL_CHECK_NULL( ( deltaXvalsSubset = createAlignedREAL8Vector( delta0valsSubset->length, 32 ) ) != NULL, XLAL_EFUNC );
218 XLAL_CHECK_NULL( ( floorDelta0valsSubset = createAlignedREAL8Vector( delta0valsSubset->length, 32 ) ) != NULL, XLAL_EFUNC );
219 XLAL_CHECK_NULL( ( floorDeltaXvalsSubset = createAlignedREAL8Vector( delta0valsSubset->length, 32 ) ) != NULL, XLAL_EFUNC );
220 XLAL_CHECK_NULL( ( roundDelta0valsSubset = createAlignedREAL8Vector( delta0valsSubset->length, 32 ) ) != NULL, XLAL_EFUNC );
221 XLAL_CHECK_NULL( ( roundDeltaXvalsSubset = createAlignedREAL8Vector( delta0valsSubset->length, 32 ) ) != NULL, XLAL_EFUNC );
222 XLAL_CHECK_NULL( ( diffFloorDeltaValsSubset = createAlignedREAL8Vector( delta0valsSubset->length, 32 ) ) != NULL, XLAL_EFUNC );
223 XLAL_CHECK_NULL( ( absDiffFloorDeltaValsSubset = createAlignedREAL8Vector( delta0valsSubset->length, 32 ) ) != NULL, XLAL_EFUNC );
224 XLAL_CHECK_NULL( ( diffRoundDeltaValsSubset = createAlignedREAL8Vector( delta0valsSubset->length, 32 ) ) != NULL, XLAL_EFUNC );
225 XLAL_CHECK_NULL( ( TwoPiGWfrequenciesTau = createAlignedREAL8Vector( GWfrequencyBins->length, 32 ) ) != NULL, XLAL_EFUNC );
226 XLAL_CHECK_NULL( ( TwoPiGWfrequenciesTau_f = XLALCreateREAL4VectorAligned( GWfrequencyBins->length, 32 ) ) != NULL, XLAL_EFUNC );
227 XLAL_CHECK_NULL( ( sinTwoPiGWfrequenciesTau = XLALCreateREAL4VectorAligned( GWfrequencyBins->length, 32 ) ) != NULL, XLAL_EFUNC );
228 XLAL_CHECK_NULL( ( cosTwoPiGWfrequenciesTau = XLALCreateREAL4VectorAligned( GWfrequencyBins->length, 32 ) ) != NULL, XLAL_EFUNC );
229
230 INT4Vector *shiftVector = NULL;
231 XLAL_CHECK_NULL( ( shiftVector = XLALCreateINT4Vector( delta0valsSubset->length ) ) != NULL, XLAL_EFUNC );
232
233 //Mid-point of the SFT
234 REAL8 Tmid = 0.5 * params->Tsft;
235
236 //Pre-compute 2*pi*{[DeltaTs(jj) - Tmid*Tdots(jj)]-[DeltaTs(0) - Tmid*Tdots(0)]}=twoPiTauVals and [Tdots(jj) - 1]=TdotMinus1s for jj=0,...,N-1
237 //Note that twoPiTauVals->data[0] is unused later in the code; it is only necessary for the intermediate calculation
238 alignedREAL8VectorArray *twoPiTauVals = NULL, *TdotMinus1s = NULL;
239 alignedREAL8Vector *Tdots = NULL, *DeltaTs = NULL;
240 XLAL_CHECK_NULL( ( twoPiTauVals = createAlignedREAL8VectorArray( multissb->length, multissb->data[0]->DeltaT->length, 32 ) ) != NULL, XLAL_EFUNC );
241 XLAL_CHECK_NULL( ( TdotMinus1s = createAlignedREAL8VectorArray( multissb->length, multissb->data[0]->Tdot->length, 32 ) ) != NULL, XLAL_EFUNC );
242 XLAL_CHECK_NULL( ( Tdots = createAlignedREAL8Vector( multissb->data[0]->Tdot->length, 32 ) ) != NULL, XLAL_EFUNC );
243 XLAL_CHECK_NULL( ( DeltaTs = createAlignedREAL8Vector( multissb->data[0]->DeltaT->length, 32 ) ) != NULL, XLAL_EFUNC );
244 for ( UINT4 ii = 0; ii < twoPiTauVals->length; ii++ ) {
245 memcpy( Tdots->data, multissb->data[ii]->Tdot->data, sizeof( REAL8 )*multissb->data[ii]->Tdot->length );
246 memcpy( DeltaTs->data, multissb->data[ii]->DeltaT->data, sizeof( REAL8 )*multissb->data[ii]->DeltaT->length );
247 XLAL_CHECK_NULL( VectorScaleREAL8( twoPiTauVals->data[ii], Tdots, -Tmid, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
248 XLAL_CHECK_NULL( VectorAddREAL8( twoPiTauVals->data[ii], twoPiTauVals->data[ii], DeltaTs, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
249 XLAL_CHECK_NULL( VectorShiftREAL8( TdotMinus1s->data[ii], Tdots, -1.0, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
250 }
251 for ( UINT4 ii = twoPiTauVals->length; ii > 0; ii-- ) {
252 XLAL_CHECK_NULL( VectorSubtractREAL8( twoPiTauVals->data[ii - 1], twoPiTauVals->data[ii - 1], twoPiTauVals->data[0], params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
253 }
254 for ( UINT4 ii = 1; ii < twoPiTauVals->length; ii++ ) {
255 XLAL_CHECK_NULL( VectorScaleREAL8( twoPiTauVals->data[ii], twoPiTauVals->data[ii], LAL_TWOPI, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
256 }
258 destroyAlignedREAL8Vector( DeltaTs );
259
260 //Pre-compute cosi and (cosi*cosi+1)/2 values
261 REAL4VectorAligned *cosiVals = NULL, *onePlusCosiSqOver2Vals = NULL;
262 XLAL_CHECK_NULL( ( cosiVals = XLALCreateREAL4VectorAligned( 21, 32 ) ) != NULL, XLAL_EFUNC );
263 XLAL_CHECK_NULL( ( onePlusCosiSqOver2Vals = XLALCreateREAL4VectorAligned( 21, 32 ) ) != NULL, XLAL_EFUNC );
264 if ( cosiSign == 1 ) for ( UINT4 ii = 0; ii < cosiVals->length; ii++ ) {
265 cosiVals->data[ii] = 1.0 - 0.05 * ii;
266 } else if ( cosiSign == -1 ) for ( UINT4 ii = 0; ii < cosiVals->length; ii++ ) {
267 cosiVals->data[ii] = -0.05 * ii;
268 } else for ( UINT4 ii = 0; ii < cosiVals->length; ii++ ) {
269 cosiVals->data[ii] = 1.0 - 0.05 * ii * 2.0;
270 }
271 XLAL_CHECK_NULL( XLALVectorMultiplyREAL4( onePlusCosiSqOver2Vals->data, cosiVals->data, cosiVals->data, cosiVals->length ) == XLAL_SUCCESS, XLAL_EFUNC );
272 XLAL_CHECK_NULL( XLALVectorShiftREAL4( onePlusCosiSqOver2Vals->data, ( REAL4 )1.0, onePlusCosiSqOver2Vals->data, onePlusCosiSqOver2Vals->length ) == XLAL_SUCCESS, XLAL_EFUNC );
273 XLAL_CHECK_NULL( XLALVectorScaleREAL4( onePlusCosiSqOver2Vals->data, ( REAL4 )0.5, onePlusCosiSqOver2Vals->data, onePlusCosiSqOver2Vals->length ) == XLAL_SUCCESS, XLAL_EFUNC );
274
275 //Pre-allocate Aplus and Across
276 REAL4VectorAligned *Aplus0s = NULL, *AplusXs = NULL, *Across0s = NULL, *AcrossXs = NULL;
277 XLAL_CHECK_NULL( ( Aplus0s = XLALCreateREAL4VectorAligned( twoPsiVec->length, 32 ) ) != NULL, XLAL_EFUNC );
278 XLAL_CHECK_NULL( ( Across0s = XLALCreateREAL4VectorAligned( twoPsiVec->length, 32 ) ) != NULL, XLAL_EFUNC );
279 XLAL_CHECK_NULL( ( AplusXs = XLALCreateREAL4VectorAligned( twoPsiVec->length, 32 ) ) != NULL, XLAL_EFUNC );
280 XLAL_CHECK_NULL( ( AcrossXs = XLALCreateREAL4VectorAligned( twoPsiVec->length, 32 ) ) != NULL, XLAL_EFUNC );
281
282 //Pre-allocate DirichletScaling0 and DirichletScaling1
283 alignedREAL8Vector *DirichletScaling0 = NULL, *DirichletScalingX = NULL, *scaling = NULL;
284 XLAL_CHECK_NULL( ( DirichletScaling0 = createAlignedREAL8Vector( delta0valsSubset->length, 32 ) ) != NULL, XLAL_EFUNC );
285 XLAL_CHECK_NULL( ( DirichletScalingX = createAlignedREAL8Vector( delta0valsSubset->length, 32 ) ) != NULL, XLAL_EFUNC );
286 XLAL_CHECK_NULL( ( scaling = createAlignedREAL8Vector( delta0valsSubset->length, 32 ) ) != NULL, XLAL_EFUNC );
287
288 //Pre-allocate the Dirichlet ratio vector and its magnitude
289 COMPLEX8Vector *Dratio = NULL;
290 REAL4VectorAligned *cabsDratio = NULL;
291 XLAL_CHECK_NULL( ( Dratio = XLALCreateCOMPLEX8Vector( delta0valsSubset->length ) ) != NULL, XLAL_EFUNC );
292 XLAL_CHECK_NULL( ( cabsDratio = XLALCreateREAL4VectorAligned( Dratio->length, 32 ) ) != NULL, XLAL_EFUNC );
293
294 //Determine new band size
295 REAL8 newFmin = SFTfrequencies->data[10];
296 REAL8 newBand = round( ( SFTfrequencies->data[SFTfrequencies->length - 11] - SFTfrequencies->data[10] ) * params->Tsft ) / params->Tsft;
297
298 //If using knowledge of a source, then we simulate it using pulsar simulation routines in order to determine the SSB GW frequency
299 REAL4TimeSeries *GWfreqInSSB = NULL;
300 if ( NSparams->assumeNSGWfreq != NULL ) {
301 PulsarParams XLAL_INIT_DECL( pulsarParams );
302 LIGOTimeGPS SSBepoch = multissb->data[0]->refTime, spinEpoch;
303 XLALGPSAdd( &SSBepoch, multissb->data[0]->DeltaT->data[0] - params->SFToverlap );
304 pulsarParams.Amp.aPlus = 0.0; //arbitrary, we only care about frequency
305 pulsarParams.Amp.aCross = 0.0; //arbitrary, we only care about frequency
306 pulsarParams.Amp.psi = 0.0; //arbitrary, we only care about frequency
307 pulsarParams.Amp.phi0 = 0.0; //arbitrary, we only care about frequency
308 if ( NSparams->assumeNSrefTime == NULL ) {
309 XLALGPSSetREAL8( &spinEpoch, params->t0 );
310 pulsarParams.Doppler.refTime = spinEpoch;
311 } else {
312 pulsarParams.Doppler.refTime = *( NSparams->assumeNSrefTime );
313 }
314 pulsarParams.Doppler.Alpha = NSparams->assumeNSpos.longitude;
315 pulsarParams.Doppler.Delta = NSparams->assumeNSpos.latitude;
316 pulsarParams.Doppler.fkdot[0] = *( NSparams->assumeNSGWfreq );
317 pulsarParams.Doppler.asini = *( NSparams->assumeNSasini );
318 pulsarParams.Doppler.period = *( NSparams->assumeNSorbitP );
319 pulsarParams.Doppler.ecc = 0.0;
320 pulsarParams.Doppler.tp = *( NSparams->assumeNSorbitTp );
321 pulsarParams.Doppler.argp = 0.0;
322 XLAL_CHECK_NULL( ( GWfreqInSSB = computeNSfreqTS( &pulsarParams, SSBepoch, params->Tsft, params->SFToverlap, params->Tobs ) ) != NULL, XLAL_EFUNC );
323 }
324
325 //Loop over the combinedSFTs vector to fill it with single or coherently combined SFTs
326 for ( UINT4 ii = 0; ii < combinedSFTs->length; ii++ ) {
327 //Loop over the interferometers, determining which ones to use and don't go off the end of a list
328 memset( whichIFOsToBeUsed->data, 0, whichIFOsToBeUsed->length * sizeof( INT4 ) ); //All index values are set to zero each time
329 BOOLEAN foundAnSFT = 0;
330 for ( UINT4 jj = 0; jj < multiSFTvector->length; jj++ ) {
331 if ( whichSFTinMultiSFTvector->data[jj] < ( INT4 )multiSFTvector->data[jj]->length && XLALGPSCmp( &( jointTimestamps->data[ii] ), &( multiSFTvector->data[jj]->data[whichSFTinMultiSFTvector->data[jj]].epoch ) ) == 0 ) {
332 whichIFOsToBeUsed->data[jj] = 1;
333 foundAnSFT = 1;
334 }
335 }
336 XLAL_CHECK_NULL( foundAnSFT == 1, XLAL_EFAILED );
337 REAL8 sftstart2 = XLALGPSGetREAL8( &( jointTimestamps->data[ii] ) );
339 INT4 fftnum2 = ( INT4 )round( ( sftstart2 - params->t0 ) / params->SFToverlap );
340
341 //Set backgroundScaling and normalizationScaling to 1 for this SFT, and then we will add to it as more SFTs are coherently summed
342 for ( UINT4 jj = 0; jj < numFbinsInBackground; jj++ ) {
343 backgroundScaling->data[fftnum2 * numFbinsInBackground + jj] = 1.0;
344 }
345
346 //If there is an assumed NS GW frequency, then determine the frequency as a function of time and scale the SFT frequency bins by -1/assumeNSGWfreq
347 if ( NSparams->assumeNSGWfreq != NULL ) {
348 for ( UINT4 jj = 0; jj < GWfrequencies->length; jj++ ) {
349 GWfrequencies->data[jj] = GWfreqInSSB->data->data[fftnum2];
350 }
351 XLAL_CHECK_NULL( VectorScaleREAL8( GWfrequencyBins, GWfrequencies, params->Tsft, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
352 XLAL_CHECK_NULL( VectorScaleREAL8( scaledSFTfrequencies, SFTfrequencies, -1.0 / GWfrequencies->data[0], params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
353 }
354
355 BOOLEAN createSFT = 1, computeAntenna0 = 1, computeDelta0vals = 1;
356 for ( UINT4 jj = 0; jj < multiSFTvector->length; jj++ ) {
357 if ( jj == 0 && whichIFOsToBeUsed->data[jj] == 1 ) {
358 //Copy the data from the multiSFTvector into the combinedSFTs vector, ignoring 10 bins on both ends of the SFT (newFmin and newBand)
359 SFTtype *thisSFT = &( combinedSFTs->data[ii] );
360 XLAL_CHECK_NULL( XLALExtractStrictBandFromSFT( &thisSFT, &( multiSFTvector->data[0]->data[whichSFTinMultiSFTvector->data[0]] ), newFmin, newBand ) == XLAL_SUCCESS, XLAL_EFUNC );
361 createSFT = 0;
362 whichSFTinMultiSFTvector->data[0]++;
363 } else if ( jj > 0 && whichIFOsToBeUsed->data[jj] == 1 ) {
364 //Create a copy of the SFT to be shifted since we will manipulate the SFT coefficients
365 SFTtype *sftcopySubset = NULL;
366 XLAL_CHECK_NULL( XLALExtractStrictBandFromSFT( &sftcopySubset, &( multiSFTvector->data[jj]->data[whichSFTinMultiSFTvector->data[jj]] ), newFmin, newBand ) == XLAL_SUCCESS, XLAL_EFUNC );
367 REAL8 sftstart = XLALGPSGetREAL8( &( sftcopySubset->epoch ) );
369 INT4 fftnum = ( INT4 )round( ( sftstart - params->t0 ) / params->SFToverlap );
370 XLAL_CHECK_NULL( fftnum == fftnum2, XLAL_EFAILED );
371
372 //First do time of arrival: 2*pi*f*tau where f is the GW frequency either assumed (assumeNSGWfreq) or unassumed (f = fk)
373 XLAL_CHECK_NULL( VectorScaleREAL8( TwoPiGWfrequenciesTau, GWfrequencies, twoPiTauVals->data[jj]->data[fftnum], params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
374
375 //If assumeNSpsi is not given compute Fplus and Fcross for detector ratio
376 if ( NSparams->assumeNSpsi == NULL ) {
377 if ( computeAntenna0 ) {
378 XLAL_CHECK_NULL( XLALVectorScaleREAL4( aVals->data, multiAMcoeffs->data[0]->a->data[fftnum], cos2psiVec->data, cos2psiVec->length ) == XLAL_SUCCESS, XLAL_EFUNC );
379 XLAL_CHECK_NULL( XLALVectorScaleREAL4( bVals->data, multiAMcoeffs->data[0]->b->data[fftnum], sin2psiVec->data, sin2psiVec->length ) == XLAL_SUCCESS, XLAL_EFUNC );
380 XLAL_CHECK_NULL( XLALVectorAddREAL4( Fplus0s->data, aVals->data, bVals->data, aVals->length ) == XLAL_SUCCESS, XLAL_EFUNC );
381 XLAL_CHECK_NULL( XLALVectorScaleREAL4( bVals->data, multiAMcoeffs->data[0]->b->data[fftnum], cos2psiVec->data, cos2psiVec->length ) == XLAL_SUCCESS, XLAL_EFUNC );
382 XLAL_CHECK_NULL( XLALVectorScaleREAL4( aVals->data, -multiAMcoeffs->data[0]->a->data[fftnum], sin2psiVec->data, sin2psiVec->length ) == XLAL_SUCCESS, XLAL_EFUNC );
383 XLAL_CHECK_NULL( XLALVectorAddREAL4( Fcross0s->data, bVals->data, aVals->data, bVals->length ) == XLAL_SUCCESS, XLAL_EFUNC );
384 computeAntenna0 = 0;
385 }
386 XLAL_CHECK_NULL( XLALVectorScaleREAL4( aVals->data, multiAMcoeffs->data[jj]->a->data[fftnum], cos2psiVec->data, cos2psiVec->length ) == XLAL_SUCCESS, XLAL_EFUNC );
387 XLAL_CHECK_NULL( XLALVectorScaleREAL4( bVals->data, multiAMcoeffs->data[jj]->b->data[fftnum], sin2psiVec->data, sin2psiVec->length ) == XLAL_SUCCESS, XLAL_EFUNC );
388 XLAL_CHECK_NULL( XLALVectorAddREAL4( FplusXs->data, aVals->data, bVals->data, aVals->length ) == XLAL_SUCCESS, XLAL_EFUNC );
389 XLAL_CHECK_NULL( XLALVectorScaleREAL4( bVals->data, multiAMcoeffs->data[jj]->b->data[fftnum], cos2psiVec->data, cos2psiVec->length ) == XLAL_SUCCESS, XLAL_EFUNC );
390 XLAL_CHECK_NULL( XLALVectorScaleREAL4( aVals->data, -multiAMcoeffs->data[jj]->a->data[fftnum], sin2psiVec->data, sin2psiVec->length ) == XLAL_SUCCESS, XLAL_EFUNC );
391 XLAL_CHECK_NULL( XLALVectorAddREAL4( FcrossXs->data, bVals->data, aVals->data, bVals->length ) == XLAL_SUCCESS, XLAL_EFUNC );
392 }
393
394 //Average of detector-signal phase relation, unless assumed NS orientation and GW polarization angle
395 //Need to compute (Aplus^X+i*Across^X)/(Aplus^0+i*Across^0) to determine an average value of magnitude and phase
396 //where exponent indicates different IFO <X>.
397 REAL4 detPhaseArg = 0.0, detPhaseMag = 0.0;
398 if ( NSparams->assumeNScosi == NULL && NSparams->assumeNSpsi == NULL ) {
399 BOOLEAN loopbroken = 0;
400 for ( UINT4 kk = 0; kk < cosiVals->length && !loopbroken; kk++ ) {
401 XLAL_CHECK_NULL( XLALVectorScaleREAL4( AplusXs->data, onePlusCosiSqOver2Vals->data[kk], FplusXs->data, FplusXs->length ) == XLAL_SUCCESS, XLAL_EFUNC );
402 XLAL_CHECK_NULL( XLALVectorScaleREAL4( AcrossXs->data, cosiVals->data[kk], FcrossXs->data, FcrossXs->length ) == XLAL_SUCCESS, XLAL_EFUNC );
403 XLAL_CHECK_NULL( XLALVectorScaleREAL4( Aplus0s->data, onePlusCosiSqOver2Vals->data[kk], Fplus0s->data, Fplus0s->length ) == XLAL_SUCCESS, XLAL_EFUNC );
404 XLAL_CHECK_NULL( XLALVectorScaleREAL4( Across0s->data, cosiVals->data[kk], Fcross0s->data, Fcross0s->length ) == XLAL_SUCCESS, XLAL_EFUNC );
405 for ( UINT4 ll = 0; ll < cos2psiVec->length && !loopbroken; ll++ ) {
406 COMPLEX8 complexnumerator = crectf( AplusXs->data[ll], AcrossXs->data[ll] );
407 COMPLEX8 complexdenominator = crectf( Aplus0s->data[ll], Across0s->data[ll] );
408 if ( cabsf( complexdenominator ) > 1.0e-6 ) {
409 COMPLEX8 complexval = complexnumerator / complexdenominator;
410 detPhaseMag += fminf( cabsf( complexval ), 10.0 ); //fmin here because sometimes the magnitude value is being divided by a small value and causing biases
411 detPhaseArg += ( REAL4 )gsl_sf_angle_restrict_pos( ( REAL8 )cargf( complexval ) );
412 } else {
413 loopbroken = 1;
414 detPhaseMag = 0.0;
415 detPhaseArg = 0.0;
416 }
417 }
418 }
419 detPhaseMag /= ( REAL4 )( cosiVals->length * cos2psiVec->length );
420 detPhaseArg /= ( REAL4 )( cosiVals->length * cos2psiVec->length );
421 } else if ( NSparams->assumeNScosi != NULL && NSparams->assumeNSpsi == NULL ) {
422 BOOLEAN loopbroken = 0;
423 REAL4 cosi = *( NSparams->assumeNScosi );
424 REAL4 onePlusCosiSqOver2 = 0.5 * ( 1.0 + cosi * cosi );
425 XLAL_CHECK_NULL( XLALVectorScaleREAL4( AplusXs->data, onePlusCosiSqOver2, FplusXs->data, FplusXs->length ) == XLAL_SUCCESS, XLAL_EFUNC );
426 XLAL_CHECK_NULL( XLALVectorScaleREAL4( AcrossXs->data, cosi, FcrossXs->data, FcrossXs->length ) == XLAL_SUCCESS, XLAL_EFUNC );
427 XLAL_CHECK_NULL( XLALVectorScaleREAL4( Aplus0s->data, onePlusCosiSqOver2, Fplus0s->data, Fplus0s->length ) == XLAL_SUCCESS, XLAL_EFUNC );
428 XLAL_CHECK_NULL( XLALVectorScaleREAL4( Across0s->data, cosi, Fcross0s->data, Fcross0s->length ) == XLAL_SUCCESS, XLAL_EFUNC );
429 for ( UINT4 kk = 0; kk < cos2psiVec->length && !loopbroken; kk++ ) {
430 COMPLEX16 complexnumerator = crect( AplusXs->data[kk], AcrossXs->data[kk] );
431 COMPLEX16 complexdenominator = crect( Aplus0s->data[kk], Across0s->data[kk] );
432 if ( cabs( complexdenominator ) > 1.0e-6 ) {
433 COMPLEX16 complexval = complexnumerator / complexdenominator;
434 detPhaseMag += fmin( cabs( complexval ), 10.0 ); //fmin here because sometimes the magnitude value is being divided by a small value and causing biases
435 detPhaseArg += gsl_sf_angle_restrict_pos( carg( complexval ) );
436 } else {
437 loopbroken = 1;
438 detPhaseMag = 0.0;
439 detPhaseArg = 0.0;
440 }
441 }
442 detPhaseMag /= ( REAL4 )( cos2psiVec->length );
443 detPhaseArg /= ( REAL4 )( cos2psiVec->length );
444 } else if ( NSparams->assumeNScosi == NULL && NSparams->assumeNSpsi != NULL ) {
445 REAL4 psi = *( NSparams->assumeNSpsi );
446 REAL4 sin2psi = 0.0, cos2psi = 0.0;
447 XLAL_CHECK_NULL( XLALSinCosLUT( &sin2psi, &cos2psi, 2.0 * psi ) == XLAL_SUCCESS, XLAL_EFUNC );
448 if ( sin2psi > 1.0 ) {
449 sin2psi = 1.0;
450 }
451 if ( sin2psi < -1.0 ) {
452 sin2psi = -1.0;
453 }
454 if ( cos2psi > 1.0 ) {
455 cos2psi = 1.0;
456 }
457 if ( cos2psi < -1.0 ) {
458 cos2psi = -1.0;
459 }
460 REAL4 Fplus0 = multiAMcoeffs->data[0]->a->data[fftnum] * cos2psi + multiAMcoeffs->data[0]->b->data[fftnum] * sin2psi;
461 REAL4 Fcross0 = multiAMcoeffs->data[0]->b->data[fftnum] * cos2psi - multiAMcoeffs->data[0]->a->data[fftnum] * sin2psi;
462 REAL4 FplusX = multiAMcoeffs->data[jj]->a->data[fftnum] * cos2psi + multiAMcoeffs->data[jj]->b->data[fftnum] * sin2psi;
463 REAL4 FcrossX = multiAMcoeffs->data[jj]->b->data[fftnum] * cos2psi - multiAMcoeffs->data[jj]->a->data[fftnum] * sin2psi;
464 BOOLEAN loopbroken = 0;
465 for ( UINT4 kk = 0; kk < cosiVals->length && !loopbroken; kk++ ) {
466 COMPLEX16 complexnumerator = crect( FplusX * onePlusCosiSqOver2Vals->data[kk], FcrossX * cosiVals->data[kk] );
467 COMPLEX16 complexdenominator = crect( Fplus0 * onePlusCosiSqOver2Vals->data[kk], Fcross0 * cosiVals->data[kk] );
468 if ( cabs( complexdenominator ) > 1.0e-6 ) {
469 COMPLEX16 complexval = complexnumerator / complexdenominator;
470 detPhaseMag += fmin( cabs( complexval ), 10.0 ); //fmin here because sometimes the magnitude value is being divided by a small value and causing biases
471 detPhaseArg += gsl_sf_angle_restrict_pos( carg( complexval ) );
472 } else {
473 loopbroken = 1;
474 detPhaseMag = 0.0;
475 detPhaseArg = 0.0;
476 }
477 }
478 detPhaseMag /= ( REAL4 )( cosiVals->length );
479 detPhaseArg /= ( REAL4 )( cosiVals->length );
480 } else {
481 REAL4 psi = *( NSparams->assumeNSpsi );
482 REAL8 sin2psi = sin( 2.0 * psi ), cos2psi = cos( 2.0 * psi );
483 REAL8 Fplus0 = multiAMcoeffs->data[0]->a->data[fftnum] * cos2psi + multiAMcoeffs->data[0]->b->data[fftnum] * sin2psi;
484 REAL8 Fcross0 = multiAMcoeffs->data[0]->b->data[fftnum] * cos2psi - multiAMcoeffs->data[0]->a->data[fftnum] * sin2psi;
485 REAL8 FplusX = multiAMcoeffs->data[jj]->a->data[fftnum] * cos2psi + multiAMcoeffs->data[jj]->b->data[fftnum] * sin2psi;
486 REAL8 FcrossX = multiAMcoeffs->data[jj]->b->data[fftnum] * cos2psi - multiAMcoeffs->data[jj]->a->data[fftnum] * sin2psi;
487 REAL4 cosi = *( NSparams->assumeNScosi );
488 REAL8 onePlusCosiSqOver2 = 0.5 * ( 1.0 + cosi * cosi );
489 COMPLEX16 complexnumerator = crect( FplusX * onePlusCosiSqOver2, FcrossX * cosi );
490 COMPLEX16 complexdenominator = crect( Fplus0 * onePlusCosiSqOver2, Fcross0 * cosi );
491 if ( cabs( complexdenominator ) > 1.0e-6 ) {
492 COMPLEX16 complexval = complexnumerator / complexdenominator;
493 detPhaseMag = ( REAL4 )fmin( cabs( complexval ), 10.0 ); //fmin here because sometimes the magnitude value is being divided by a small value and causing biases
494 detPhaseArg = ( REAL4 )gsl_sf_angle_restrict_pos( carg( complexval ) );
495 } else {
496 detPhaseMag = 0.0;
497 detPhaseArg = 0.0;
498 }
499 }
500
501 //When reference detector is off, don't add data that is less sensitive than if the reference detector data would be when it is on
502 if ( detPhaseMag != 0.0 && whichIFOsToBeUsed->data[0] == 0 && detPhaseMag < 1.0 ) {
503 detPhaseMag = 0.0;
504 detPhaseArg = 0.0;
505 }
506
507 //Compute delta values for the frequency difference for the Dirichlet kernel.
508 //This is done either with no assumption on GW frequency: fk*Tsft*[Tdot(jj) - 1] where jj=0,...,N-1
509 //Or it is done with assumption of GW frequency: f*Tsft*[Tdot(jj) - fk/f] where jj=0,...,N-1
510 //Also do shifting and computing floor(delta) and round(delta)
511 //Compute scaling values delta*[delta^2 - 1], and Dirichlet scaling = scaling0/scalingX
512 if ( computeDelta0vals ) {
513 if ( NSparams->assumeNSGWfreq == NULL ) {
514 XLAL_CHECK_NULL( VectorScaleREAL8( delta0vals, GWfrequencyBins, TdotMinus1s->data[0]->data[fftnum], params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
515 } else {
516 XLAL_CHECK_NULL( VectorShiftREAL8( delta0vals, scaledSFTfrequencies, multissb->data[0]->Tdot->data[fftnum], params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
517 XLAL_CHECK_NULL( VectorScaleREAL8( delta0vals, delta0vals, GWfrequencyBins->data[0], params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
518 }
519 memcpy( delta0valsSubset->data, &( delta0vals->data[10] ), sizeof( REAL8 )*delta0valsSubset->length );
520 XLAL_CHECK_NULL( VectorFloorREAL8( floorDelta0valsSubset, delta0valsSubset, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
521 XLAL_CHECK_NULL( VectorRoundREAL8( roundDelta0valsSubset, delta0valsSubset, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
522 XLAL_CHECK_NULL( VectorMultiplyREAL8( DirichletScaling0, delta0valsSubset, delta0valsSubset, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
523 XLAL_CHECK_NULL( VectorShiftREAL8( DirichletScaling0, DirichletScaling0, -1.0, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
524 XLAL_CHECK_NULL( VectorMultiplyREAL8( DirichletScaling0, DirichletScaling0, delta0valsSubset, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
525 }
526 if ( computeDelta0vals ) {
527 computeDelta0vals = 0;
528 }
529 if ( NSparams->assumeNSGWfreq == NULL ) {
530 XLAL_CHECK_NULL( VectorScaleREAL8( deltaXvals, GWfrequencyBins, TdotMinus1s->data[jj]->data[fftnum], params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
531 } else {
532 XLAL_CHECK_NULL( VectorShiftREAL8( deltaXvals, scaledSFTfrequencies, multissb->data[jj]->Tdot->data[fftnum], params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
533 XLAL_CHECK_NULL( VectorScaleREAL8( deltaXvals, deltaXvals, GWfrequencyBins->data[0], params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
534 }
535 memcpy( deltaXvalsSubset->data, &( deltaXvals->data[10] ), sizeof( REAL8 )*deltaXvalsSubset->length );
536 XLAL_CHECK_NULL( VectorRoundREAL8( roundDeltaXvalsSubset, deltaXvalsSubset, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
537 XLAL_CHECK_NULL( VectorSubtractREAL8( diffRoundDeltaValsSubset, roundDelta0valsSubset, roundDeltaXvalsSubset, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
538 XLAL_CHECK_NULL( VectorAddREAL8( deltaXvalsSubset, deltaXvalsSubset, diffRoundDeltaValsSubset, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
539 for ( UINT4 kk = 0; kk < diffRoundDeltaValsSubset->length; kk++ ) {
540 INT4 shiftVal = ( INT4 )diffRoundDeltaValsSubset->data[kk];
541 shiftVector->data[kk] = shiftVal;
542 }
543 XLAL_CHECK_NULL( VectorFloorREAL8( floorDeltaXvalsSubset, deltaXvalsSubset, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
544 XLAL_CHECK_NULL( VectorSubtractREAL8( diffFloorDeltaValsSubset, floorDelta0valsSubset, floorDeltaXvalsSubset, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
545 XLAL_CHECK_NULL( VectorMultiplyREAL8( DirichletScalingX, deltaXvalsSubset, deltaXvalsSubset, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
546 XLAL_CHECK_NULL( VectorShiftREAL8( DirichletScalingX, DirichletScalingX, -1.0, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
547 XLAL_CHECK_NULL( VectorMultiplyREAL8( DirichletScalingX, DirichletScalingX, deltaXvalsSubset, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
548
549 for ( UINT4 kk = 0; kk < scaling->length; kk++ ) {
550 scaling->data[kk] = DirichletScaling0->data[kk] / DirichletScalingX->data[kk];
551 }
552 XLAL_CHECK_NULL( DirichletRatioVector( Dratio, delta0valsSubset, deltaXvalsSubset, scaling, params ) == XLAL_SUCCESS, XLAL_EFUNC );
553
554 //Compute values outside inner loop
555 //If no assumed GW frequency, then normalize Dirichlet kernel ratio values
556 XLAL_CHECK_NULL( VectorCabsfCOMPLEX8( cabsDratio, Dratio, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
557 if ( NSparams->assumeNSGWfreq == NULL ) {
558 for ( UINT4 kk = 0; kk < Dratio->length; kk++ ) if ( cabsDratio->data[kk] != 0.0 ) {
559 Dratio->data[kk] /= cabsDratio->data[kk];
560 }
561 } else {
562 for ( UINT4 kk = 0; kk < Dratio->length; kk++ ) {
563 if ( cabsDratio->data[kk] != 0.0 ) {
564 Dratio->data[kk] /= cabsDratio->data[kk];
565 }
566 }
567 }
568 //compute the complex detector-phase relationship
569 COMPLEX8 detPhaseVal = cpolarf( detPhaseMag, detPhaseArg );
570 //If no assumed GW frequency, compute absDiffFloorDeltaValsSubset and phase shift by pi for when signals are on different "sides" of a bin, to reduce error and improve detection efficiency
571 if ( NSparams->assumeNSGWfreq == NULL ) {
572 XLAL_CHECK_NULL( VectorAbsREAL8( absDiffFloorDeltaValsSubset, diffFloorDeltaValsSubset, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
573 for ( UINT4 kk = 0; kk < absDiffFloorDeltaValsSubset->length; kk++ ) if ( absDiffFloorDeltaValsSubset->data[kk] >= 1.0 ) {
574 TwoPiGWfrequenciesTau->data[kk + 10] -= LAL_PI;
575 }
576 }
577 //Cast TwoPiGWfrequenciesTau into a REAL4VectorAligned and take sin and cos
578 //We want to compute exp[-i*(2*pi*f*tau)] = cos[-2*pi*f*tau]+i*sin[-2*pi*f*tau] = cos[2*pi*f*tau]-i*sin[2*pi*f*tau]
579 //Sometimes this is exp[-i*(2*pi*f*tau)+i*pi] = exp[-i*(2*pi*f*tau-pi)] = cos[2*pi*f*tau-pi]-i*sin[2*pi*f*tau-pi]
580 for ( UINT4 kk = 0; kk < TwoPiGWfrequenciesTau->length; kk++ ) {
581 TwoPiGWfrequenciesTau_f->data[kk] = ( REAL4 )TwoPiGWfrequenciesTau->data[kk];
582 }
583 XLAL_CHECK_NULL( XLALVectorSinCosREAL4( sinTwoPiGWfrequenciesTau->data, cosTwoPiGWfrequenciesTau->data, TwoPiGWfrequenciesTau_f->data, TwoPiGWfrequenciesTau_f->length ) == XLAL_SUCCESS, XLAL_EFUNC );
584
585 //Now finish the computation with the Dirichlet kernel ratio and final correction
586 for ( UINT4 kk = 0; kk < sftcopySubset->data->length; kk++ ) {
587 //The complex coefficient to scale SFT bins
588 COMPLEX8 complexfactor = detPhaseVal * crectf( cosTwoPiGWfrequenciesTau->data[kk + 10], -sinTwoPiGWfrequenciesTau->data[kk + 10] ) * conj( Dratio->data[kk] );
589
590 REAL4 noiseWeighting = 1.0;
591 if ( kk >= numSFTbins2skip && kk < numFbinsInBackground + numSFTbins2skip && createSFT ) {
592 if ( NSparams->assumeNSGWfreq == NULL ) {
593 backgroundScaling->data[fftnum2 * numFbinsInBackground + ( kk - numSFTbins2skip )] = detPhaseMag * detPhaseMag;
594 } else {
595 backgroundScaling->data[fftnum2 * numFbinsInBackground + ( kk - numSFTbins2skip )] = detPhaseMag * detPhaseMag;
596 }
597 } else if ( kk >= numSFTbins2skip && kk < numFbinsInBackground + numSFTbins2skip && !createSFT ) {
598 noiseWeighting = backgroundRatio->data[jj]->data[fftnum];
599 if ( NSparams->assumeNSGWfreq == NULL ) {
600 backgroundScaling->data[fftnum2 * numFbinsInBackground + ( kk - numSFTbins2skip )] += noiseWeighting * detPhaseMag * detPhaseMag;
601 } else {
602 backgroundScaling->data[fftnum2 * numFbinsInBackground + ( kk - numSFTbins2skip )] += noiseWeighting * detPhaseMag * detPhaseMag;
603 }
604 }
605
606 //Apply the corrections and shift to correct bin with the shiftVector
607 sftcopySubset->data->data[kk] = multiSFTvector->data[jj]->data[whichSFTinMultiSFTvector->data[jj]].data->data[kk + 10 - shiftVector->data[kk]] * noiseWeighting * complexfactor;
608 }
609
610 if ( createSFT ) {
611 XLAL_CHECK_NULL( XLALCopySFT( &( combinedSFTs->data[ii] ), sftcopySubset ) == XLAL_SUCCESS, XLAL_EFUNC );
612 createSFT = 0;
613 } else {
614 XLAL_CHECK_NULL( XLALSFTAdd( &( combinedSFTs->data[ii] ), sftcopySubset ) == XLAL_SUCCESS, XLAL_EFUNC );
615 }
616 XLALDestroySFT( sftcopySubset );
617 whichSFTinMultiSFTvector->data[jj]++;
618 }
619 } //loop over detectors
620 } //loop over SFT times
621
622 XLALDestroyINT4Vector( whichSFTinMultiSFTvector );
623 XLALDestroyINT4Vector( whichIFOsToBeUsed );
625 XLALDestroyREAL4VectorAligned( onePlusCosiSqOver2Vals );
627 XLALDestroyREAL4VectorAligned( sin2psiVec );
628 XLALDestroyREAL4VectorAligned( cos2psiVec );
635 XLALDestroyREAL4VectorAligned( TwoPiGWfrequenciesTau_f );
636 XLALDestroyREAL4VectorAligned( sinTwoPiGWfrequenciesTau );
637 XLALDestroyREAL4VectorAligned( cosTwoPiGWfrequenciesTau );
638 destroyAlignedREAL8Vector( GWfrequencies );
639 destroyAlignedREAL8Vector( GWfrequencyBins );
640 destroyAlignedREAL8Vector( SFTfrequencies );
641 if ( NSparams->assumeNSGWfreq != NULL ) {
642 destroyAlignedREAL8Vector( scaledSFTfrequencies );
643 XLALDestroyREAL4TimeSeries( GWfreqInSSB );
644 }
645 destroyAlignedREAL8VectorArray( twoPiTauVals );
646 destroyAlignedREAL8VectorArray( TdotMinus1s );
647 destroyAlignedREAL8Vector( delta0vals );
648 destroyAlignedREAL8Vector( deltaXvals );
649 destroyAlignedREAL8Vector( delta0valsSubset );
650 destroyAlignedREAL8Vector( deltaXvalsSubset );
651 destroyAlignedREAL8Vector( floorDelta0valsSubset );
652 destroyAlignedREAL8Vector( floorDeltaXvalsSubset );
653 destroyAlignedREAL8Vector( roundDelta0valsSubset );
654 destroyAlignedREAL8Vector( roundDeltaXvalsSubset );
655 destroyAlignedREAL8Vector( diffFloorDeltaValsSubset );
656 destroyAlignedREAL8Vector( absDiffFloorDeltaValsSubset );
657 destroyAlignedREAL8Vector( diffRoundDeltaValsSubset );
658 XLALDestroyINT4Vector( shiftVector );
659 destroyAlignedREAL8Vector( TwoPiGWfrequenciesTau );
664 destroyAlignedREAL8Vector( DirichletScaling0 );
665 destroyAlignedREAL8Vector( DirichletScalingX );
668 XLALDestroyREAL4VectorAligned( cabsDratio );
669
670 fprintf( stderr, "done\n" );
671
672 REAL4VectorAligned *tfdata = NULL;
673 XLAL_CHECK_NULL( ( tfdata = convertSFTdataToPowers( combinedSFTs, params, 2.0 / params->Tsft / ( multiNoiseFloor.sqrtSn[0] * multiNoiseFloor.sqrtSn[0] ) ) ) != NULL, XLAL_EFUNC );
674 XLALDestroySFTVector( combinedSFTs );
675
676 return tfdata;
677
678} // coherentlyAddSFTs()
679
680/**
681 * Compute the GW frequency time series (in the SSB) for a set of assumed parameters. This is supposed to be the same as what is done for Makefakedata
682 * \param [in] pulsarParams Pointer to PulsarParams struct
683 * \param [in] tepoch Start time in SSB
684 * \param [in] Tsft SFT length (s)
685 * \param [in] SFToverlap Overlap of the SFTs (s)
686 * \param [in] duration Total duration of the search
687 * \return Pointer to REAL4TimeSeries containing the resulting frequency as a function of time
688 */
690{
691 XLAL_CHECK_NULL( pulsarParams != NULL && duration > 0, XLAL_EINVAL );
692
694 sourceParams.psi = pulsarParams->Amp.psi;
695 sourceParams.aPlus = pulsarParams->Amp.aPlus;
696 sourceParams.aCross = pulsarParams->Amp.aCross;
697 sourceParams.phi0 = pulsarParams->Amp.phi0;
698 sourceParams.f0 = pulsarParams->Doppler.fkdot[0];
699 sourceParams.position.latitude = pulsarParams->Doppler.Delta;
700 sourceParams.position.longitude = pulsarParams->Doppler.Alpha;
701 sourceParams.position.system = COORDINATESYSTEM_EQUATORIAL;
702
703 sourceParams.orbitEpoch = pulsarParams->Doppler.tp;
704 sourceParams.omega = pulsarParams->Doppler.argp;
705 sourceParams.oneMinusEcc = 1.0 - pulsarParams->Doppler.ecc;
706 sourceParams.rPeriNorm = pulsarParams->Doppler.asini * sourceParams.oneMinusEcc;
707 sourceParams.angularSpeed = ( LAL_TWOPI / pulsarParams->Doppler.period ) * sqrt( ( 1.0 + pulsarParams->Doppler.ecc ) / pow( sourceParams.oneMinusEcc, 3.0 ) );
708
709 sourceParams.spinEpoch = pulsarParams->Doppler.refTime;
710
711 sourceParams.deltaT = 5;
712
713 sourceParams.epoch = tepoch;
714 sourceParams.length = ( UINT4 )round( duration / sourceParams.deltaT );
715
716 PulsarCoherentGW XLAL_INIT_DECL( sourceSignal );
717
718 XLAL_CHECK_NULL( XLALGenerateSpinOrbitCW( &sourceSignal, &sourceParams ) == XLAL_SUCCESS, XLAL_EFUNC );
719
720 UINT4 downsampledSeriesLength = ( UINT4 )floor( duration / SFToverlap - 1 );
721
722 REAL4TimeSeries *out = NULL;
723 XLAL_CHECK_NULL( ( out = XLALCreateREAL4TimeSeries( "Time series", &( sourceParams.spinEpoch ), 0, SFToverlap, &lalHertzUnit, downsampledSeriesLength ) ) != NULL, XLAL_EFUNC );
724
725 REAL4VectorAligned *tempVector = NULL;
726 XLAL_CHECK_NULL( ( tempVector = XLALCreateREAL4VectorAligned( ( UINT4 )round( Tsft / 5 ), 32 ) ) != NULL, XLAL_EFUNC );
727 UINT4 stepsize = tempVector->length / 2;
728 XLAL_CHECK_NULL( ( downsampledSeriesLength - 1 )*stepsize + tempVector->length <= sourceSignal.f->data->length, XLAL_EBADLEN ); //Check to make sure we don't go past the end of the array
729
730 //compute the mean value of the frequency during an SFT
731 for ( UINT4 ii = 0; ii < downsampledSeriesLength; ii++ ) {
732 memcpy( tempVector->data, &( sourceSignal.f->data->data[ii * stepsize] ), sizeof( REAL4 )*tempVector->length );
733 out->data->data[ii] = calcMean( tempVector );
734 }
735
736 //Destroy things
737 XLALDestroyREAL4VectorAligned( tempVector );
738 XLALDestroyREAL4VectorSequence( sourceSignal.a->data );
739 XLALFree( sourceSignal.a );
740 XLALDestroyREAL4TimeSeries( sourceSignal.f );
741 XLALDestroyREAL8TimeSeries( sourceSignal.phi );
742
743 return out;
744}
745
746/**
747 * Convert a SFTVector of sfts into powers
748 * \param [in] sfts Pointer to the SFTVector
749 * \param [in] params Pointer to the UserInput_t
750 * \param [in] normalization Normalization value to prevent underflow
751 * \return Pointer to REAL4VectorAligned containing powers
752 */
753REAL4VectorAligned *convertSFTdataToPowers( const SFTVector *sfts, const UserInput_t *params, const REAL8 normalization )
754{
755
756 XLAL_CHECK_NULL( sfts != NULL && params != NULL, XLAL_EINVAL );
757
758 fprintf( LOG, "Converting band to powers... " );
759 fprintf( stderr, "Converting band to powers... " );
760
761 UINT4 numffts = ( UINT4 )floor( params->Tobs / ( params->Tsft - params->SFToverlap ) - 1 );
762 UINT4 sftlength = sfts->data[0].data->length;
763 XLAL_CHECK_NULL( sftlength != 0, XLAL_EINVAL, "SFT has length of 0!\n" );
764
765 //Allocate output
766 REAL4VectorAligned *tfdata = NULL;
767 XLAL_CHECK_NULL( ( tfdata = XLALCreateREAL4VectorAligned( numffts * sftlength, 32 ) ) != NULL, XLAL_EFUNC );
768
769 //Non-existant SFT counter
770 UINT4 nonexistantsft = 0;
771
772 REAL8 starttime = params->t0;
773
774 //Allocate temporary normalized abs[sftcoeff]^2 vector
775 alignedREAL8Vector *normAbsSFTcoeffSq = NULL;
776 XLAL_CHECK_NULL( ( normAbsSFTcoeffSq = createAlignedREAL8Vector( sftlength, 32 ) ) != NULL, XLAL_EFUNC );
777
778 //Load the data into the output vector, roughly normalizing as we go along from the input value
779 //REAL8 sqrtnorm = sqrt(normalization);
780 for ( UINT4 ii = 0; ii < numffts; ii++ ) {
781 if ( ii - nonexistantsft < sfts->length ) {
782 SFTtype *sft = &( sfts->data[ii - nonexistantsft] );
783 if ( sft->epoch.gpsSeconds == ( INT4 )round( ii * ( params->Tsft - params->SFToverlap ) + starttime ) ) {
784 XLAL_CHECK_NULL( VectorCabsCOMPLEX8( normAbsSFTcoeffSq, sft->data, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
785 XLAL_CHECK_NULL( VectorMultiplyREAL8( normAbsSFTcoeffSq, normAbsSFTcoeffSq, normAbsSFTcoeffSq, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
786 XLAL_CHECK_NULL( VectorScaleREAL8( normAbsSFTcoeffSq, normAbsSFTcoeffSq, normalization, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
787 for ( UINT4 jj = 0; jj < sftlength; jj++ ) {
788 tfdata->data[ii * sftlength + jj] = ( REAL4 )normAbsSFTcoeffSq->data[jj];
789 } /* for jj < sftLength */
790 } else {
791 memset( &( tfdata->data[ii * sftlength] ), 0, sizeof( REAL4 )*sftlength );
792 nonexistantsft++; //increment the nonexistantsft counter
793 }
794 } else {
795 memset( &( tfdata->data[ii * sftlength] ), 0, sizeof( REAL4 )*sftlength );
796 nonexistantsft++; //increment the nonexistantsft counter
797 }
798 } /* for ii < numffts */
799
800 destroyAlignedREAL8Vector( normAbsSFTcoeffSq );
801
802 fprintf( LOG, "done\n" );
803 fprintf( stderr, "done\n" );
804
805 fprintf( LOG, "Duty factor = %f\n", 1.0 - ( REAL4 )nonexistantsft / ( REAL4 )numffts );
806 fprintf( stderr, "Duty factor = %f\n", 1.0 - ( REAL4 )nonexistantsft / ( REAL4 )numffts );
807
808 //REAL4 meanTFdata = calcMean(tfdata);
809 //REAL4 stddev = 0.0;
810 //XLAL_CHECK_NULL( calcStddev(&stddev, tfdata) == XLAL_SUCCESS, XLAL_EFUNC );
811 //fprintf(LOG, "TF before weighting, mean subtraction: mean = %g, std. dev. = %g\n", meanTFdata, stddev);
812 //fprintf(stderr, "TF before weighting, mean subtraction: mean = %g, std. dev. = %g\n", meanTFdata, stddev);
813
814 return tfdata;
815
816} // convertSFTdataToPowers()
817
818/**
819 * Read in the data SFTs in one function
820 * \param [in] params Pointer to UserInput_t
821 * \param [in] normalization Normalization value determined from expected noise background
822 * \param [in] minfbin Frequency value of the minimum frequency bin
823 * \param [in] maxfbin Frequency value of the maximum frequency bin
824 * \return REAL4VectorAligned of SFT powers
825 */
826REAL4VectorAligned *readInSFTs( UserInput_t *params, const REAL8 normalization, const REAL8 minfbin, const REAL8 maxfbin )
827{
828
830
831 MultiSFTVector *sftvector = NULL;
832 XLAL_CHECK_NULL( ( sftvector = getMultiSFTVector( params, minfbin, maxfbin ) ) != NULL, XLAL_EFUNC );
833
834 REAL4VectorAligned *tfdata = NULL;
835 XLAL_CHECK_NULL( ( tfdata = convertSFTdataToPowers( sftvector->data[0], params, normalization ) ) != NULL, XLAL_EFUNC );
836
837 //Destroy stuff
838 XLALDestroyMultiSFTVector( sftvector );
839
840 return tfdata;
841
842} // readInSFTs()
843
844/**
845 * Create a list of timestamps from an SFTCatalog
846 * \param [in] catalog Pointer to an SFTCatalog
847 * \return Pointer to a list of GPS timestamps in a MultiLIGOTimeGPSVector
848 */
850{
851
852 XLAL_CHECK_NULL( catalog != NULL, XLAL_EINVAL );
853
854 //Get the MultiSFTCatalogView
855 MultiSFTCatalogView *catalogView = NULL;
856 XLAL_CHECK_NULL( ( catalogView = XLALGetMultiSFTCatalogView( catalog ) ) != NULL, XLAL_EFUNC );
857
860
861 XLALDestroyMultiSFTCatalogView( catalogView );
862
863 return multiTimestamps;
864
865} // getMultiTimeStampsFromSFTCatalog()
866
867/**
868 * Create a list of timestamps from SFTs that might be a subset from those in an SFTCatalog, applying KS/Kuipers test if desired
869 * \param [in] multiSFTvector Pointer to a MultiSFTVector
870 * \param [in] params Pointer to UserInput_t
871 * \return Pointer to a list of GPS timestamps in a MultiLIGOTimeGPSVector
872 */
874{
875
877
879
880 //parse noise list
881 MultiNoiseFloor multiNoiseFloor;
882 XLAL_CHECK_NULL( XLALParseMultiNoiseFloor( &multiNoiseFloor, params->avesqrtSh, params->IFO->length ) == XLAL_SUCCESS, XLAL_EFUNC );
883
884 if ( params->markBadSFTs && !params->signalOnly ) {
885 //REAL8 tfnormval = 2.0/(params->Tsft*(params->avesqrtSh*params->avesqrtSh));
886
887 XLAL_CHECK_NULL( ( multiTimestamps = XLALCalloc( 1, sizeof( *multiTimestamps ) ) ) != NULL, XLAL_ENOMEM );
888 XLAL_CHECK_NULL( ( multiTimestamps->data = XLALCalloc( multiSFTvector->length, sizeof( *multiTimestamps->data ) ) ) != NULL, XLAL_ENOMEM );
889 multiTimestamps->length = multiSFTvector->length;
890
891 for ( UINT4 ii = 0; ii < multiSFTvector->length; ii++ ) {
892 REAL8 tfnormval = 2.0 / ( params->Tsft * ( multiNoiseFloor.sqrtSn[ii] * multiNoiseFloor.sqrtSn[ii] ) );
893
894 REAL4VectorAligned *tfdata = NULL;
895 XLAL_CHECK_NULL( ( tfdata = convertSFTdataToPowers( multiSFTvector->data[ii], params, tfnormval ) ) != NULL, XLAL_EFUNC );
896
897 INT4Vector *removeTheseSFTs = NULL;
898 XLAL_CHECK_NULL( ( removeTheseSFTs = markBadSFTs( tfdata, params ) ) != NULL, XLAL_EFUNC );
899
900 INT4 numberofsfts = 0, sftlength = ( INT4 )multiSFTvector->data[ii]->data[0].data->length;
901 for ( UINT4 jj = 0; jj < removeTheseSFTs->length; jj++ ) if ( removeTheseSFTs->data[jj] == 0 && tfdata->data[jj * sftlength] != 0.0 ) {
902 numberofsfts++;
903 }
904
905 XLAL_CHECK_NULL( ( multiTimestamps->data[ii] = XLALCreateTimestampVector( numberofsfts ) ) != NULL, XLAL_EFUNC );
906
907 INT4 kk = 0;
908 for ( UINT4 jj = 0; jj < removeTheseSFTs->length; jj++ ) {
909 if ( removeTheseSFTs->data[jj] == 0 && tfdata->data[jj * sftlength] != 0.0 ) {
910 XLALGPSSetREAL8( &( multiTimestamps->data[ii]->data[kk] ), params->t0 + jj * ( params->Tsft - params->SFToverlap ) );
912 kk++;
913 }
914 }
915
916 multiTimestamps->data[ii]->deltaT = params->Tsft;
917
918 XLALDestroyINT4Vector( removeTheseSFTs );
920 }
921 } else {
923 }
924
925 return multiTimestamps;
926
927} // getMultiTimeStampsFromSFTs()
928
929/**
930 * Create a list of timestamps from a segment list
931 * \param [in] filenames LALStringVector of filenames
932 * \param [in] t0 Start time of the search
933 * \param [in] Tsft SFT timespan in seconds
934 * \param [in] SFToverlap Overlap of the SFTs in seconds
935 * \param [in] dur Duration of the search in seconds
936 * \return Pointer to a list of GPS timestamps in a MultiLIGOTimeGPSVector
937 */
938MultiLIGOTimeGPSVector *getMultiTimeStampsFromSegmentsFile( const LALStringVector *filenames, const REAL8 t0, const REAL8 Tsft, const REAL8 SFToverlap, const REAL8 dur )
939{
940
941 XLAL_CHECK_NULL( filenames != NULL, XLAL_EINVAL );
942
944 XLAL_CHECK_NULL( ( multiTimestamps = XLALCalloc( 1, sizeof( *multiTimestamps ) ) ) != NULL, XLAL_ENOMEM );
945 multiTimestamps->length = filenames->length;
946 XLAL_CHECK_NULL( ( multiTimestamps->data = XLALCalloc( multiTimestamps->length, sizeof( *( multiTimestamps->data ) ) ) ) != NULL, XLAL_ENOMEM );
947
948 for ( UINT4 ii = 0; ii < filenames->length; ii++ ) {
950 XLAL_CHECK_NULL( ( timestamps = XLALTimestampsFromSegmentFile( filenames->data[ii], Tsft, SFToverlap, 0, 1 ) ) != NULL, XLAL_EFUNC );
951
952 //Shift to find the correct range of start to end SFTs
953 //startSFTindex will be at the correct first SFT
954 //endSFTindex will be one past the last SFT
955 //Thus, we don't need to be inclusive when taking the difference between these two index values
956 UINT4 startSFTindex = 0, endSFTindex = 0;
957 for ( UINT4 jj = 0; jj < timestamps->length; jj++ ) {
958 REAL8 timestamp = XLALGPSGetREAL8( &( timestamps->data[jj] ) );
959 if ( timestamp < t0 ) {
960 startSFTindex++;
961 } else if ( timestamp >= t0 && timestamp + Tsft <= t0 + dur ) {
962 endSFTindex++;
963 } else if ( timestamp + Tsft > t0 + dur ) {
964 break;
965 }
966 }
967
968 XLAL_CHECK_NULL( ( multiTimestamps->data[ii] = XLALCreateTimestampVector( endSFTindex - startSFTindex ) ) != NULL, XLAL_EFUNC );
969 multiTimestamps->data[ii]->deltaT = Tsft;
970 for ( UINT4 jj = 0; jj < multiTimestamps->data[ii]->length; jj++ ) {
971 multiTimestamps->data[ii]->data[jj] = timestamps->data[startSFTindex + jj];
972 }
973
975 }
976
977 return multiTimestamps;
978
979} // getMultiTimeStampsFromSegmentsFile
980
982{
984
985 MultiLIGOTimeGPSVector *squeezedMultiTS = NULL;
986 XLAL_CHECK_NULL( ( squeezedMultiTS = squeezeMultiTimestamps( multiTimestamps ) ) != NULL, XLAL_EFUNC );
987 XLAL_CHECK_NULL( squeezedMultiTS->length > 0, XLAL_EFAILED, "Squeezed multiTimestamps to zero length" );
988
990 if ( squeezedMultiTS->length < 2 ) {
991 XLAL_CHECK_NULL( ( output = XLALCreateTimestampVector( squeezedMultiTS->data[0]->length ) ) != NULL, XLAL_EFUNC );
992 memcpy( output->data, squeezedMultiTS->data[0]->data, sizeof( LIGOTimeGPS )*output->length );
993 output->deltaT = squeezedMultiTS->data[0]->deltaT;
994 XLALDestroyMultiTimestamps( squeezedMultiTS );
995 return output;
996 }
997
998 //First determine the earliest of the timestamps and the last of the timestamps, adding one timestride to the last of the timestamps
999 LIGOTimeGPS earliestTimestamp, lastTimesample, currentTimestamp;
1000 LIGOTimeGPS smallestGPS = LIGOTIMEGPSZERO, largestGPS = LIGOTIMEGPSZERO;
1001 INT4 ifoWithSmallestGPS = -1, ifoWithLargestGPS = -1;
1002 for ( UINT4 ii = 0; ii < squeezedMultiTS->length; ii++ ) {
1003 if ( ifoWithSmallestGPS < 0 ) {
1004 smallestGPS = squeezedMultiTS->data[ii]->data[0];
1005 ifoWithSmallestGPS = ii;
1006 } else if ( XLALGPSCmp( &( squeezedMultiTS->data[ifoWithSmallestGPS]->data[0] ), &( squeezedMultiTS->data[ii]->data[0] ) ) > 0 ) {
1007 smallestGPS = squeezedMultiTS->data[ii]->data[0];
1008 ifoWithSmallestGPS = ii;
1009 }
1010 if ( ifoWithLargestGPS < 0 ) {
1011 largestGPS = squeezedMultiTS->data[ii]->data[squeezedMultiTS->data[ii]->length - 1];
1012 ifoWithLargestGPS = ii;
1013 } else if ( XLALGPSCmp( &( squeezedMultiTS->data[ifoWithLargestGPS]->data[squeezedMultiTS->data[ifoWithLargestGPS]->length - 1] ), &( squeezedMultiTS->data[ii]->data[squeezedMultiTS->data[ii]->length - 1] ) ) < 0 ) {
1014 largestGPS = squeezedMultiTS->data[ii]->data[squeezedMultiTS->data[ii]->length - 1];
1015 ifoWithLargestGPS = ii;
1016 }
1017 }
1018 earliestTimestamp = smallestGPS;
1019 lastTimesample = largestGPS;
1020 XLAL_CHECK_NULL( XLALGPSAdd( &lastTimesample, squeezedMultiTS->data[0]->deltaT ) != NULL, XLAL_EFUNC );
1021 currentTimestamp = earliestTimestamp;
1022
1023 //Now determine the number of unique timestamps for the length of the output vector
1024 UINT4Vector *indexInTimestampVector = NULL;
1025 XLAL_CHECK_NULL( ( indexInTimestampVector = XLALCreateUINT4Vector( squeezedMultiTS->length ) ) != NULL, XLAL_EFUNC );
1026 memset( indexInTimestampVector->data, 0, sizeof( UINT4 )*indexInTimestampVector->length );
1027 UINT4 totalnumTimestamps = 0;
1028 while ( XLALGPSCmp( &currentTimestamp, &lastTimesample ) < 0 ) {
1029 BOOLEAN foundTimestamp = 0;
1030 for ( UINT4 ii = 0; ii < squeezedMultiTS->length; ii++ ) {
1031 if ( indexInTimestampVector->data[ii] < squeezedMultiTS->data[ii]->length && XLALGPSCmp( &currentTimestamp, &( squeezedMultiTS->data[ii]->data[indexInTimestampVector->data[ii]] ) ) == 0 && !foundTimestamp ) {
1032 foundTimestamp = 1;
1033 totalnumTimestamps++;
1034 ( indexInTimestampVector->data[ii] )++;
1035 } else if ( indexInTimestampVector->data[ii] < squeezedMultiTS->data[ii]->length && XLALGPSCmp( &currentTimestamp, &( squeezedMultiTS->data[ii]->data[indexInTimestampVector->data[ii]] ) ) == 0 && foundTimestamp ) {
1036 ( indexInTimestampVector->data[ii] )++;
1037 }
1038 }
1039 XLAL_CHECK_NULL( XLALGPSAdd( &currentTimestamp, 0.5 * squeezedMultiTS->data[0]->deltaT ) != NULL, XLAL_EFUNC );
1040 }
1041
1042 //Reset index vector
1043 memset( indexInTimestampVector->data, 0, sizeof( UINT4 )*indexInTimestampVector->length );
1044
1045 //Populate the output vector
1046 XLAL_CHECK_NULL( ( output = XLALCreateTimestampVector( totalnumTimestamps ) ) != NULL, XLAL_EFUNC );
1047 output->deltaT = squeezedMultiTS->data[0]->deltaT;
1048 for ( UINT4 ii = 0; ii < output->length; ii++ ) {
1049 smallestGPS.gpsSeconds = 0;
1050 smallestGPS.gpsNanoSeconds = 0;
1051 ifoWithSmallestGPS = -1;
1052 for ( UINT4 jj = 0; jj < squeezedMultiTS->length; jj++ ) {
1053 if ( indexInTimestampVector->data[jj] < squeezedMultiTS->data[jj]->length && ifoWithSmallestGPS < 0 ) {
1054 smallestGPS = squeezedMultiTS->data[jj]->data[indexInTimestampVector->data[jj]];
1055 ifoWithSmallestGPS = jj;
1056 } else if ( indexInTimestampVector->data[jj] < squeezedMultiTS->data[jj]->length && XLALGPSCmp( &( squeezedMultiTS->data[ifoWithSmallestGPS]->data[indexInTimestampVector->data[ifoWithSmallestGPS]] ), &( squeezedMultiTS->data[jj]->data[indexInTimestampVector->data[jj]] ) ) > 0 ) {
1057 smallestGPS = squeezedMultiTS->data[jj]->data[indexInTimestampVector->data[jj]];
1058 ifoWithSmallestGPS = jj;
1059 }
1060 }
1061 output->data[ii] = smallestGPS;
1062 for ( UINT4 jj = 0; jj < squeezedMultiTS->length; jj++ ) {
1063 if ( indexInTimestampVector->data[jj] < squeezedMultiTS->data[jj]->length && XLALGPSCmp( &( squeezedMultiTS->data[jj]->data[indexInTimestampVector->data[jj]] ), &smallestGPS ) == 0 ) {
1064 ( indexInTimestampVector->data[jj] )++;
1065 }
1066 }
1067 }
1068
1069 XLALDestroyUINT4Vector( indexInTimestampVector );
1070 XLALDestroyMultiTimestamps( squeezedMultiTS );
1071
1072 return output;
1073}
1074
1076{
1078 UINT4 nonZeroLength = 0;
1079 for ( UINT4 ii = 0; ii < multiTS->length; ii++ ) if ( multiTS->data[ii]->length > 0 ) {
1080 nonZeroLength++;
1081 }
1082 XLAL_CHECK_NULL( nonZeroLength != 0, XLAL_EFAILED );
1083 MultiLIGOTimeGPSVector *ret = NULL;
1084 XLAL_CHECK_NULL( ( ret = XLALCalloc( 1, sizeof( *ret ) ) ) != NULL, XLAL_ENOMEM );
1085 XLAL_CHECK_NULL( ( ret->data = XLALCalloc( nonZeroLength, sizeof( *ret->data ) ) ) != NULL, XLAL_ENOMEM );
1086 ret->length = nonZeroLength;
1087 UINT4 whichTS = 0;
1088 for ( UINT4 ii = 0; ii < nonZeroLength; ii++ ) {
1089 while ( multiTS->data[whichTS]->length == 0 ) {
1090 whichTS++;
1091 }
1092 XLAL_CHECK_NULL( ( ret->data[ii] = XLALCreateTimestampVector( multiTS->data[whichTS]->length ) ) != NULL, XLAL_EFUNC );
1093 memcpy( ret->data[ii]->data, multiTS->data[whichTS]->data, sizeof( LIGOTimeGPS )*ret->data[ii]->length );
1094 ret->data[ii]->deltaT = multiTS->data[whichTS]->deltaT;
1095 whichTS++;
1096 }
1097 return ret;
1098}
1099
1100/**
1101 * \param [in] multiSFTvector Pointer to MultiSFTVector of the SFTs
1102 * \param [in] params Pointer to the user input
1103 * \return Status value
1104 */
1106{
1107 XLAL_CHECK( multiSFTvector != NULL && params != NULL, XLAL_EINVAL );
1108
1109 MultiLIGOTimeGPSVector *GPStimes = NULL;
1110 XLAL_CHECK( ( GPStimes = getMultiTimeStampsFromSFTs( multiSFTvector, params ) ) != NULL, XLAL_EFUNC );
1111
1112 CHARVector *outputfile = NULL;
1113 XLAL_CHECK( ( outputfile = XLALCreateCHARVector( strlen( params->outdirectory ) + 25 ) ) != NULL, XLAL_EFUNC );
1114
1115 for ( UINT4 ii = 0; ii < GPStimes->length; ii++ ) {
1116 memset( outputfile->data, 0, sizeof( CHAR )*outputfile->length );
1117 sprintf( outputfile->data, "%s/%s-%s", params->outdirectory, params->IFO->data[ii], "timestamps.dat" );
1118
1119 FILE *INSFTTIMES = NULL;
1120 XLAL_CHECK( ( INSFTTIMES = fopen( outputfile->data, "w" ) ) != NULL, XLAL_EIO, "Couldn't fopen %s", outputfile->data );
1121
1122 for ( UINT4 jj = 0; jj < GPStimes->data[ii]->length; jj++ ) {
1123 fprintf( INSFTTIMES, "%d 0\n", GPStimes->data[ii]->data[jj].gpsSeconds );
1124 }
1125
1126 fclose( INSFTTIMES );
1127 }
1128
1129 XLALDestroyMultiTimestamps( GPStimes );
1130 XLALDestroyCHARVector( outputfile );
1131
1132 return XLAL_SUCCESS;
1133}
1134
1135/* Critical values of KS test (from Bickel and Doksum). Does not apply directly (mean determined from distribution)
1136alpha=0.01
1137n 10 20 30 40 50 60 80 n>80
1138 .489 .352 .290 .252 .226 .207 .179 1.628/(sqrt(n)+0.12+0.11/sqrt(n))
1139
1140alpha=0.05
1141n 10 20 30 40 50 60 80 n>80
1142 .409 .294 .242 .210 .188 .172 .150 1.358/(sqrt(n)+0.12+0.11/sqrt(n))
1143
1144alpha=0.1 (E.G derived using root finding)
1145n 10 20 30 40 50 60 80 n>80
1146 .369 .265 .218 .189 .170 .155 .135 1.224/(sqrt(n)+0.12+0.11/sqrt(n))
1147
1148alpha=0.2 (E.G derived using root finding)
1149n 10 20 30 40 50 60 80 n>80
1150 1.073/(sqrt(n)+0.12+0.11/sqrt(n))
1151
1152alpha=0.5 (E.G derived using root finding)
1153n 10 20 30 40 50 60 80 n>80
1154 .249 .179 .147 .128 .115 .105 .091 0.828/(sqrt(n)+0.12+0.11/sqrt(n))
1155
1156alpha=0.9 (E.G derived using root finding)
1157n n>80
1158 0.571/(sqrt(n)+0.12+0.11/sqrt(n))
1159
1160Critical values of Kuiper's test using root finding by E.G.
1161alpha=0.05
1162n n>80
1163 1.747/(sqrt(n)+0.155+0.24/sqrt(n))
1164
1165alpha=0.1
1166n n>80
1167 1.620/(sqrt(n)+0.155+0.24/sqrt(n))
1168
1169alpha=0.2
1170n n>80
1171 1.473/(sqrt(n)+0.155+0.24/sqrt(n))
1172*/
1173/**
1174 * Mark the non-Gaussian SFTs using K-S and Kuiper's tests
1175 * \param [in] tfdata Pointer to REAL4VectorAligned of SFT powers
1176 * \param [in] params Pointer to UserInput_t
1177 * \return Pointer to an INT4Vector with marked SFTs to be removed with 1 and keep SFT with 0
1178 */
1180{
1181
1182 XLAL_CHECK_NULL( tfdata != NULL && params != NULL, XLAL_EINVAL );
1183
1184 fprintf( stderr, "Marking bad SFTs... " );
1185
1186 INT4 numffts = ( INT4 )floor( params->Tobs / ( params->Tsft - params->SFToverlap ) - 1 ); //Number of FFTs
1187 //INT4 numfbins = (INT4)(round(params->fspan*params->Tsft+2.0*params->dfmax*params->Tsft)+12+1)+2*params->maxbinshift+params->blksize-1; //Number of frequency bins
1188 INT4 numfbins = ( INT4 )tfdata->length / numffts;
1189
1190 //Allocate output data vector and a single SFT data vector
1191 INT4Vector *output = NULL;
1192 XLAL_CHECK_NULL( ( output = XLALCreateINT4Vector( numffts ) ) != NULL, XLAL_EFUNC );
1193 memset( output->data, 0, sizeof( INT4 )*output->length );
1194 REAL4VectorAligned *tempvect = NULL;
1195 XLAL_CHECK_NULL( ( tempvect = XLALCreateREAL4VectorAligned( numfbins, 32 ) ) != NULL, XLAL_EFUNC );
1196
1197 //Do the KS and Kuiper test on each SFT
1198 REAL8 ksthreshold = 1.358 / ( sqrt( numfbins ) + 0.12 + 0.11 / sqrt( numfbins ) );
1199 //REAL8 ksthreshold = 1.224/(sqrt(numfbins)+0.12+0.11/sqrt(numfbins)); //This is a tighter restriction
1200 //REAL8 ksthreshold = 1.073/(sqrt(numfbins)+0.12+0.11/sqrt(numfbins)); //This is an even tighter restriction
1201 REAL8 kuiperthreshold = 1.747 / ( sqrt( numfbins ) + 0.155 + 0.24 / sqrt( numfbins ) );
1202 //REAL8 kuiperthreshold = 1.620/(sqrt(numfbins)+0.155+0.24/sqrt(numfbins)); //This is a tighter restriction
1203 //REAL8 kuiperthreshold = 1.473/(sqrt(numfbins)+0.155+0.24/sqrt(numfbins)); //This is an even tighter restriction
1204 INT4 badsfts = 0, totalsfts = 0;
1205 REAL8 kstest = 0.0, kuipertest = 0.0;
1206 for ( INT4 ii = 0; ii < numffts; ii++ ) {
1207 if ( tfdata->data[ii * numfbins] != 0.0 ) {
1208 totalsfts++;
1209 memcpy( tempvect->data, &( tfdata->data[ii * numfbins] ), sizeof( REAL4 )*tempvect->length );
1210 XLAL_CHECK_NULL( ks_test_exp( &kstest, tempvect ) == XLAL_SUCCESS, XLAL_EFUNC );
1211 XLAL_CHECK_NULL( kuipers_test_exp( &kuipertest, tempvect ) == XLAL_SUCCESS, XLAL_EFUNC );
1212 if ( kstest > ksthreshold || kuipertest > kuiperthreshold ) {
1213 output->data[ii] = 1;
1214 badsfts++;
1215 }
1216 }
1217 }
1218
1219 //Destroy stuff
1221
1222 fprintf( stderr, "done\n" );
1223 fprintf( stderr, "Fraction excluded in K-S and Kuiper's tests = %f\n", ( REAL4 )badsfts / ( REAL4 )totalsfts );
1224
1225 return output;
1226
1227} // markBadSFTs()
1228
1229/**
1230 * Remove the marked SFTs as bad by setting values to 0
1231 * \param [in,out] tfdata Pointer to REAL4VectorAligned of SFT powers
1232 * \param [in] badsfts Poienter to INT4Vector of bad SFTs
1233 */
1234void removeBadSFTs( REAL4VectorAligned *tfdata, const INT4Vector *badsfts )
1235{
1236
1237 XLAL_CHECK_VOID( tfdata != NULL && badsfts != NULL, XLAL_EINVAL );
1238 fprintf( stderr, "Removing bad SFTs... " );
1239 UINT4 numfbins_tfdata = tfdata->length / badsfts->length;
1240 for ( UINT4 ii = 0; ii < badsfts->length; ii++ ) if ( badsfts->data[ii] == 1 ) {
1241 memset( &( tfdata->data[ii * numfbins_tfdata] ), 0, sizeof( REAL4 )*numfbins_tfdata );
1242 }
1243 fprintf( stderr, "done.\n" );
1244
1245} // removeBadSFTs()
1246
1247MultiSFTVector *generateSFTdata( UserInput_t *uvar, const MultiLALDetector *detectors, const EphemerisData *edat, const INT4 maxbinshift, const gsl_rng *rng )
1248{
1249
1250 XLAL_CHECK_NULL( uvar != NULL && detectors != NULL && edat != NULL && rng != NULL, XLAL_EINVAL );
1251
1252 MultiSFTVector *multiSFTvector = NULL;
1253
1254 //parse noise list
1255 MultiNoiseFloor multiNoiseFloor;
1256 XLAL_CHECK_NULL( XLALParseMultiNoiseFloor( &multiNoiseFloor, uvar->avesqrtSh, uvar->IFO->length ) == XLAL_SUCCESS, XLAL_EFUNC );
1257
1258 //Determine band size to get the SFT data (remember to get extra bins because of the running median and the bin shifts due to detector velocity) with nudge of 0.1/Tsft for rounding issues
1259 REAL8 minfbin = round( uvar->fmin * uvar->Tsft - uvar->dfmax * uvar->Tsft - 0.5 * ( uvar->blksize - 1 ) - ( REAL8 )( maxbinshift ) - 6.0 ) / uvar->Tsft + 0.1 / uvar->Tsft;
1260 REAL8 maxfbin = round( ( uvar->fmin + uvar->fspan ) * uvar->Tsft + uvar->dfmax * uvar->Tsft + 0.5 * ( uvar->blksize - 1 ) + ( REAL8 )( maxbinshift ) + 6.0 ) / uvar->Tsft - 0.1 / uvar->Tsft;
1261
1262 //Start by getting timestamps or creating them
1264 if ( XLALUserVarWasSet( &uvar->inputSFTs ) ) {
1265 SFTCatalog *catalog = NULL;
1266 XLAL_CHECK_NULL( ( catalog = findSFTdata( uvar ) ) != NULL, XLAL_EFUNC );
1267 MultiSFTVector *tmpsftvector = NULL;
1268 XLAL_CHECK_NULL( ( tmpsftvector = extractSFTband( catalog, minfbin, maxfbin ) ) != NULL, XLAL_EFUNC );
1269 XLAL_CHECK_NULL( ( multiTimestamps = getMultiTimeStampsFromSFTs( tmpsftvector, uvar ) ) != NULL, XLAL_EFUNC );
1270 XLALDestroyMultiSFTVector( tmpsftvector );
1271 XLALDestroySFTCatalog( catalog );
1272 } else if ( XLALUserVarWasSet( &uvar->timestampsFile ) ) {
1273 XLAL_CHECK_NULL( ( multiTimestamps = XLALReadMultiTimestampsFiles( uvar->timestampsFile ) ) != NULL, XLAL_EFUNC );
1274 for ( UINT4 ii = 0; ii < multiTimestamps->length; ii++ ) {
1275 multiTimestamps->data[ii]->deltaT = uvar->Tsft;
1276 }
1277 } else if ( XLALUserVarWasSet( &uvar->segmentFile ) ) {
1278 XLAL_CHECK_NULL( ( multiTimestamps = getMultiTimeStampsFromSegmentsFile( uvar->segmentFile, uvar->t0, uvar->Tsft, uvar->SFToverlap, uvar->Tobs ) ) != NULL, XLAL_EFUNC );
1279 } else {
1280 LIGOTimeGPS tStart;
1281 XLALGPSSetREAL8( &tStart, uvar->t0 );
1282 XLAL_CHECK_NULL( xlalErrno == 0, XLAL_EFUNC, "XLALGPSSetREAL8 failed\n" );
1283 XLAL_CHECK_NULL( ( multiTimestamps = XLALMakeMultiTimestamps( tStart, uvar->Tobs, uvar->Tsft, uvar->SFToverlap, detectors->length ) ) != NULL, XLAL_EFUNC );
1284 }
1285
1286 //TwoSpect to analyze:
1287 REAL8 TwoSpectFmin = round( uvar->fmin * uvar->Tsft - uvar->dfmax * uvar->Tsft - 0.5 * ( uvar->blksize - 1 ) - ( REAL8 )maxbinshift - 6.0 ) / uvar->Tsft;
1288 REAL8 TwoSpectBand = round( uvar->fspan * uvar->Tsft + 2.0 * uvar->dfmax * uvar->Tsft + ( uvar->blksize - 1 ) + ( REAL8 )( 2.0 * maxbinshift ) + 12.0 ) / uvar->Tsft;
1289
1290 //Setup the MFD data parameters
1291 CWMFDataParams XLAL_INIT_DECL( DataParams );
1292 if ( XLALUserVarWasSet( &uvar->injFmin ) && XLALUserVarWasSet( &uvar->injBand ) && uvar->injFmin <= TwoSpectFmin && uvar->injFmin + uvar->injBand >= TwoSpectFmin + TwoSpectBand ) {
1293 DataParams.fMin = uvar->injFmin;
1294 DataParams.Band = uvar->injBand;
1295 } else {
1296 DataParams.fMin = TwoSpectFmin;
1297 DataParams.Band = TwoSpectBand;
1298 }
1299 DataParams.multiIFO.length = detectors->length;
1300 for ( UINT4 ii = 0; ii < detectors->length; ii++ ) {
1301 DataParams.multiIFO.sites[ii] = detectors->sites[ii];
1302 }
1303 DataParams.multiNoiseFloor.length = detectors->length;
1304 for ( UINT4 ii = 0; ii < detectors->length; ii++ ) {
1305 DataParams.multiNoiseFloor.sqrtSn[ii] = 0.0;
1306 }
1307 DataParams.multiTimestamps = multiTimestamps;
1308 DataParams.randSeed = uvar->injRandSeed;
1309 DataParams.SFTWindowType = "Hann";
1310 DataParams.SFTWindowParam = 0;
1311
1312 MultiSFTVector *signalSFTs = NULL;
1314 //If injection sources then read them and make signal sfts
1315 if ( XLALUserVarWasSet( &uvar->injectionSources ) ) {
1316 XLAL_CHECK_NULL( ( injectionSources = XLALPulsarParamsFromUserInput( uvar->injectionSources, NULL ) ) != NULL, XLAL_EFUNC );
1317
1318 fprintf( stderr, "Generating signal SFTs with %d signals... ", ( INT4 )injectionSources->length );
1319
1320 if ( !uvar->signalOnly ) {
1321 XLAL_CHECK_NULL( XLALCWMakeFakeMultiData( &signalSFTs, NULL, injectionSources, &DataParams, edat ) == XLAL_SUCCESS, XLAL_EFUNC );
1322 } else {
1323 XLAL_CHECK_NULL( XLALCWMakeFakeMultiData( &multiSFTvector, NULL, injectionSources, &DataParams, edat ) == XLAL_SUCCESS, XLAL_EFUNC );
1324 }
1325
1326 if ( DataParams.fMin != TwoSpectFmin ) {
1327 XLAL_CHECK_NULL( XLALMultiSFTVectorResizeBand( signalSFTs, TwoSpectFmin, TwoSpectBand ) == XLAL_SUCCESS, XLAL_EFUNC );
1328 }
1329
1330 fprintf( stderr, "done\n" );
1331 } // if there are injections
1332
1333 //If not signal only, create sfts that include noise or extract a band from real data
1334 if ( !uvar->signalOnly ) {
1335 if ( uvar->gaussNoiseWithSFTgaps || XLALUserVarWasSet( &uvar->timestampsFile ) || XLALUserVarWasSet( &uvar->segmentFile ) || !XLALUserVarWasSet( &uvar->inputSFTs ) ) {
1336 for ( UINT4 ii = 0; ii < detectors->length; ii++ ) {
1337 DataParams.multiNoiseFloor.sqrtSn[ii] = multiNoiseFloor.sqrtSn[ii];
1338 }
1339
1340 fprintf( stderr, "Generating noise SFTs... " );
1341
1342 XLAL_CHECK_NULL( XLALCWMakeFakeMultiData( &multiSFTvector, NULL, NULL, &DataParams, edat ) == XLAL_SUCCESS, XLAL_EFUNC );
1343
1344 if ( DataParams.fMin != TwoSpectFmin ) {
1345 XLAL_CHECK_NULL( XLALMultiSFTVectorResizeBand( multiSFTvector, TwoSpectFmin, TwoSpectBand ) == XLAL_SUCCESS, XLAL_EFUNC );
1346 }
1347
1348 fprintf( stderr, "done\n" );
1349 } else {
1350 SFTCatalog *catalog = NULL;
1351 XLAL_CHECK_NULL( ( catalog = findSFTdata( uvar ) ) != NULL, XLAL_EFUNC );
1352 MultiSFTVector *tmpmultiSFTvector = NULL;
1353 XLAL_CHECK_NULL( ( tmpmultiSFTvector = extractSFTband( catalog, minfbin, maxfbin ) ) != NULL, XLAL_EFUNC );
1354 XLAL_CHECK_NULL( ( multiSFTvector = XLALExtractMultiSFTVectorWithMultiTimestamps( tmpmultiSFTvector, multiTimestamps ) ) != NULL, XLAL_EFUNC );
1355 XLALDestroyMultiSFTVector( tmpmultiSFTvector );
1356 XLALDestroySFTCatalog( catalog );
1357 }
1358 } // if not signal only SFTs
1359
1360 //Add the SFT vectors together
1361 if ( XLALUserVarWasSet( &uvar->injectionSources ) && !uvar->signalOnly ) {
1362 XLAL_CHECK_NULL( XLALMultiSFTVectorAdd( multiSFTvector, signalSFTs ) == XLAL_SUCCESS, XLAL_EFUNC );
1363 XLALDestroyMultiSFTVector( signalSFTs );
1364 }
1365
1366 //If printing the data outputs, then do that here
1367 if ( ( XLALUserVarWasSet( &uvar->printSignalData ) || XLALUserVarWasSet( &uvar->printMarginalizedSignalData ) ) && XLALUserVarWasSet( &uvar->injectionSources ) ) {
1368 for ( UINT4 ii = 0; ii < detectors->length; ii++ ) {
1369 DataParams.multiNoiseFloor.sqrtSn[ii] = 0.0;
1370 }
1371 PulsarParamsVector *oneSignal = NULL;
1372 XLAL_CHECK_NULL( ( oneSignal = XLALCreatePulsarParamsVector( 1 ) ) != NULL, XLAL_EFUNC );
1373
1374 FILE *SIGNALOUT = NULL, *MARGINALIZEDSIGNALOUT = NULL;
1375 if ( XLALUserVarWasSet( &uvar->printSignalData ) ) {
1376 XLAL_CHECK_NULL( ( SIGNALOUT = fopen( uvar->printSignalData, "w" ) ) != NULL, XLAL_EIO, "Failed to open %s for writing\n", uvar->printSignalData );
1377 }
1378 if ( XLALUserVarWasSet( &uvar->printMarginalizedSignalData ) ) {
1379 XLAL_CHECK_NULL( ( MARGINALIZEDSIGNALOUT = fopen( uvar->printMarginalizedSignalData, "w" ) ) != NULL, XLAL_EIO, "Failed to open %s for writing\n", uvar->printMarginalizedSignalData );
1380 }
1381
1382 for ( UINT4 ii = 0; ii < injectionSources->length; ii++ ) {
1383 memcpy( oneSignal->data, &( injectionSources->data[ii] ), sizeof( injectionSources->data[0] ) );
1384 if ( XLALUserVarWasSet( &uvar->printSignalData ) ) {
1385 MultiSFTVector *oneSignalSFTs = NULL;
1386 XLAL_CHECK_NULL( XLALCWMakeFakeMultiData( &oneSignalSFTs, NULL, oneSignal, &DataParams, edat ) == XLAL_SUCCESS, XLAL_EFUNC );
1387 if ( DataParams.fMin != TwoSpectFmin ) {
1388 XLAL_CHECK_NULL( XLALMultiSFTVectorResizeBand( oneSignalSFTs, TwoSpectFmin, TwoSpectBand ) == XLAL_SUCCESS, XLAL_EFUNC );
1389 }
1390
1391 alignedREAL8Vector *aveSFTsPower = NULL;
1392 XLAL_CHECK_NULL( ( aveSFTsPower = createAlignedREAL8Vector( multiSFTvector->data[0]->data->data->length, 32 ) ) != NULL, XLAL_EFUNC );
1393 memset( aveSFTsPower->data, 0, sizeof( REAL8 )*aveSFTsPower->length );
1394
1395 alignedREAL8Vector *SFTpower = NULL;
1396 XLAL_CHECK_NULL( ( SFTpower = createAlignedREAL8Vector( aveSFTsPower->length, 32 ) ) != NULL, XLAL_EFUNC );
1397
1398 for ( UINT4 jj = 0; jj < oneSignalSFTs->data[0]->length; jj++ ) {
1399 SFTtype *sft = &( oneSignalSFTs->data[0]->data[jj] );
1400 XLAL_CHECK_NULL( VectorCabsCOMPLEX8( SFTpower, sft->data, uvar->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
1401 XLAL_CHECK_NULL( VectorMultiplyREAL8( SFTpower, SFTpower, SFTpower, uvar->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
1402 XLAL_CHECK_NULL( VectorScaleREAL8( SFTpower, SFTpower, 2.0 / uvar->Tsft, uvar->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
1403 //for (UINT4 kk=0; kk<aveSFTsPower->length; kk++) {
1404 // REAL8 powerval = 2.0*(creal(sft->data->data[kk])*creal(sft->data->data[kk]) + cimag(sft->data->data[kk])*cimag(sft->data->data[kk]))/uvar->Tsft;
1405 // SFTpower->data[kk] = powerval;
1406 // aveSFTsPower->data[kk] += powerval;
1407 //}
1408 XLAL_CHECK_NULL( VectorAddREAL8( aveSFTsPower, aveSFTsPower, SFTpower, uvar->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
1409 //UINT4 max_index = max_index_double((REAL8Vector*)SFTpower);
1410 //fprintf(SIGNALOUT, "%d %.6g\n", max_index, SFTpower->data[max_index]);
1411 }
1412 for ( UINT4 jj = 0; jj < aveSFTsPower->length; jj++ ) {
1413 fprintf( SIGNALOUT, "%.9g %.9g\n", DataParams.fMin + jj / uvar->Tsft, aveSFTsPower->data[jj] / multiSFTvector->data[0]->length );
1414 }
1415 XLALDestroyMultiSFTVector( oneSignalSFTs );
1416 destroyAlignedREAL8Vector( aveSFTsPower );
1417 destroyAlignedREAL8Vector( SFTpower );
1418 }
1419 if ( XLALUserVarWasSet( &uvar->printMarginalizedSignalData ) ) {
1420 REAL8Vector *marginalizedSignalData = NULL;
1421 XLAL_CHECK_NULL( ( marginalizedSignalData = XLALCreateREAL8Vector( multiSFTvector->data[0]->data->data->length ) ) != NULL, XLAL_EFUNC );
1422 memset( marginalizedSignalData->data, 0, sizeof( REAL8 )*marginalizedSignalData->length );
1423 for ( UINT4 jj = 0; jj < 300; jj++ ) {
1424 REAL8 cosi = 2.0 * gsl_rng_uniform( rng ) - 1.0;
1425 REAL8 h0 = 1.0;
1426 oneSignal->data[0].Amp.aPlus = 0.5 * h0 * ( 1.0 + cosi * cosi );
1427 oneSignal->data[0].Amp.aCross = h0 * cosi;
1428 oneSignal->data[0].Amp.psi = LAL_TWOPI * gsl_rng_uniform( rng );
1429 oneSignal->data[0].Amp.phi0 = LAL_TWOPI * gsl_rng_uniform( rng );
1430 oneSignal->data[0].Doppler.argp = LAL_TWOPI * gsl_rng_uniform( rng );
1431 MultiSFTVector *oneSignalSFTs = NULL;
1432 XLAL_CHECK_NULL( XLALCWMakeFakeMultiData( &oneSignalSFTs, NULL, oneSignal, &DataParams, edat ) == XLAL_SUCCESS, XLAL_EFUNC );
1433 if ( DataParams.fMin != TwoSpectFmin ) {
1434 XLAL_CHECK_NULL( XLALMultiSFTVectorResizeBand( oneSignalSFTs, TwoSpectFmin, TwoSpectBand ) == XLAL_SUCCESS, XLAL_EFUNC );
1435 }
1436
1437 for ( UINT4 kk = 0; kk < oneSignalSFTs->data[0]->length; kk++ ) {
1438 SFTtype *sft = &( oneSignalSFTs->data[0]->data[kk] );
1439 for ( UINT4 ll = 0; ll < marginalizedSignalData->length; ll++ ) {
1440 marginalizedSignalData->data[ll] += ( 2.0 * ( creal( sft->data->data[ll] ) * creal( sft->data->data[ll] ) + cimag( sft->data->data[ll] ) * cimag( sft->data->data[ll] ) ) / uvar->Tsft );
1441 }
1442 }
1443 XLALDestroyMultiSFTVector( oneSignalSFTs );
1444 } //Loop over trials
1445 for ( UINT4 jj = 0; jj < marginalizedSignalData->length; jj++ ) {
1446 marginalizedSignalData->data[jj] /= 300.0 * multiSFTvector->data[0]->length;
1447 fprintf( MARGINALIZEDSIGNALOUT, "%.9g %.9g\n", DataParams.fMin + jj / uvar->Tsft, marginalizedSignalData->data[jj] );
1448 }
1449 XLALDestroyREAL8Vector( marginalizedSignalData );
1450 } //If printing marginalized data
1451 } //loop over the number of injected sources
1452 memset( oneSignal->data, 0, sizeof( injectionSources->data[0] ) );
1453 XLALDestroyPulsarParamsVector( oneSignal );
1454 if ( XLALUserVarWasSet( &uvar->printSignalData ) ) {
1455 fclose( SIGNALOUT );
1456 }
1457 if ( XLALUserVarWasSet( &uvar->printMarginalizedSignalData ) ) {
1458 fclose( MARGINALIZEDSIGNALOUT );
1459 }
1460 } //end printing data
1461
1462 if ( XLALUserVarWasSet( &uvar->injectionSources ) ) {
1464 }
1466
1467 return multiSFTvector;
1468
1469} // generateSFTdata()
1470
1471/**
1472 * Slide the time-frequency data to account for detector motion
1473 * \param [out] output Pointer to REAL4VectorAligned of SFT powers that have been corrected
1474 * \param [in] params Pointer to UserInput_t
1475 * \param [in] tfdata Pointer to REAL4VectorAligned of SFT powers
1476 * \param [in] binshifts Pointer to INT4Vector of bin shift values
1477 * \return Status value
1478 */
1480{
1481
1482 XLAL_CHECK( output != NULL && params != NULL && tfdata != NULL && binshifts != NULL, XLAL_EINVAL );
1483
1484 fprintf( stderr, "Sliding TF data... " );
1485
1486 UINT4 numffts = ( UINT4 )floor( params->Tobs / ( params->Tsft - params->SFToverlap ) - 1 );
1487 UINT4 numfbins = output->length / numffts;
1488 UINT4 maxbinshift = ( tfdata->length / numffts - numfbins ) / 2;
1489
1490 for ( UINT4 ii = 0; ii < numffts; ii++ ) {
1491 XLAL_CHECK( abs( binshifts->data[ii] ) < ( INT4 )maxbinshift, XLAL_EFAILED, "SFT slide value %d is greater than maximum value predicted (%d)", binshifts->data[ii], maxbinshift );
1492 memcpy( &( output->data[ii * numfbins] ), &( tfdata->data[ii * ( numfbins + 2 * maxbinshift ) + maxbinshift + binshifts->data[ii]] ), sizeof( REAL4 )*numfbins );
1493 }
1494
1495 //fprintf(stderr, "Mean = %g ", calcMean(output));
1496
1497 fprintf( stderr, "done\n" );
1498
1499 return XLAL_SUCCESS;
1500
1501} // slideTFdata()
1502
1503/**
1504 * Determine the running mean of each SFT
1505 * \param [out] output Pointer to REAL4VectorAligned of running mean values of each SFT
1506 * \param [in] tfdata Pointer to REAL4VectorAligned of SFT powers
1507 * \param [in] numffts Number of SFTs in the observation time
1508 * \param [in] numfbins Number of frequency bins
1509 * \param [in] blksize Number of bins in the running median
1510 * \return Status value
1511 */
1512INT4 tfRngMeans( REAL4VectorAligned *output, const REAL4VectorAligned *tfdata, const UINT4 numffts, const UINT4 numfbins, const UINT4 blksize )
1513{
1514
1515 XLAL_CHECK( output != NULL && tfdata != NULL && numffts > 0 && numfbins > 0 && blksize > 0, XLAL_EINVAL );
1516
1517 fprintf( LOG, "Assessing SFT background... " );
1518 fprintf( stderr, "Assessing SFT background... " );
1519
1521 REAL8 bias;
1522 UINT4 totalfbins = numfbins + blksize - 1;
1523
1524 //Blocksize of running median
1525 LALRunningMedianPar block = {blksize};
1526
1527 //Running median bias calculation
1528 if ( blksize < 1000 ) {
1529 bias = XLALRngMedBias( blksize );
1531 } else {
1532 bias = LAL_LN2;
1533 }
1534 //REAL8 invbias = 1.0/(bias*1.0099993480677538); //StackSlide normalization for 101 bins
1535 REAL8 invbias = 1.0 / bias;
1536
1537 //Allocate for a single SFT data and the medians out of each SFT
1538 REAL4VectorAligned *inpsd = NULL, *mediansout = NULL;
1539 XLAL_CHECK( ( inpsd = XLALCreateREAL4VectorAligned( totalfbins, 32 ) ) != NULL, XLAL_EFUNC );
1540 XLAL_CHECK( ( mediansout = XLALCreateREAL4VectorAligned( numfbins, 32 ) ) != NULL, XLAL_EFUNC );
1541
1542 //Now do the running median
1543 for ( UINT4 ii = 0; ii < numffts; ii++ ) {
1544 //If the SFT values were not zero, then compute the running median
1545 if ( tfdata->data[ii * totalfbins] != 0.0 ) {
1546 //Determine running median value, convert to mean value
1547 memcpy( inpsd->data, &( tfdata->data[ii * inpsd->length] ), sizeof( REAL4 )*inpsd->length );
1548
1549 //calculate running median
1550 LALSRunningMedian2( &status, ( REAL4Vector * )mediansout, ( REAL4Vector * )inpsd, block );
1551 XLAL_CHECK( status.statusCode == 0, XLAL_EFUNC );
1552
1553 //Now make the output medians into means by multiplying by 1/bias
1554 for ( UINT4 jj = 0; jj < mediansout->length; jj++ ) {
1555 output->data[ii * numfbins + jj] = ( REAL4 )( mediansout->data[jj] * invbias );
1556 }
1557 } else {
1558 //Otherwise, set means to zero
1559 //for (UINT4 jj=0; jj<mediansout->length; jj++) output->data[ii*numfbins + jj] = 0.0;
1560 memset( &( output->data[ii * numfbins] ), 0, sizeof( REAL4 )*mediansout->length );
1561 }
1562 } /* for ii < numffts */
1563
1564 //Destroy stuff
1566 XLALDestroyREAL4VectorAligned( mediansout );
1567
1568 fprintf( LOG, "done\n" );
1569 fprintf( stderr, "done\n" );
1570 fprintf( stderr, "Mean of running means = %g\n", calcMean( output ) );
1571
1572 return XLAL_SUCCESS;
1573
1574} // tfRngMeans()
1575
1576/**
1577 * \param [in,out] tfdataarray Pointer to a REAL4VectorAlignedArray that has the time-frequency data of powers
1578 * \param [in] numffts Number of FFTs in the total observation time
1579 * \return Status value
1580 */
1582{
1583 XLAL_CHECK( tfdataarray != NULL && numffts > 0, XLAL_EINVAL );
1584 if ( tfdataarray->length == 1 ) {
1585 return XLAL_SUCCESS;
1586 }
1587 UINT4 sftlength = tfdataarray->data[0]->length / numffts;
1588 for ( UINT4 ii = 0; ii < numffts; ii++ ) {
1589 if ( tfdataarray->data[0]->data[ii * sftlength] == 0.0 ) {
1590 for ( UINT4 jj = 1; jj < tfdataarray->length; jj++ ) {
1591 if ( tfdataarray->data[jj]->data[ii * sftlength] != 0.0 ) {
1592 memcpy( &( tfdataarray->data[0]->data[ii * sftlength] ), &( tfdataarray->data[jj]->data[ii * sftlength] ), sizeof( REAL4 )*sftlength );
1593 break;
1594 }
1595 }
1596 }
1597 }
1598 return XLAL_SUCCESS;
1599}
1600
1601/**
1602 * Subtract running mean values from the SFT data, modifying input time-frequency data
1603 * \param [in,out] tfdata Pointer to REAL4VectorAligned time-frequency data (modified by this function!)
1604 * \param [in] rngMeans Pointer to REAL4VectorAligned of running mean values
1605 * \param [in] backgroundScaling Pointer to REAL4VectorAligned of background scaling values
1606 * \param [in] numffts Number of SFTs from observation time
1607 * \param [in] numfbins Number of SFT frequency bins
1608 * \return Status value
1609 */
1610INT4 tfMeanSubtract( REAL4VectorAligned *tfdata, const REAL4VectorAligned *rngMeans, const REAL4VectorAligned *backgroundScaling, const UINT4 numffts, const UINT4 numfbins )
1611{
1612
1613 XLAL_CHECK( tfdata != NULL && rngMeans != NULL && backgroundScaling != NULL && numffts > 0 && numfbins > 0, XLAL_EINVAL );
1614 fprintf( stderr, "Subtracting expected background... " );
1615 for ( UINT4 ii = 0; ii < numffts; ii++ ) if ( rngMeans->data[ii * numfbins] != 0.0 ) for ( UINT4 jj = 0; jj < numfbins; jj++ ) {
1616 tfdata->data[ii * numfbins + jj] -= rngMeans->data[ii * numfbins + jj] * backgroundScaling->data[ii * numfbins + jj];
1617 }
1618 //fprintf(stderr, "TF mean after subtraction = %g ", calcMean(tfdata));
1619 fprintf( stderr, "done\n" );
1620 return XLAL_SUCCESS;
1621
1622} // tfMeanSubtract()
1623
1624/**
1625 * Weight the SFTs based on antenna pattern and noise variance (Equation 11, assuming the input time-frequency data is already mean subtracted)
1626 * \param [out] output Pointer to REAL4VectorAligned of mean subtracted, noise and antenna pattern weighted SFTs
1627 * \param [in] tfdata Pointer to REAL4VectorAligned of mean subtracted SFTs
1628 * \param [in] rngMeans Pointer to REAL4VectorAligned of running mean values
1629 * \param [in] antPatternWeights Pointer to REAL4VectorAligned of antenna pattern weights
1630 * \param [in] backgroundScaling Pointer to REAL4VectorAligned of background scaling values
1631 * \param [in] indexValuesOfExistingSFTs Pointer to INT4Vector of the index values of the existing SFTs
1632 * \param [in] params Pointer to UserInput_t
1633 * \return Status value
1634 */
1635INT4 tfWeight( REAL4VectorAligned *output, const REAL4VectorAligned *tfdata, REAL4VectorAligned *rngMeans, REAL4VectorAligned *antPatternWeights, const REAL4VectorAligned *backgroundScaling, const INT4Vector *indexValuesOfExistingSFTs, const UserInput_t *params )
1636{
1637
1638 XLAL_CHECK( output != NULL && tfdata != NULL && rngMeans != NULL && antPatternWeights != NULL && backgroundScaling != NULL && indexValuesOfExistingSFTs != NULL && params != NULL, XLAL_EINVAL );
1639
1640 fprintf( stderr, "Applying weighting to SFT data... " );
1641
1642 UINT4 numffts = antPatternWeights->length;
1643 UINT4 numfbins = tfdata->length / numffts;
1644
1645 //Initially set output to zero
1646 memset( output->data, 0, sizeof( REAL4 )*output->length );
1647
1648 REAL4VectorAligned *antweightssq0 = NULL, *antweightssq = NULL, *rngMeanssq = NULL, *backgroundScalingSq = NULL;
1649 XLAL_CHECK( ( antweightssq0 = XLALCreateREAL4VectorAligned( numffts, 32 ) ) != NULL, XLAL_EFUNC );
1650 XLAL_CHECK( ( antweightssq = XLALCreateREAL4VectorAligned( numffts, 32 ) ) != NULL, XLAL_EFUNC );
1651 XLAL_CHECK( ( rngMeanssq = XLALCreateREAL4VectorAligned( numffts, 32 ) ) != NULL, XLAL_EFUNC );
1652 XLAL_CHECK( ( backgroundScalingSq = XLALCreateREAL4VectorAligned( numffts, 32 ) ) != NULL, XLAL_EFUNC );
1653
1654 //If user specified that SFTs contain signal only, then the values of the backgrnd vector will be zeros.
1655 //We must set them equal to 1.0 here, then return them to zeros at the end
1656 if ( params->signalOnly ) {
1657 for ( UINT4 ii = 0; ii < indexValuesOfExistingSFTs->length; ii++ ) {
1658 for ( UINT4 jj = 0; jj < numfbins; jj++ ) {
1659 rngMeans->data[numfbins * indexValuesOfExistingSFTs->data[ii] + jj] = 1.0;
1660 }
1661 }
1662 }
1663
1664 //Determine antenna pattern weights squared
1665 XLAL_CHECK( XLALVectorMultiplyREAL4( antweightssq0->data, antPatternWeights->data, antPatternWeights->data, antPatternWeights->length ) == XLAL_SUCCESS, XLAL_EFUNC );
1666
1667
1668 //Loop through the SFT frequency bins and weight the data and normalize
1669 for ( UINT4 ii = 0; ii < numfbins; ii++ ) {
1670
1671 //scale antenna pattern weights by the right SFT frequency bin value in backgroundScaling**2
1672 fastSSVectorMultiply_with_stride_and_offset( backgroundScalingSq, backgroundScaling, backgroundScaling, numfbins, numfbins, ii, ii );
1673 memcpy( antweightssq->data, antweightssq0->data, sizeof( REAL4 )*antweightssq0->length );
1674 XLAL_CHECK( XLALVectorMultiplyREAL4( antweightssq->data, antweightssq->data, backgroundScalingSq->data, antweightssq->length ) == XLAL_SUCCESS, XLAL_EFUNC );
1675
1676 //get the background squared
1677 fastSSVectorMultiply_with_stride_and_offset( rngMeanssq, rngMeans, rngMeans, numfbins, numfbins, ii, ii );
1678
1679 //If noiseWeightOff is given, then set all the noise weights to be 1.0
1680 if ( params->noiseWeightOff ) for ( UINT4 jj = 0; jj < rngMeanssq->length; jj++ ) if ( rngMeanssq->data[jj] != 0.0 ) {
1681 rngMeanssq->data[jj] = 1.0;
1682 }
1683
1684 //Get sum of antenna pattern weight/variances for each frequency bin as a function of time (only for existant SFTs)
1685 REAL8 sumofweights = determineSumOfWeights( antweightssq, rngMeanssq );
1686 REAL8 invsumofweights = 1.0 / sumofweights;
1687
1688 //Now do noise weighting, antenna pattern weighting
1689 for ( UINT4 jj = 0; jj < indexValuesOfExistingSFTs->length; jj++ ) {
1690 output->data[indexValuesOfExistingSFTs->data[jj]*numfbins + ii] = ( REAL4 )( invsumofweights * antPatternWeights->data[indexValuesOfExistingSFTs->data[jj]] * tfdata->data[indexValuesOfExistingSFTs->data[jj] * numfbins + ii] / rngMeanssq->data[indexValuesOfExistingSFTs->data[jj]] );
1691 } /* for jj < indexValuesOfExisitingSFTs->length */
1692 } /* for ii < numfbins */
1693
1694 //Remember to reset the backgrnd vector to zero
1695 if ( params->signalOnly ) {
1696 memset( rngMeans->data, 0, sizeof( REAL4 )*rngMeans->length );
1697 }
1698
1699 //Destroy stuff
1700 XLALDestroyREAL4VectorAligned( antweightssq0 );
1701 XLALDestroyREAL4VectorAligned( antweightssq );
1702 XLALDestroyREAL4VectorAligned( rngMeanssq );
1703 XLALDestroyREAL4VectorAligned( backgroundScalingSq );
1704
1705 fprintf( stderr, "done\n" );
1706
1707 return XLAL_SUCCESS;
1708
1709} // tfWeight()
1710
1711/**
1712 * Determine the sum of the weights
1713 * \param [in] antweightssq Antenna pattern weights squared
1714 * \param [in] rngMeanssq Running mean values squared
1715 * \return Sum of the weights
1716 */
1717REAL8 determineSumOfWeights( const REAL4VectorAligned *antweightssq, const REAL4VectorAligned *rngMeanssq )
1718{
1719
1720 XLAL_CHECK_REAL8( antweightssq != NULL && rngMeanssq != NULL, XLAL_EINVAL );
1721
1722 REAL8 sumofweights = 0.0;
1723 for ( UINT4 ii = 0; ii < antweightssq->length; ii++ ) if ( rngMeanssq->data[ii] != 0.0 ) {
1724 sumofweights += antweightssq->data[ii] / rngMeanssq->data[ii];
1725 }
1726
1727 return sumofweights;
1728
1729} /* determineSumOfWeights */
1730
1731
1732/**
1733 * Determine if the SFTs are existing based on whether the first data point of power in each SFT is 0
1734 * \param [in] tfdata Pointer to REAL4VectorAligned of time-frequency data
1735 * \param [in] numffts Number of SFTs from the observation time
1736 * \return Pointer to INT4Vector containing 0 for non-present SFT or 1 for present SFT
1737 */
1738INT4Vector *existingSFTs( const REAL4VectorAligned *tfdata, const UINT4 numffts )
1739{
1740
1741 XLAL_CHECK_NULL( tfdata != NULL && numffts > 0, XLAL_EINVAL );
1742
1743 fprintf( stderr, "Determining existing SFTs... " );
1744
1745 UINT4 numfbins = tfdata->length / numffts; //Number of frequency bins
1746
1747 INT4Vector *sftexist = NULL;
1748 XLAL_CHECK_NULL( ( sftexist = XLALCreateINT4Vector( numffts ) ) != NULL, XLAL_EFUNC );
1749
1750 for ( UINT4 ii = 0; ii < numffts; ii++ ) {
1751 if ( tfdata->data[ii * numfbins] == 0.0 ) {
1752 sftexist->data[ii] = 0;
1753 } else {
1754 sftexist->data[ii] = 1;
1755 }
1756 }
1757
1758 fprintf( stderr, "done\n" );
1759
1760 return sftexist;
1761
1762} /* existingSFTs() */
1763
1764/**
1765 * Go through the backgroundScaling vector and zero out if the SFTexistVector has a 0 in that SFT location instead of 1
1766 * \param [in,out] background Pointer to REAL4VectorAligned of background data
1767 * \param [in,out] backgroundScaling Pointer to REAL4VectorAligned of background scaling data
1768 * \param [in] SFTexistVector Pointer to INT4Vector of 0 or 1 values representing if an SFT is present (1) or not (0)
1769 * \return Status value
1770 */
1771INT4 checkBackgroundScaling( const REAL4VectorAligned *background, const REAL4VectorAligned *backgroundScaling, const INT4Vector *SFTexistVector )
1772{
1773 XLAL_CHECK( backgroundScaling != NULL && SFTexistVector != NULL, XLAL_EINVAL );
1774
1775 UINT4 numfbins = backgroundScaling->length / SFTexistVector->length;
1776 for ( UINT4 ii = 0; ii < SFTexistVector->length; ii++ ) {
1777 if ( SFTexistVector->data[ii] == 0 && backgroundScaling->data[ii * numfbins] != 0.0 ) {
1778 memset( &( backgroundScaling->data[ii * numfbins] ), 0, sizeof( REAL4 )*numfbins );
1779 }
1780 if ( SFTexistVector->data[ii] == 0 && background->data[ii * numfbins] != 0.0 ) {
1781 memset( &( background->data[ii * numfbins] ), 0, sizeof( REAL4 )*numfbins );
1782 }
1783 }
1784 return XLAL_SUCCESS;
1785}
const REAL8 scaling[15]
REAL4VectorAligned * readInSFTs(UserInput_t *params, const REAL8 normalization, const REAL8 minfbin, const REAL8 maxfbin)
Read in the data SFTs in one function.
Definition: SFTfunctions.c:826
MultiLIGOTimeGPSVector * getMultiTimeStampsFromSFTs(const MultiSFTVector *multiSFTvector, const UserInput_t *params)
Create a list of timestamps from SFTs that might be a subset from those in an SFTCatalog,...
Definition: SFTfunctions.c:873
REAL8 determineSumOfWeights(const REAL4VectorAligned *antweightssq, const REAL4VectorAligned *rngMeanssq)
Determine the sum of the weights.
INT4 tfRngMeans(REAL4VectorAligned *output, const REAL4VectorAligned *tfdata, const UINT4 numffts, const UINT4 numfbins, const UINT4 blksize)
Determine the running mean of each SFT.
void removeBadSFTs(REAL4VectorAligned *tfdata, const INT4Vector *badsfts)
Remove the marked SFTs as bad by setting values to 0.
LIGOTimeGPSVector * jointTimestampsFromMultiTimestamps(const MultiLIGOTimeGPSVector *multiTimestamps)
Definition: SFTfunctions.c:981
REAL4TimeSeries * computeNSfreqTS(const PulsarParams *pulsarParams, LIGOTimeGPS tepoch, REAL8 Tsft, REAL8 SFToverlap, REAL8 duration)
Compute the GW frequency time series (in the SSB) for a set of assumed parameters.
Definition: SFTfunctions.c:689
INT4 replaceTFdataWithSubsequentTFdata(REAL4VectorAlignedArray *tfdataarray, const UINT4 numffts)
INT4Vector * markBadSFTs(const REAL4VectorAligned *tfdata, const UserInput_t *params)
Mark the non-Gaussian SFTs using K-S and Kuiper's tests.
SFTCatalog * findSFTdata(const UserInput_t *params)
Find the SFT data specified by user input.
Definition: SFTfunctions.c:35
INT4 tfMeanSubtract(REAL4VectorAligned *tfdata, const REAL4VectorAligned *rngMeans, const REAL4VectorAligned *backgroundScaling, const UINT4 numffts, const UINT4 numfbins)
Subtract running mean values from the SFT data, modifying input time-frequency data.
INT4 printSFTtimestamps2File(const MultiSFTVector *multiSFTvector, const UserInput_t *params)
MultiSFTVector * extractSFTband(const SFTCatalog *catalog, const REAL8 minfbin, const REAL8 maxfbin)
Extract the SFT coefficients from the band of interest.
Definition: SFTfunctions.c:76
REAL4VectorAligned * convertSFTdataToPowers(const SFTVector *sfts, const UserInput_t *params, const REAL8 normalization)
Convert a SFTVector of sfts into powers.
Definition: SFTfunctions.c:753
INT4 tfWeight(REAL4VectorAligned *output, const REAL4VectorAligned *tfdata, REAL4VectorAligned *rngMeans, REAL4VectorAligned *antPatternWeights, const REAL4VectorAligned *backgroundScaling, const INT4Vector *indexValuesOfExistingSFTs, const UserInput_t *params)
Weight the SFTs based on antenna pattern and noise variance (Equation 11, assuming the input time-fre...
MultiSFTVector * generateSFTdata(UserInput_t *uvar, const MultiLALDetector *detectors, const EphemerisData *edat, const INT4 maxbinshift, const gsl_rng *rng)
INT4 slideTFdata(REAL4VectorAligned *output, const UserInput_t *params, const REAL4VectorAligned *tfdata, const INT4Vector *binshifts)
Slide the time-frequency data to account for detector motion.
MultiLIGOTimeGPSVector * squeezeMultiTimestamps(const MultiLIGOTimeGPSVector *multiTS)
MultiSFTVector * getMultiSFTVector(UserInput_t *params, const REAL8 minfbin, const REAL8 maxfbin)
Get a MultiSFTVector from the user-input values.
Definition: SFTfunctions.c:102
INT4 checkBackgroundScaling(const REAL4VectorAligned *background, const REAL4VectorAligned *backgroundScaling, const INT4Vector *SFTexistVector)
Go through the backgroundScaling vector and zero out if the SFTexistVector has a 0 in that SFT locati...
MultiLIGOTimeGPSVector * getMultiTimeStampsFromSegmentsFile(const LALStringVector *filenames, const REAL8 t0, const REAL8 Tsft, const REAL8 SFToverlap, const REAL8 dur)
Create a list of timestamps from a segment list.
Definition: SFTfunctions.c:938
MultiLIGOTimeGPSVector * getMultiTimeStampsFromSFTCatalog(const SFTCatalog *catalog)
Create a list of timestamps from an SFTCatalog.
Definition: SFTfunctions.c:849
INT4Vector * existingSFTs(const REAL4VectorAligned *tfdata, const UINT4 numffts)
Determine if the SFTs are existing based on whether the first data point of power in each SFT is 0.
REAL4VectorAligned * coherentlyAddSFTs(const MultiSFTVector *multiSFTvector, const MultiSSBtimes *multissb, const MultiAMCoeffs *multiAMcoeffs, const LIGOTimeGPSVector *jointTimestamps, const REAL4VectorAlignedArray *backgroundRatio, const INT4 cosiSign, const assumeNSparams *NSparams, const UserInput_t *params, REAL4VectorAligned *backgroundScaling)
Add SFTs together from a MultiSFTVector.
Definition: SFTfunctions.c:147
LIGOTimeGPSVector * timestamps
#define fprintf
FILE * LOG
Definition: TwoSpect.c:55
INT4 kuipers_test_exp(REAL8 *kuipervalue, const REAL4VectorAligned *vector)
Kuiper's test of data against an expected exponential distribution.
INT4 ks_test_exp(REAL8 *ksvalue, const REAL4VectorAligned *vector)
KS test of data against an expected exponential distribution.
REAL4 calcMean(const REAL4VectorAligned *vector)
Compute the mean value of a REAL4VectorAligned, computed via recursion like in GSL.
void XLALDestroyREAL4VectorSequence(REAL4VectorSequence *vecseq)
PulsarParamsVector * XLALCreatePulsarParamsVector(UINT4 numPulsars)
Create zero-initialized PulsarParamsVector for numPulsars.
PulsarParamsVector * XLALPulsarParamsFromUserInput(const LALStringVector *UserInput, const LIGOTimeGPS *refTimeDef)
Function to determine the PulsarParamsVector input from a user-input defining CW sources.
void XLALDestroyPulsarParamsVector(PulsarParamsVector *ppvect)
Destructor for PulsarParamsVector type.
int XLALCWMakeFakeMultiData(MultiSFTVector **multiSFTs, MultiREAL8TimeSeries **multiTseries, const PulsarParamsVector *injectionSources, const CWMFDataParams *dataParams, const EphemerisData *edat)
Generate fake 'CW' data, returned either as SFTVector or REAL4Timeseries or both, for given CW-signal...
int XLALParseMultiNoiseFloor(MultiNoiseFloor *multiNoiseFloor, const LALStringVector *sqrtSX, UINT4 numDetectors)
Parse string-vectors (typically input by user) of N detector noise-floors for detectors ,...
int XLALGenerateSpinOrbitCW(PulsarCoherentGW *sourceSignal, SpinOrbitCWParamStruc *sourceParams)
FIXME: Temporary XLAL-wapper to LAL-function LALGenerateSpinOrbitCW()
#define LAL_LN2
#define LAL_PI
#define LAL_TWOPI
#define crect(re, im)
unsigned char BOOLEAN
#define LIGOTIMEGPSZERO
double complex COMPLEX16
#define cpolarf(r, th)
double REAL8
#define XLAL_INIT_DECL(var,...)
#define crectf(re, im)
char CHAR
uint32_t UINT4
float complex COMPLEX8
int32_t INT4
float REAL4
void * XLALCalloc(size_t m, size_t n)
void XLALFree(void *p)
void LALSRunningMedian2(LALStatus *status, REAL4Sequence *medians, const REAL4Sequence *input, LALRunningMedianPar param)
REAL8 XLALRngMedBias(INT4 blkSize)
void XLALDestroySFTVector(SFTVector *vect)
XLAL interface to destroy an SFTVector.
Definition: SFTtypes.c:300
void XLALDestroySFTCatalog(SFTCatalog *catalog)
Free an 'SFT-catalogue'.
Definition: SFTcatalog.c:329
MultiSFTVector * XLALLoadMultiSFTs(const SFTCatalog *inputCatalog, REAL8 fMin, REAL8 fMax)
Function to load a catalog of SFTs from possibly different detectors.
Definition: SFTfileIO.c:416
void XLALDestroyMultiSFTCatalogView(MultiSFTCatalogView *multiView)
Destroys a MultiSFTCatalogView, without freeing the original catalog that the 'view' was referring to...
Definition: SFTcatalog.c:496
void XLALDestroyMultiSFTVector(MultiSFTVector *multvect)
Destroy a multi SFT-vector.
Definition: SFTtypes.c:424
MultiLIGOTimeGPSVector * XLALMakeMultiTimestamps(LIGOTimeGPS tStart, REAL8 Tspan, REAL8 Tsft, REAL8 Toverlap, UINT4 numDet)
Same as XLALMakeTimestamps() just for several detectors, additionally specify the number of detectors...
int XLALExtractStrictBandFromSFT(SFTtype **outSFT, const SFTtype *inSFT, REAL8 fMin, REAL8 Band)
Return a copy of an SFT containing only the bins in [fMin, fMin+Band).
Definition: SFTtypes.c:601
MultiSFTCatalogView * XLALGetMultiSFTCatalogView(const SFTCatalog *catalog)
Return a MultiSFTCatalogView generated from an input SFTCatalog.
Definition: SFTcatalog.c:380
int XLALMultiSFTVectorResizeBand(MultiSFTVector *multiSFTs, REAL8 f0, REAL8 Band)
Resize the frequency-band of a given multi-SFT vector to [f0, f0+Band].
Definition: SFTtypes.c:946
MultiLIGOTimeGPSVector * XLALTimestampsFromMultiSFTCatalogView(const MultiSFTCatalogView *multiView)
Given a multi-SFTCatalogView, return a MultiLIGOTimeGPSVector holding the SFT timestamps.
MultiSFTVector * XLALExtractMultiSFTVectorWithMultiTimestamps(const MultiSFTVector *multiSFTs, const MultiLIGOTimeGPSVector *multiTimestamps)
Extract a MultiSFTVector from another MultiSFTVector but only those timestamps matching.
Definition: SFTtypes.c:1116
LIGOTimeGPSVector * XLALCreateTimestampVector(UINT4 len)
Allocate a LIGOTimeGPSVector.
Definition: SFTtimestamps.c:47
int XLALSFTAdd(SFTtype *a, const SFTtype *b)
Adds SFT-data from SFT 'b' to SFT 'a'.
Definition: SFTtypes.c:775
int XLALCopySFT(SFTtype *dest, const SFTtype *src)
Copy an entire SFT-type into another.
Definition: SFTtypes.c:202
LIGOTimeGPSVector * XLALTimestampsFromSegmentFile(const char *filename, REAL8 Tsft, REAL8 Toverlap, BOOLEAN adjustSegExtraTime, BOOLEAN synchronize)
Extract timestamps-vector from a segment file, with functionality based on MakeSFTDAG The filename sh...
SFTCatalog * XLALSFTdataFind(const CHAR *file_pattern, const SFTConstraints *constraints)
Find the list of SFTs matching the file_pattern and satisfying the given constraints,...
Definition: SFTcatalog.c:71
void XLALDestroySFT(SFTtype *sft)
Destructor for one SFT.
Definition: SFTtypes.c:176
int XLALReorderMultiSFTVector(MultiSFTVector *multiSFTs, const LALStringVector *IFOs)
Reorder the MultiSFTVector with specified list of IFOs.
Definition: SFTtypes.c:738
MultiLIGOTimeGPSVector * XLALExtractMultiTimestampsFromSFTs(const MultiSFTVector *multiSFTs)
Given a multi-SFT vector, return a MultiLIGOTimeGPSVector holding the SFT timestamps.
int XLALMultiSFTVectorAdd(MultiSFTVector *a, const MultiSFTVector *b)
Adds SFT-data from MultiSFTvector 'b' to elements of MultiSFTVector 'a'.
Definition: SFTtypes.c:842
void XLALDestroyMultiTimestamps(MultiLIGOTimeGPSVector *multiTS)
Destroy a MultiLIGOTimeGPSVector timestamps vector.
SFTVector * XLALCreateSFTVector(UINT4 numSFTs, UINT4 numBins)
XLAL function to create an SFTVector of numSFT SFTs with SFTlen frequency-bins (which will be allocat...
Definition: SFTtypes.c:230
void XLALDestroyTimestampVector(LIGOTimeGPSVector *vect)
De-allocate a LIGOTimeGPSVector.
Definition: SFTtimestamps.c:69
MultiLIGOTimeGPSVector * XLALReadMultiTimestampsFiles(const LALStringVector *fnames)
backwards compatible wrapper to XLALReadMultiTimestampsFilesConstrained() without GPS-time constraint...
int XLALSinCosLUT(REAL4 *sinx, REAL4 *cosx, REAL8 x)
Calculate sin(x) and cos(x) to roughly 1e-7 precision using a lookup-table and Taylor-expansion.
Definition: SinCosLUT.c:83
COORDINATESYSTEM_EQUATORIAL
REAL4TimeSeries * XLALCreateREAL4TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
void XLALDestroyREAL8TimeSeries(REAL8TimeSeries *series)
void XLALDestroyREAL4TimeSeries(REAL4TimeSeries *series)
const LALUnit lalHertzUnit
int XLALUserVarWasSet(const void *cvar)
void XLALDestroyUINT4Vector(UINT4Vector *vector)
INT4Vector * XLALCreateINT4Vector(UINT4 length)
COMPLEX8Vector * XLALCreateCOMPLEX8Vector(UINT4 length)
void XLALDestroyINT4Vector(INT4Vector *vector)
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
void XLALDestroyREAL8Vector(REAL8Vector *vector)
UINT4Vector * XLALCreateUINT4Vector(UINT4 length)
void XLALDestroyCHARVector(CHARVector *vector)
CHARVector * XLALCreateCHARVector(UINT4 length)
void XLALDestroyCOMPLEX8Vector(COMPLEX8Vector *vector)
int XLALVectorSinCosREAL4(REAL4 *out1, REAL4 *out2, const REAL4 *in, const UINT4 len)
int XLALVectorScaleREAL4(REAL4 *out, REAL4 scalar, const REAL4 *in, const UINT4 len)
void XLALDestroyREAL4VectorAligned(REAL4VectorAligned *in)
REAL4VectorAligned * XLALCreateREAL4VectorAligned(const UINT4 length, const UINT4 align)
int XLALVectorSinCos2PiREAL4(REAL4 *out1, REAL4 *out2, const REAL4 *in, const UINT4 len)
int XLALVectorAddREAL4(REAL4 *out, const REAL4 *in1, const REAL4 *in2, const UINT4 len)
int XLALVectorMultiplyREAL4(REAL4 *out, const REAL4 *in1, const REAL4 *in2, const UINT4 len)
int XLALVectorShiftREAL4(REAL4 *out, REAL4 scalar, const REAL4 *in, const UINT4 len)
#define xlalErrno
#define XLAL_CHECK(assertion,...)
#define XLAL_CHECK_VOID(assertion,...)
#define XLAL_CHECK_REAL8(assertion,...)
#define XLAL_CHECK_NULL(assertion,...)
XLAL_EBADLEN
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EIO
XLAL_EINVAL
XLAL_EFAILED
LIGOTimeGPS * XLALGPSSetREAL8(LIGOTimeGPS *epoch, REAL8 t)
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
int XLALGPSCmp(const LIGOTimeGPS *t0, const LIGOTimeGPS *t1)
REAL8 XLALGPSGetREAL8(const LIGOTimeGPS *epoch)
end
out
LALDetector detectors[LAL_NUM_DETECTORS]
double t0
REAL4Vector * b
(weighted) per-SFT antenna-pattern function
Definition: LALComputeAM.h:65
REAL4Vector * a
(weighted) per-SFT antenna-pattern function
Definition: LALComputeAM.h:64
CHAR * data
UINT4 length
COMPLEX8Sequence * data
A vector of COMPLEX8FrequencySeries.
COMPLEX8FrequencySeries * data
Pointer to the data array.
UINT4 length
Number of elements in array.
COMPLEX8 * data
Struct controlling all the aspects of the fake data (time-series + SFTs) to be produced by XLALCWMake...
This structure contains all information about the center-of-mass positions of the Earth and Sun,...
INT4 * data
UINT4 length
INT4 gpsNanoSeconds
A vector of 'timestamps' of type LIGOTimeGPS.
Definition: SFTfileIO.h:188
REAL8 deltaT
'length' of each timestamp (e.g.
Definition: SFTfileIO.h:194
LIGOTimeGPS * data
array of timestamps
Definition: SFTfileIO.h:193
UINT4 length
number of timestamps
Definition: SFTfileIO.h:192
Multi-IFO container for antenna-pattern coefficients and atenna-pattern matrix .
Definition: LALComputeAM.h:137
AMCoeffs ** data
noise-weighted AM-coeffs , and
Definition: LALComputeAM.h:142
array of detectors definitions 'LALDetector'
A collection of (multi-IFO) LIGOTimeGPSVector time-stamps vectors.
Definition: SFTfileIO.h:198
UINT4 length
number of timestamps vectors or ifos
Definition: SFTfileIO.h:202
LIGOTimeGPSVector ** data
timestamps vector for each ifo
Definition: SFTfileIO.h:203
array of detector-specific 'noise floors' (ie PSD values), assumed constant over the frequency-band o...
REAL8 sqrtSn[PULSAR_MAX_DETECTORS]
per-IFO sqrt(PSD) values , where
A multi-SFT-catalogue "view": a multi-IFO vector of SFT-catalogs.
Definition: SFTfileIO.h:255
A collection of SFT vectors – one for each IFO in a multi-IFO search.
Definition: SFTfileIO.h:179
UINT4 length
number of ifos
Definition: SFTfileIO.h:183
SFTVector ** data
sftvector for each ifo
Definition: SFTfileIO.h:184
Multi-IFO container for SSB timings.
Definition: SSBtimes.h:67
SSBtimes ** data
array of SSBtimes (pointers)
Definition: SSBtimes.h:72
UINT4 length
number of IFOs
Definition: SSBtimes.h:71
REAL8 aCross
Signal amplitude (cross)
REAL8 psi
polarization angle psi
REAL8 aPlus
Signal amplitude (plus)
REAL8 phi0
initial signal-phase (at some reference time)
This structure stores a representation of a plane gravitational wave propagating from a particular po...
REAL8 period
Binary: orbital period (sec)
LIGOTimeGPS tp
Binary: time of observed periapsis passage (in SSB)
PulsarSpins fkdot
Intrinsic spins: [Freq, f1dot, f2dot, ... ] where fkdot = d^kFreq/dt^k.
REAL8 Delta
Sky position: DEC (latitude) in equatorial coords and radians.
LIGOTimeGPS refTime
Reference time of pulsar parameters (in SSB!)
REAL8 Alpha
Sky position: RA (longitude) in equatorial coords and radians.
REAL8 ecc
Binary: orbital eccentricity.
REAL8 asini
Binary: projected, normalized orbital semi-major axis (s).
REAL8 argp
Binary: argument of periapsis (radians)
Type defining the parameters of a pulsar-source of CW Gravitational waves.
PulsarAmplitudeParams Amp
'Amplitude-parameters': h0, cosi, phi0, psi
PulsarDopplerParams Doppler
'Phase-evolution parameters': {skypos, fkdot, orbital params }
Straightforward vector type of N PulsarParams structs.
PulsarParams * data
array of pulsar-signal parameters
REAL4Sequence * data
REAL4VectorAligned ** data
REAL4 * data
REAL8 * data
An "SFT-catalogue": a vector of SFTdescriptors, as returned by XLALSFTdataFind()
Definition: SFTfileIO.h:238
'Constraints' for SFT-matching: which detector, within which time-stretch and which timestamps exactl...
Definition: SFTfileIO.h:212
REAL8Vector * Tdot
dT/dt : time-derivative of SSB-time wrt local time for SFT-alpha
Definition: SSBtimes.h:63
REAL8Vector * DeltaT
Time-difference of SFT-alpha - tau0 in SSB-frame.
Definition: SSBtimes.h:62
LIGOTimeGPS refTime
reference-time 'tau0'
Definition: SSBtimes.h:61
REAL8 longitude
REAL8 latitude
This structure stores the parameters for constructing a gravitational waveform with both a Taylor-pol...
UINT4 * data
alignedREAL8Vector ** data
REAL8 * assumeNSpsi
Definition: SFTfunctions.h:35
REAL8 * assumeNScosi
Definition: SFTfunctions.h:34
LIGOTimeGPS * assumeNSorbitTp
Definition: SFTfunctions.h:39
LIGOTimeGPS * assumeNSrefTime
Definition: SFTfunctions.h:40
REAL8 * assumeNSasini
Definition: SFTfunctions.h:38
SkyPosition assumeNSpos
Definition: SFTfunctions.h:33
REAL8 * assumeNSGWfreq
Definition: SFTfunctions.h:36
REAL8 * assumeNSorbitP
Definition: SFTfunctions.h:37
INT4 VectorShiftREAL8(alignedREAL8Vector *output, alignedREAL8Vector *input, REAL8 shift, INT4 vectorMath)
Definition: vectormath.c:269
void destroyAlignedREAL8Vector(alignedREAL8Vector *vector)
Definition: vectormath.c:67
INT4 DirichletRatioVector(COMPLEX8Vector *output, alignedREAL8Vector *delta0, alignedREAL8Vector *delta1, alignedREAL8Vector *scaling, const UserInput_t *params)
Definition: vectormath.c:122
INT4 VectorCabsfCOMPLEX8(REAL4VectorAligned *output, COMPLEX8Vector *input, INT4 vectorMath)
Definition: vectormath.c:533
void destroyAlignedREAL8VectorArray(alignedREAL8VectorArray *array)
Definition: vectormath.c:88
INT4 VectorRoundREAL8(alignedREAL8Vector *output, alignedREAL8Vector *input, INT4 vectorMath)
Definition: vectormath.c:396
INT4 VectorAddREAL8(alignedREAL8Vector *output, alignedREAL8Vector *input1, alignedREAL8Vector *input2, INT4 vectorMath)
Definition: vectormath.c:282
alignedREAL8VectorArray * createAlignedREAL8VectorArray(const UINT4 length, const UINT4 vectorLength, const size_t align)
Definition: vectormath.c:77
INT4 VectorFloorREAL8(alignedREAL8Vector *output, alignedREAL8Vector *input, INT4 vectorMath)
Definition: vectormath.c:334
INT4 VectorScaleREAL8(alignedREAL8Vector *output, alignedREAL8Vector *input, REAL8 scale, INT4 vectorMath)
Definition: vectormath.c:256
INT4 VectorAbsREAL8(alignedREAL8Vector *output, alignedREAL8Vector *input, INT4 vectorMath)
Definition: vectormath.c:480
INT4 fastSSVectorMultiply_with_stride_and_offset(REAL4VectorAligned *output, const REAL4VectorAligned *input1, const REAL4VectorAligned *input2, const INT4 stride1, const INT4 stride2, const INT4 offset1, const INT4 offset2)
Computes a multiplication of two vectors with a stride and initial offset.
Definition: vectormath.c:221
alignedREAL8Vector * createAlignedREAL8Vector(UINT4 length, const size_t align)
Definition: vectormath.c:55
INT4 VectorSubtractREAL8(alignedREAL8Vector *output, alignedREAL8Vector *input1, alignedREAL8Vector *input2, INT4 vectorMath)
Definition: vectormath.c:295
INT4 VectorMultiplyREAL8(alignedREAL8Vector *output, alignedREAL8Vector *input1, alignedREAL8Vector *input2, INT4 vectorMath)
Definition: vectormath.c:308
INT4 VectorCabsCOMPLEX8(alignedREAL8Vector *output, COMPLEX8Vector *input, INT4 vectorMath)
Definition: vectormath.c:591