cw-data-analysis (OctApps Reference Manual)

2.5 `cw-data-analysis`

Function File: Amu = amplitudeParams2Vect ( Amp )

compute the amplitude-vector {A^mu} for given amplitude-params, which can follow

either the MLDC convention {Amplitude, Inclination, Polarization, IntialPhase },
or in LIGO convention {h0, cosi, psi, phi0}: this will be auto-detected and properly converted.

multiple signals must correspond to different *lines* in those fields, i.e. column-vectors! the output consists of 4D line vectors Amu(,1:4), multiple lines corresponding to multiple signals

Examples

p0 = struct("h0", 1e-24, "cosi", 0, "psi", 0.33 * pi/4, "phi0", pi/5);
p = amplitudeVect2Params(amplitudeParams2Vect(p0));
assert(p.h0, p0.h0, 1e-3);
assert(p.cosi, p0.cosi, 1e-3);
assert(p.psi, p0.psi, 1e-3);
assert(p.phi0, p0.phi0, 1e-3);

Function File: Amp = amplitudeVect2Params ( Amu, convention )

compute amplitude-vector {A^mu} from (MLDC) amplitudes {Amplitude, Inclination, Polarization, InitialPhase }

Amu is a row-vector for each signal, multiple signals being stored in multiple rows , the resulting fields in Amp are also column-vectors for multiple signals

if convention == "LIGO", return {h0,cosi,psi.iota} and {aPlus,aCross},
if convention == "MLDC" return {Amplitude,Inclination,Polarization, InitialPhase} using MLDC conventions
the default = "LIGO" if not specified

Note

Adapted from algorithm in LALEstimatePulsarAmplitudeParams()

Examples

A0 = 1e-24 * randn(1, 4);
sert(amplitudeParams2Vect(amplitudeVect2Params(A0, "LIGO")), A0, 1e-3)
sert(amplitudeParams2Vect(amplitudeVect2Params(A0, "MLDC")), A0, 1e-3)

Function File: [ convention, numSignals ] = checkAmplitudeParams ( Amp )

check syntactic correctness of amplitude-parameter struct, and determine its convention: "LIGO" || "MLDC", depending on whether fields {Amplitude, Inclination, Polarization, InitialPhase}, or {h0, cosi, psi, phi0} are present. The presence of both types of fields is an error.

Examples

assert(checkAmplitudeParams(struct("h0", 1e-24, "cosi", 0, "psi", pi/4, "phi0", pi/5)), "LIGO")

Function File: results = DoFstatInjections ( opt, val, … )

Perform full software injections in generated SFTs using LALPulsar functions.

Arguments

result: results structure

Options (starred options are returned in results)

ref_time: reference time in GPS seconds
start_time: start time in GPS seconds (default: ref_time - 0.5*time_span)
time_span: observation time-span in seconds
detectors: comma-separated list of detector names
det_sqrt_sqrtSX: sqrt(single-sided noise sqrtSX) to assume for each detector
ephemerides: Earth/Sun ephemerides from loadEphemerides()
sft_time_span: SFT time-span in seconds (default: 1800)
sft_overlap: SFT overlap in seconds (default: 0)
sft_noise_window: number of bins used when estimating SFT noise (default: 50)
inj_sqrt_sqrtSX: inject Gaussian random noise with sqrt(single-sided noise sqrtSX) for each detector
*inj_h0: injected h0 strain amplitude (default: 1.0)
*inj_cosi: injected cosine of inclination angle (default: random)
*inj_psi: injected polarisation angle (default: random)
*inj_phi0: injected initial phase (default: random)
*inj_alpha: injected right ascension (default: random)
*inj_delta: injected declination (default: random)
*inj_fndot: injected frequency/spindowns (default: 100 Hz)
*sch_alpha: searched right ascension (default: same as injected)
*sch_delta: searched declination (default: same as injected)
*sch_fndot: searched frequency/spindowns (default: same as injected)
OrbitParams: option that needs to be set to ’true’ to be able to specify the orbital parameters (default: false)
*inj_orbitasini: injected orbital projected semi-major axis (normalised by the speed of light) in seconds (default: random)
*inj_orbitEcc: injected orbital eccentricity (default: random)
*inj_orbitTpSSB: injected (SSB) time of periapsis passage (in seconds) (default: random)
*inj_orbitPeriod: injected orbital period (seconds) (default: random)
*inj_orbitArgp: injected argument of periapsis (radians) (default: random)
*sch_orbitasini: searched orbital projected semi-major axis (normalised by the speed of light) in seconds (default: same as injected)
*sch_orbitEcc: searched orbital eccentricity (default: same as injected)
*sch_orbitTpSSB: searched (SSB) time of periapsis passage (in seconds) (default: same as injected)
*sch_orbitPeriod: searched orbital period (seconds) (default: same as injected)
*sch_orbitArgp: searched argument of periapsis (radians) (default: same as injected)
Dterms: number of Dirichlet terms to use in ComputeFstat() (default: number used by optimised hotloops)

Examples

See the tutorial on DoFstatInjections().

Function File: [ p, v, sp, sv, op, ov ] = getDetectorPosVel ( opt, val, … )

Compute the spin, orbital, and total components of a detector’s position and velocity at a list of GPS times, using LALPulsar

Arguments

p: detector positions, in equatorial coordinates
v: detector velocities, in equatorial coordinates
sp: spin components of detector positions
sv: spin components of detector velocities
op: orbital components of detector positions
ov: orbital components of detector velocities

Options

gps_times: list of GPS times
detector: name of detector (default: H1)
motion: type of motion (default: spin+orbit)
ephemerides: Earth/Sun ephemerides from loadEphemerides()

Examples

try
  lal; lalpulsar;
catch
  disp("skipping test: LALSuite bindings not available"); return;
end_try_catch
[p, v, sp, sv, op, ov] = getDetectorPosVel("gps_times", 800000000 + 100*(0:5));
assert(p, sp + op, 1e-3);
assert(v, sv + ov, 1e-3);

Function File: results = injectionRecoveryGCT ( opt, val, … )

Perform signal injection and (area-search) recovery using HierarchSearchGCT

Arguments

results: structure containing various histograms of measured statistics and mismatches from the injection+recovery runs, and

Options

Ntrials: (optional) number of repeated injection+recovery trials to perform [default: 1]
timestampsFiles: CSV list of SFT timestamp filenames
IFOs: CSV list of IFO names (eg "H1,L1,...")
segmentList: filename of segment list (containing lines of the form "startGPS endGPS\n")
inj_sqrtSX: injections: (optional) CSV list of per-detector noise-floor sqrt(PSD) to generate
inj_h0: injections: signal amplitude ’h0’ of signals
inj_SNR: injections: alternative: signal-to-noise ratio ’SNR’ of signals
inj_AlphaRange: injections: range of sky-position alpha to (isotropically) draw from [default: [0, 2pi]]
inj_DeltaRange: injections: range of sky-position delta to (isotropically) draw from [default: [-pi/2, pi/2]]
inj_FreqRange: injections: range of signal frequencies to draw from
inj_fkdotRange: injections: [numSpindowns x 2] ranges of spindown-values to draw from [default: []]
dFreq: search: frequency resolution
dfkdot: search: numSpindowns vector of spindown resolutions to use in search
gammaRefine: search: numSpindowns vector of ’gammeRefine[s]’ refinement factors to use
skyGridFile: search: sky-grid file to use
sch_Nsky: search-box: number of nearest-neighbor skygrid points to use around injection
sch_Nfreq: search-box: number of frequency bins to use around injection frequency
sch_Nfkdot: search-box: number of spindown-bins to use around injection spindown-value
FstatMethod: search: F-statistic method to use: "DemodBest", "ResampBest", ...
computeBSGL: search: additionally compute and histogram B_S/GL statistic values
Fstar0: search: BSGL parameter ’Fstar0sc’
nCand: search: number of toplist candidates to keep
GCT_binary: which GCT executable to use for searching
cleanup: boolean: remove intermediate output files at the end or not

Examples

if isempty(file_in_path(getenv("PATH"), "lalapps_HierarchSearchGCT"))
  disp("skipping test: LALApps programs not available"); return;
endif
output = nthargout(2, @system, "lalapps_HierarchSearchGCT --version");
LALApps_version = versionstr2hex(nthargout(5, @regexp, output, "LALApps: ([0-9.]+)"){1}{1,1});
if LALApps_version <= 0x06210000
  disp("cannot run test as version of lalapps_HierarchSearchGCT is too old"); return;
endif
oldpwd = pwd;
basedir = mkpath(tempname(tempdir));
unwind_protect
  cd(basedir);
  args = struct;
  args.timestampsFiles = "H1.txt";
  args.IFOs = "H1";
  args.segmentList = "segs.txt";
  args.inj_sqrtSX = 1.0;
  args.inj_h0 = 1.0;
  args.inj_AlphaRange = [0, 2*pi];
  args.inj_DeltaRange = [-pi/2, pi/2];
  args.inj_FreqRange = [100, 100.01];
  args.inj_fkdotRange = [-1e-8, 0];
  args.dFreq = 1e-7;
  args.dfkdot = 1e-11;
  args.gammaRefine = 100;
  args.skyGridFile = "sky.txt";
  args.sch_Nsky = 5;
  args.sch_Nfreq = 5;
  args.sch_Nfkdot = 5;
  args.FstatMethod = "DemodBest";
  args.cleanup = true;
  fid = fopen(args.timestampsFiles, "w");
  for i = 1:10
    fprintf(fid, "%i\n", 800000000 + 1800*i);
  endfor
  fclose(fid);
  fid = fopen(args.segmentList, "w");
  fprintf(fid, "%i %i\n", 800000000 + 1800*[0, 5]);
  fprintf(fid, "%i %i\n", 800000000 + 1800*[5, 10]);
  fclose(fid);
  fid = fopen(args.skyGridFile, "w");
  for i = 1:50
    fprintf(fid, "%.5f %.5f\n", octforge_unifrnd(0+eps, 2*pi-eps), octforge_unifrnd(-pi/2+eps, pi/2-eps));
  endfor
  fclose(fid);
  fevalstruct(@injectionRecoveryGCT, args);
unwind_protect_cleanup
  cd(oldpwd);
end_unwind_protect

Function File: loadCandidateFile ( fname )

loads a ’candidate-file’ from lalapps_ComputeFStatistic_v2 --outputLoudest=cand.file and returns a struct containing the data

Examples

if isempty(file_in_path(getenv("PATH"), "lalapps_ComputeFstatistic_v2"))
  disp("skipping test: LALApps programs not available"); return;
endif
output = nthargout(2, @system, "lalapps_ComputeFstatistic_v2 --version");
LALApps_version = versionstr2hex(nthargout(5, @regexp, output, "LALApps: ([0-9.]+)"){1}{1,1});
if LALApps_version <= 0x06210000
  disp("cannot run test as version of lalapps_ComputeFstatistic_v2 is too old"); return;
endif
args = struct;
args.Alpha = 1.2;
args.Delta = 3.4;
args.Freq = 100;
args.f1dot = -1e-8;
args.Tsft = 1800;
args.IFOs = "H1";
args.injectSqrtSX = 1.0;
args.timestampsFiles = tempname(tempdir);
args.outputLoudest = tempname(tempdir);
fid = fopen(args.timestampsFiles, "w");
fprintf(fid, "800000000\n800001800\n800003600\n");
fclose(fid);
runCode(args, "lalapps_ComputeFstatistic_v2");
cand_file = loadCandidateFile(args.outputLoudest);
assert(isstruct(cand_file));

Function File: ephemerides = loadEphemerides ( opt, val, … )

Load Earth and Sun ephemerides from LALPulsar.

Arguments

ephemerides: structure containing ephemerides

Options

earth_file: Earth ephemerides file (default: earth00-40-DE405.dat.gz)
sun_file: Sun ephemerides file (default: sun00-40-DE405.dat.gz)

Examples

try
  lal; lalpulsar;
catch
  disp("skipping test: LALSuite bindings not available"); return;
end_try_catch
ephemerides = loadEphemerides();
ephemerides = loadEphemerides();
ephemerides = loadEphemerides();

Function File: [ resampInfo, demodInfo ] = predictFstatTimeAndMemory ( varargin )

Predict single-segment F-statistic computation time per frequency bin per detector (tauF_core and tauF_buffer) and corresponding memory requirements (MBWorkspace, MBDataPerDetSeg) for both Resampling and Demod Fstat methods.

See the F-stat timing notes at https://dcc.ligo.org/LIGO-T1600531-v4 for a detailed description of the F-statistic timing model and notation.

Note

The estimate is for one coherent segment of length Tcoh, whileTspan is only used to correctly deal with the GCT code’s handling of multi-segment searches, which can affect timing and memory requirements for each segment. In the case of Weave, however, use Tspan = Tcoh [default].

Input parameters

Tcoh: coherent segment length
Tspan: total data time-span [Default: Tcoh]. Note: this is used to estimate the total memory in the case of the GCT search code, which load the full frequency band of data over all segments. In the case of the Weave code, the SFT frequency band is computed for each segment separately, so in this case one should use Tspan==Tcoh to correctly estimate the timing and memory!
Freq0: start search frequency
FreqBand: search frequency band
dFreq: search frequency spacing

Optional arguments

f1dot0
f1dotBand: first-order spindown range [f1dot0, f1dot0+f1dotBand] [Default: 0]
f2dot0
f2dotBand: 2nd-order spindown range [f2dot0,f2dot0+f2dotBand] [Default: 0]
Dterms: number of Dterms used in sinc-interpolation [Default: 8]
Nsft: number of SFTs (for single segment, single detector) [Default: Nsft=Tcoh/Tsft]
Tsft: SFT length [Default: 1800]
refTimeShift: offset of reference time from starttime, measured in units of Tspan, ie refTimeShift = (refTime - startTime)/Tspan [Default: 0.5]
binaryMaxAsini: Maximum projected semi-major axis a*sini/c (= 0 for isolated sources) [Default: 0]
binaryMinPeriod: Minimum orbital period (s); must be 0 for isolated signals [Default: 0]
binaryMaxEcc: Maximal binary eccentricity: must be 0 for isolated signals [Default: 0]
resampFFTPowerOf2: enforce FFT length to be a power of two (by rounding up) [Default: true]

Resampling timing model coefficients

tau0_Fbin: Resampling timing coefficient for contributions scaling with output frequency bins NFbin
tau0_FFT: Resampling timing coefficient for FFT performance. Can be 2-element vector [t1, t2]: use t1 if log2(NsampFFT) <= 18, t2 otherwise
tau0_spin: Resampling timing coefficient for applying spin-down corrections
tau0_bary: Resampling timing coefficient (buffered) barycentering (contributes to tauF_buffer)

Demod timing model coefficients

tau0_coreLD: Demod timing coefficient for core F-stat time
tau0_bufferLD: Demod timing coefficient for computation of buffered quantities

Return values

Two structs resampInfo and demodInfo with fields:

tauF_core

Fstat time per frequency bin per detector excluding time to compute buffered quantities (eg barycentering) [in seconds]

tauF_buffer

Fstat time per frequency bin per detector for computing all the buffered quantities once [in seconds] The effective F-stat time per frequency bin per detector is therefore tauF_eff = tauF_core + b * tauF_buffer, where b = 1/N_{f1,f2,...} is the fraction of calls to XLALComputeFstat() where the buffer can be re-used

l2NsampFFT

resamp only log_2(NsampFFT) where NsampFFT is the number of FFT samples

MBWorkspace

resamp only memory for (possibly shared) workspace [in MBytes]

MBDataPerDetSeg

memory to hold all data per detector, per-segment (original+buffered) [in MBytes] ie total data memory would be

  memData[all] = Nseg * Ndet * memDataPerDetSeg

Examples

try
  lal; lalpulsar;
catch
  disp("skipping test: LALSuite bindings not available"); return;
end_try_catch
[resampInfo, demodInfo] = predictFstatTimeAndMemory("Tcoh", 86400, "Freq0", 100, "FreqBand", 1e-2, "dFreq", 1e-7);
assert(isstruct(resampInfo));
assert(isfield(resampInfo, "tauF_core"));
assert(resampInfo.tauF_core > 0);
assert(isfield(resampInfo, "tauF_buffer"));
assert(resampInfo.tauF_buffer > 0);
assert(isfield(resampInfo, "MBDataPerDetSeg"));
assert(resampInfo.MBDataPerDetSeg > 0);
assert(isstruct(demodInfo));
assert(isfield(demodInfo, "tauF_core"));
assert(demodInfo.tauF_core > 0);
assert(isfield(demodInfo, "tauF_buffer"));
assert(demodInfo.tauF_buffer > 0);
assert(isfield(demodInfo, "MBDataPerDetSeg"));
assert(demodInfo.MBDataPerDetSeg > 0);

Function File: ret = randSignalParams ( ranges, [ numSignals ] )

generate random-parameters for ’numSignals’ (default=1, returns a colunm-vector) signals within given ranges and return the signal-parameters in a struct sigparams = [h0, cosi, psi, phi0, alpha, delta, f, f1dot, f2dot, f3dot]

Examples

assert(isstruct(randSignalParams(struct("h0", 1e-24, "cosi", 0, "psi", pi/4, "phi0", pi/5))))

Function File: ret = readSFT ( fname )

read a given SFT-file and return its meta-info (header) and data as a struct: ret = {version; epoch; Tsft; f0; Band; SFTdata }

C-type of SFTs

typedef struct tagSFTHeader {
   REAL8  version;              /* SFT version-number (currently only 1.0 allowed )*/
   INT4   gpsSeconds;           /* gps start-time */
   INT4   gpsNanoSeconds;
   REAL8  timeBase;             /* length of data-stretch in seconds */
   INT4   fminBinIndex;         /* first frequency-index contained in SFT */
   INT4   length;               /* number of frequency bins */

   /* v2-specific part: */
   INT8 crc64;                  /* 64 bits */
   CHAR detector[2];
   CHAR padding[2];
   INT comment_length;
} SFTHeader;
CHAR[comment_length] comment;

Examples

sft = readSFT(fullfile(fileparts(file_in_loadpath("readSFT.m")), "SFT-good"));

2.5 cw-data-analysis

Examples

Note

Examples

Examples

Arguments

Options (starred options are returned in results)

Examples

Arguments

Options

Examples

Arguments

Options

Examples

Examples

Arguments

Options

Examples

Note

Input parameters

Optional arguments

Resampling timing model coefficients

Demod timing model coefficients

Return values

Examples

Examples

C-type of SFTs

Examples

2.5 `cw-data-analysis`