Processing math: 25%
LALApps 10.1.0.1-3a66518
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
inspinj.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2007 Chad Hanna, Alexander Dietz, Duncan Brown, Gareth Jones, Jolien Creighton, Nickolas Fotopoulos, Patrick Brady, Stephen Fairhurst, Tania Regimbau, Salvatore Vitale
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
21
22/*-----------------------------------------------------------------------
23 *
24 * File Name: inspinj.c
25 *
26 * Author: Brown, D. A., Creighton, J. D. E. and Dietz A. IPN contributions from Predoi, V.
27 *
28 *
29 *-----------------------------------------------------------------------
30 */
31
32/**
33 * \file
34 * \ingroup lalapps_inspiral
35 *
36 * <dl>
37 * <dt>Name</dt><dd>
38 * \c lalapps_inspinj --- produces inspiral injection data files.</dd>
39 *
40 * <dt>Synopsis</dt><dd>
41 * \c lalapps_inspinj
42 *
43 * [<tt>--help</tt>]
44 * <tt>--source-file</tt> \c sfile
45 * <tt>--mass-file</tt> \c mfile
46 *
47 * [<tt>--gps-start-time</tt> \c tstart]
48 * [<tt>--gps-end-time</tt> \c tend]
49 *
50 * [<tt>--time-step</tt> \c tstep]
51 * [<tt>--time-interval</tt> \c tinterval]
52 *
53 * [<tt>--seed</tt> \c seed]
54 * [<tt>--waveform</tt> \c wave]
55 * [<tt>--lal-eff-dist</tt>]
56 * [<tt>--usertag</tt> \c tag]
57 *
58 * [<tt>--tama-output</tt>]
59 * [<tt>--write-eff-dist</tt>]
60 *
61 * [<tt>--ilwd</tt>]</dd>
62 *
63 * <dt>Description</dt><dd>
64 * \c lalapps_inspinj
65 * generates a number of inspiral parameters suitable for using in a Monte
66 * Carlo injection to test the efficiency of a inspiral search. The various
67 * parameters (detailed below) are randomly chosen and are appropriate for a
68 * particular population of binary neutron stars whose spatial distribution
69 * includes the Milky Way and a number of extragalactic objects that are input
70 * in a datafile. The possible mass pairs for the binary neutron star com-
71 * panions are also specified in a (different) datafile.
72 *
73 * The output of this program is a list of the injected events, starting
74 * at the specified start time and ending at the specified end time. One
75 * injection with random inspiral parameters will be made every specified time
76 * step, and will be randomly placed within the specified time interval.
77 * The output is written to a file name in the standard inspiral pipeline format:
78 *
79 * \code
80 * HL-INJECTIONS_USERTAG_SEED-GPSSTART-DURATION.xml
81 * \endcode
82 *
83 * where \c USERTAG is \c tag as specfied on the command line,
84 * \c SEED is the value of the random number seed chosen and
85 * \c GPSSTART and \c DURATION describes the GPS time interval that
86 * the file covers. The file is in the standard LIGO lightweight XML format
87 * containing a \c sim_inspiral table that describes the injections.
88 * In addition, an ASCII log file called <tt>injlog.txt</tt> is also written.
89 * If a <tt>--user-tag</tt> is not specified on the command line, the
90 * \c _USERTAG part of the filename will be omitted.</dd>
91 *
92 * <dt>Options</dt><dd>
93 * <ul>
94 * <li><tt>--help</tt>: Print a help message.</li>
95 *
96 * <li><tt>--source-file</tt> \c sfile:
97 * Optional. Data file containing spatial distribution of extragalactic objects.
98 * Default is the file <tt>inspsrcs.dat</tt> provided by LALApps. If that file is
99 * empty, all signals are in the Milky Way.</li>
100 *
101 * <li><tt>--mass-file</tt> \c mfile:
102 * Optional. Data file containing mass pairs for the binary neutron star
103 * companions. Default is the file <tt>BNSMasses.dat</tt> provided by LALApps.</li>
104 *
105 * <li><tt>--gps-start-time</tt> \c tstart:
106 * Optional. Start time of the injection data to be created. Defaults to the
107 * start of S2, Feb 14 2003 16:00:00 UTC (GPS time 729273613)</li>
108 *
109 * <li><tt>--gps-end-time</tt> \c tend:
110 * Optional. End time of the injection data to be created. Defaults to the end of
111 * S2, Apr 14 2003 15:00:00 UTC (GPS time 734367613).</li>
112 *
113 * <li><tt>--time-step</tt> \c tstep:
114 * Optional. Sets the time step interval between injections. The injections will
115 * occur with an average spacing of \c tstep seconds. Defaults to
116 * \f$2630/\pi\f$.</li>
117 *
118 * <li><tt>--time-interval</tt> \c tinterval:
119 * Optional. Sets the time interval during which an injection can occur.
120 * Injections are uniformly distributed over the interval. Setting \c tstep
121 * to \f$6370\f$ and \c tinterval to 600 guarantees there will be one injection
122 * into each playground segment and they will be randomly distributed within the
123 * playground times - taken the fact that your gps start time coincides with start of a playground segment.</li>
124 *
125 * <li><tt>--seed</tt> \c seed:
126 * Optional. Seed the random number generator with the integer \c seed.
127 * Defaults to \f$1\f$.</li>
128 *
129 * <li><tt>--waveform</tt> \c wave:
130 * Optional. The string \c wave will be written into the \c waveform
131 * column of the \c sim_inspiral table output. This is used by the
132 * inspiral code to determine which type of waveforms it should inject into the
133 * data. Defaults is \c GeneratePPNtwoPN.</li>
134 *
135 * <li><tt>--lal-eff-dist</tt>:
136 * Optional. If this option is specified, the effective distance will be
137 * calculated using routines from LAL. Otherwise, the default behaviour is to
138 * use an independent method contained in inspinj.c. There is good agreement
139 * between these two methods, see below for more details.</li>
140 *
141 * <li><tt>--user-tag</tt> \c string: Optional. Set the user tag for this
142 * job to be \c string. May also be specified on the command line as
143 * <tt>-userTag</tt> for LIGO database compatibility.</li>
144 *
145 * <li><tt>--tama-output</tt>:
146 * Optional. If this option is given, \c lalapps_inspinj also produces a
147 * text output file:
148 *
149 * \code
150 * HLT-INJECTIONS_USERTAG_SEED-GPSSTART-DURATION.txt
151 * \endcode
152 *
153 * which contains the following fields:
154 *
155 * <ul>
156 * <li> geocentric end time</li>
157 * <li> Hanford end time</li>
158 * <li> Livingston end time</li>
159 * <li> TAMA end time</li>
160 * <li> total mass, \f$M_{\mathrm{TOT}}\f$</li>
161 * <li> mass ratio, \f$\eta\f$</li>
162 * <li> distance to source (in kpc)</li>
163 * <li> longitude</li>
164 * <li> latitude</li>
165 * <li> inclination</li>
166 * <li> coalescence phase</li>
167 * <li> polarization</li>
168 * <li> TAMA polarization</li>
169 * <li> end time GMST</li>
170 * </ul>
171 *
172 * In the above, all times are recorded as double precision real numbers and all
173 * angles are in radians. The TAMA polarization is calculated using
174 *
175 * \f{equation}{
176 * \tan( \psi_{T} ) = \frac{ x \cdot T_{z} }{ y \cdot T_{z} } \, .
177 * \f}
178 *
179 * Here x and y are the x,y axes of the radiation frame expressed in earth fixed
180 * coordinates \eqref{xrad}, \eqref{yrad}. \f$T_{z}\f$ is a unit vector in earth fixed
181 * coordinates which is orthogonal to the two arms of the TAMA detector
182 * \eqref{tarm}. It is given by
183 *
184 * \f{equation}{
185 * T_{z} = ( -0.6180, +0.5272, +0.5832 )
186 * \f}</li>
187 *
188 * <li><tt>--write-eff-dist</tt>: Optional. If this option is given, three extra
189 * columns are added to the TAMA output file described above. They are
190 * <ul>
191 * <li> Hanford effective distance (kpc)</li>
192 * <li> Livingston effective distance (kpc)</li>
193 * <li> TAMA effective distance (kpc)</li>
194 * </ul>
195 *
196 * These entries are added to the list immediately after TAMA end time and before
197 * total mass.</li>
198 *
199 * <li><tt>--ilwd</tt>: Optional. If this option is given,
200 * \c lalapps_inspinj also produces two ILWD-format files, injepochs.ilwd and
201 * injparams.ilwd, that contain, respectively, the GPS times suitable for
202 * inspiral injections, and the intrinsic inspiral signal parameters to be used
203 * for those injections.
204 *
205 * The file injepochs.ilwd contains a sequence of integer pairs representing
206 * the injection GPS time in seconds and residual nano-seconds. The file
207 * injparams.ilwd contains the intrinsic binary parameters for each injection,
208 * which is a sequence of eight real numbers representing (in order) (1) the
209 * total mass of the binary system (in solar masses), (2) the dimensionless
210 * reduced mass --- reduced mass per unit total mass --- in the range from 0
211 * (extreme mass ratio) to 0.25 (equal masses), (3) the distance to the system
212 * in meters, (4) the inclination of the binary system orbit to the plane of
213 * the sky in radians, (5) the coalescence phase in radians, (6) the longitude
214 * to the direction of the source in radians, (7) the latitude to the
215 * direction of the source in radians, (8) and the polar- ization angle of the
216 * source in radians.</li>
217 * </ul></dd>
218 *
219 * <dt>Example</dt><dd>
220 * \code
221 * lalapps_inspinj --seed 45\
222 * --source-file inspsrcs.dat --mass-file BNSMasses.dat
223 * \endcode</dd>
224 *
225 * <dt>Algorithm</dt><dd>
226 *
227 * The algorithm for computing the effective distance will be described in some
228 * detail below. The method is to compute both the strain due to the inspiral
229 * and the detector response in the earth fixed frame. This frame is such that
230 * the z-axis points from the earth's centre to the North Pole, the x-axis points
231 * from the centre to the intersection of the equator and the prime meridian and
232 * the y-axis is chosen to complete the orthonormal basis. The coordinates of
233 * the injection are specified by longitude (or right ascension) \f$\alpha\f$ and
234 * latitude (or declination) \f$\delta\f$. The polarization is appropriate for
235 * transferring from the radiation to earth fixed frame. These are then
236 * converted to the earth fixed frame by
237 *
238 * \f{eqnarray}{
239 * \theta &=& \frac{\pi}{2} - \delta \\
240 * \phi &=& \alpha - \textrm{gmst} \, .
241 * \f}
242 *
243 * Here, gmst is the Greenwich Mean sidereal time of the injection. The axes of
244 * the radiation frame (x,y,z) can be expressed in terms of the earth fixed
245 * coordinates as:
246 *
247 * \f{eqnarray}{
248 * x(1) &=& +( \sin( \phi ) \cos( \psi ) - \sin( \psi ) \cos( \phi )
249 * \cos( \theta ) ) \nonumber \\
250 * x(2) &=& -( \cos( \phi ) \cos( \psi ) + \sin( \psi ) \sin( \phi )
251 * \cos( \theta ) ) \nonumber \\
252 * x(3) &=& \sin( \psi ) \sin( \theta ) \label{xrad}\\
253 * y(1) &=& -( \sin( \phi ) \sin( \psi ) + \cos( \psi ) \cos( \phi )
254 * \cos( \theta ) ) \nonumber\\
255 * y(2) &=& +( \cos( \phi ) \sin( \psi ) - \cos( \psi ) \sin( \phi )
256 * \cos( \theta ) ) \nonumber \\
257 * y(3) &=& \cos( \psi ) \sin( \theta ) \label{yrad}
258 * \f}
259 *
260 * Making use of these expressions, we can express the gravitational wave strain in
261 * earth fixed coordinates as
262 *
263 * \f{equation}{\label{hij}
264 * h_{ij} = ( h^{+}(t) e^{+}_{ij} ) + (h^{\times}(t) e^{\times}_{ij})
265 * \f}
266 *
267 * where
268 *
269 * \f{equation}{
270 * e^{+}_{ij} = x_{i} * x_{j} - y_{i} * y_{j} \qquad \mathrm{and} \qquad
271 * e^{\times}_{ij} = x_{i} * y_{j} + y_{i} * x_{j}.
272 * \f}
273 *
274 * For the case of a binary inspiral signal, the two polarizations \f$h^{+}\f$
275 * and \f$h^{\times}\f$ of the gravitational wave are given by
276 *
277 * \f{eqnarray}{
278 * h^{+}(t) &=& \frac{A}{r} ( 1 + \cos^2 ( \iota ) ) * \cos( \Phi(t) ) \\
279 * h^{\times}(t) &=& \frac{A}{r} * ( 2 \cos( \iota ) ) * \sin( \Phi(t) )
280 * \f}
281 *
282 * where \f$A\f$ is a mass and frequency dependent amplitude factor, \f$r\f$ is the
283 * physical distance at which the injection is located and \f$\iota\f$ is the
284 * inclination angle.
285 *
286 * Next, we can write the detector response function as
287 *
288 * \f{equation}{
289 * d^{ij} = \left(\frac{1}{2} \right) \left( n_{x}^{i} n_{x}^{j}
290 * - n_{y}^{i} n_{y}^{j} \right) \, .
291 * \f}
292 *
293 * Here, \f$n_{x}\f$ and \f$n_{y}\f$ are unit vectors directed along the arms of the
294 * detector. Specifically, for the Hanford, Livingston, GEO, TAMA and Virgo
295 * detectors we use:
296 *
297 * \f{eqnarray}{
298 * H_{x} &=& ( -0.2239, +0.7998, +0.5569 ) \nonumber \\
299 * H_{y} &=& ( -0.9140, +0.0261, -0.4049 ) \\
300 * L_{x} &=& ( -0.9546, -0.1416, -0.2622 ) \nonumber \\
301 * L_{y} &=& ( +0.2977, -0.4879, -0.8205 ) \\
302 * G_{x} &=& ( -0.6261, -0.5522, +0.5506 ) \nonumber \\
303 * G_{y} &=& ( -0.4453, +0.8665, +0.2255 ) \\
304 * T_{x} &=& ( +0.6490, +0.7608, +0.0000 ) \nonumber \\
305 * T_{y} &=& ( -0.4437, +0.3785, -0.8123 ) \label{tarm} \\
306 * V_{x} &=& ( -0.7005, +0.2085, +0.6826 ) \nonumber \\
307 * V_{y} &=& ( -0.0538, -0.9691, +0.2408 )
308 * \f}
309 *
310 * The response of an interferometric detector with arm locations given by \f$n_{x}\f$
311 * and \f$n_{y}\f$ to an inspiralling binary system described by \eqref{hij} is
312 *
313 * \f{eqnarray}{
314 * h(t) &=& h^{+}(t) ( d^{ij} e^{+}_{ij} )
315 * + h^{\times}(t) ( d^{ij} e^{\times}_{ij} ) \nonumber \\
316 * &=&
317 * \left(\frac{A}{r}\right) \left[
318 * ( 1 + \cos^2 ( \iota ) ) F_{+} \cos( \Phi(t)) +
319 * 2 \cos( \iota ) F_{\times} \sin( \Phi(t) ) \right] \, ,
320 * \f}
321 *
322 * where we have introduced
323 *
324 * \f{equation}{
325 * F_{+} = d^{ij} e^{+}_{ij} \qquad \mathrm{and} \qquad
326 * F_{\times} = d^{ij} e^{\times}_{ij}
327 * \f}
328 *
329 * Finally, to calculate the effective distance, we note that the two contributions
330 * to \f$h(t)\f$ are \f$\pi/2\f$ radians out of phase, and hence orthogonal. Thus, we can
331 * compute the effective distance to be:
332 *
333 * \f{equation}{
334 * D_{\mathrm{eff}} = r / \left( \frac{ (1 + \cos^2(\iota))^2 }{4} F_{+}^{2} +
335 * cos^{2}(\iota) F_{\times}^{2} \right)
336 * \f}
337 *
338 * \anchor eff_dist_comparison
339 * \image html effective_distance_comparison.png "Comparison of effective distance computed by inspinj.c and LAL routines"
340 *
341 * The algorithm to calculate effective distances described above is completely
342 * contained within inspinj.c. There is an independent method of computing
343 * effective distances can also be called by inspinj. It is contained in the LAL
344 * function <tt>LALPopulateSimInspiralSiteInfo()</tt>. This function populates
345 * the site end time and effective distance for all the interferomter sites. It
346 * makes use of LAL functionality in the tools and date packages. These same
347 * functions are used when generating the injection waveform which is added to
348 * the data stream (in lalapps_inspiral). As a check that these two
349 * calculations produce the same effective distance, lalapps_inspinj was run
350 * twice, once with the <tt>--lal-eff-dist</tt> option and once without.
351 * \ref eff_dist_comparison "This figure" shows the fractional difference in effective
352 * distance between the two methods for a set of injections. We see that the
353 * distances agree within 1\
354 * occuring for the largest effective disances, i.e. close to the dead spot of
355 * the instrument. For injections which initial LIGO is sensitive to, the
356 * accuracy is few \f$\times 10^{-4}\f$. </dd>
357 *
358 * <dt>Environment</dt><dd>
359 * <ul>
360 * <li>LALAPPS_DATA_PATH: Directory to look for the default mass
361 * file <tt>BNSMasses.dat</tt> and the default source file <tt>inspsrcs.dat</tt>.</li>
362 * </ul></dd>
363 *
364 * <dt>Author</dt><dd>
365 * Jolien Creighton, Patrick Brady, Duncan Brown</dd>
366 * </dl>
367 */
368
369#include "config.h"
370
371#include <ctype.h>
372#include <lal/Date.h>
373#include <lal/LALgetopt.h>
374#include <lal/LIGOLwXML.h>
375#include <lal/LIGOLwXMLRead.h>
376#include <lal/LIGOMetadataTables.h>
377#include <lal/LIGOMetadataUtils.h>
378#include <lal/LIGOMetadataInspiralUtils.h>
379#include <lal/Random.h>
380#include <lal/AVFactories.h>
381#include <lal/InspiralInjectionParams.h>
382#include <lal/LALDetectors.h>
383#include <lal/LALSimulation.h>
384#include <lal/RingUtils.h>
385#include <LALAppsVCSInfo.h>
386#include <lal/LALDatatypes.h>
387#include <lal/FrequencySeries.h>
388#include "inspiral.h"
389#include "LALExtTriggerTableFromLIGOLw.h"
390
391#define CVS_REVISION "$Revision$"
392#define CVS_ID_STRING "$Id$"
393#define CVS_SOURCE "$Source$"
394#define CVS_DATE "$Date$"
395#define CVS_NAME_STRING "$Name$"
396#define PROGRAM_NAME "inspinj"
397
398#define ADD_PROCESS_PARAM( pptype, format, ppvalue ) \
399 this_proc_param = this_proc_param->next = (ProcessParamsTable *) \
400calloc( 1, sizeof(ProcessParamsTable) ); \
401snprintf( this_proc_param->program, LIGOMETA_PROGRAM_MAX, "%s", \
402 PROGRAM_NAME ); \
403snprintf( this_proc_param->param, LIGOMETA_PARAM_MAX, "--%s", \
404 long_options[option_index].name ); \
405snprintf( this_proc_param->type, LIGOMETA_TYPE_MAX, "%s", pptype ); \
406snprintf( this_proc_param->value, LIGOMETA_VALUE_MAX, format, ppvalue );
407
408#ifdef __GNUC__
409#define UNUSED __attribute__ ((unused))
410#else
411#define UNUSED
412#endif
413
414#define AXIS_MAX 12
415
416/*
417 * *********************************
418 * Definition of the prototypes
419 * *********************************
420 */
421extern int vrbflg;
422ProcessParamsTable *next_process_param( const char *name, const char *type,
423 const char *fmt, ... );
424void read_mass_data( char *filename );
425void read_time_data( char *filename );
426void read_nr_data( char* filename );
427void read_source_data( char* filename );
428void sourceComplete(void);
429void drawFromSource( REAL8 *rightAscension,
430 REAL8 *declination,
431 REAL8 *distance,
432 CHAR name[LIGOMETA_SOURCE_MAX] );
433void read_IPN_grid_from_file( char *fname );
434void drawFromIPNsim( REAL8 *rightAscension,
435 REAL8 *declination );
436void drawLocationFromExttrig( SimInspiralTable* table );
437void drawMassFromSource( SimInspiralTable* table );
438void drawMassSpinFromNR( SimInspiralTable* table );
439void drawMassSpinFromNRNinja2( SimInspiralTable* table );
440
441void adjust_snr(SimInspiralTable *inj, REAL8 target_snr, const char *ifos);
442void adjust_snr_real8(SimInspiralTable *inj, REAL8 target_snr, const char *ifos);
443REAL8 network_snr(const char *ifos, SimInspiralTable *inj);
444REAL8 snr_in_ifo(const char *ifo, SimInspiralTable *inj);
445REAL8 network_snr_real8(const char *ifos, SimInspiralTable *inj);
446REAL8 snr_in_ifo_real8(const char *ifo, SimInspiralTable *inj);
447
448REAL8 snr_in_psd_real8(const char *ifo, REAL8FrequencySeries *psd, REAL8 start_freq, SimInspiralTable *inj);
449REAL8 network_snr_with_psds_real8(int num_ifos, const char **ifo_list, REAL8FrequencySeries **psds, REAL8 *start_freqs, SimInspiralTable *inj);
450void adjust_snr_with_psds_real8(SimInspiralTable *inj, REAL8 target_snr, int num_ifos, const char **ifo_list, REAL8FrequencySeries **psds, REAL8 *start_freqs);
451
452REAL8 probability_redshift(REAL8 rshift);
453REAL8 luminosity_distance(REAL8 rshift);
454REAL8 mean_time_step_sfr(REAL8 zmax, REAL8 rate_local);
455REAL8 drawRedshift(REAL8 zmin, REAL8 zmax, REAL8 pzmax);
456REAL8 redshift_mass(REAL8 mass, REAL8 z);
457static void scale_lalsim_distance(SimInspiralTable *inj,char ** IFOnames, REAL8FrequencySeries **psds,REAL8 *start_freqs,LoudnessDistribution dDistr);
458static REAL8 draw_uniform_snr(REAL8 snrmin,REAL8 snrmax);
459static REAL8 draw_log10_snr(REAL8 snrmin,REAL8 snrmax);
460static REAL8 draw_volume_snr(REAL8 snrmin,REAL8 snrmax);
461/*
462 * *************************************
463 * Defining of the used global variables
464 * *************************************
465 */
466
467lalinspiral_time_distribution tDistr;
468LoudnessDistribution dDistr;
469SkyLocationDistribution lDistr;
470MassDistribution mDistr;
471InclDistribution iDistr;
472SpinDistribution spinDistr = uniformSpinDist;
473
474SimInspiralTable *simTable;
475SimRingdownTable *simRingTable;
476
477char *massFileName = NULL;
478char *nrFileName = NULL;
479char *sourceFileName = NULL;
480char *outputFileName = NULL;
481char *injtimesFileName = NULL;
482char *exttrigFileName = NULL;
483char *IPNSkyPositionsFile = NULL;
484
485INT4 outCompress = 0;
486INT4 ninjaMass = 0;
487INT4 real8Ninja2 = 0;
488INT4 ninjaSNR = 0;
489INT4 haveLoudness= 0;
490
491REAL4 mwLuminosity = -1;
492REAL4 minD = -1;
493REAL4 maxD = -1;
494REAL4 minZ = -1;
495REAL4 maxZ = -1;
496REAL4 localRate = -1.0;
497REAL4 minSNR = -1;
498REAL4 maxSNR = -1;
499char *ifos = NULL;
500
501REAL4 minMass1 = -1;
502REAL4 maxMass1 = -1;
503REAL4 minMass2 = -1;
504REAL4 maxMass2 = -1;
505REAL4 minMtotal = -1;
506REAL4 maxMtotal = -1;
507REAL4 meanMass1 = -1.0;
508REAL4 meanMass2 = -1.0;
509REAL4 massStdev1= -1.0;
510REAL4 massStdev2= -1.0;
511REAL4 minMassRatio=-1.0;
512REAL4 maxMassRatio=-1.0;
513
514REAL4 inclStd=-1.0;
515REAL4 fixed_inc=-1.0;
516REAL4 max_inc=LAL_PI/2.0;
517int coaPhaseFixed = 0;
518REAL4 fixedCoaPhase = 0;
519REAL4 psi=-1.0;
520REAL4 longitude=181.0;
521REAL4 latitude=91.0;
522REAL4 epsAngle=1e-7;
523int spinInjections=-1;
524int spinAligned=-1;
525REAL4 minSpin1=-1.0;
526REAL4 maxSpin1=-1.0;
527REAL4 meanSpin1=0.0;
528REAL4 Spin1Std=0.0;
529REAL4 minSpin2=-1.0;
530REAL4 maxSpin2=-1.0;
531REAL4 meanSpin2=0.0;
532REAL4 Spin2Std=0.0;
533REAL4 minKappa1=-1.0;
534REAL4 maxKappa1=1.0;
535REAL4 minabsKappa1=0.0;
536REAL4 maxabsKappa1=1.0;
537REAL4 fixedMass1=-1.0;
538REAL4 fixedMass2=-1.0;
539INT4 pntMass1=1;
540INT4 pntMass2=1;
541REAL4 deltaMass1=-1;
542REAL4 deltaMass2=-1;
543INT4 bandPassInj = 0;
544LALSimInspiralApplyTaper taperInj = LAL_SIM_INSPIRAL_TAPER_NONE;
545AlignmentType alignInj = notAligned;
546REAL8 redshift;
547
548REAL8 single_IFO_SNR_threshold=0.0;
549char ** ifonames=NULL;
550int numifos=0;
551
552static LALStatus status;
553static RandomParams* randParams=NULL;
554INT4 numExtTriggers = 0;
555ExtTriggerTable *exttrigHead = NULL;
556
557int num_source;
558int numSkyPoints;
559int galaxynum;
560struct {
561 char name[LIGOMETA_SOURCE_MAX];
562 REAL8 ra;
563 REAL8 dec;
564 REAL8 dist;
565 REAL8 lum;
566 REAL8 fudge;
567} *source_data, *old_source_data,*temparray, *skyPoints;
568
569char MW_name[LIGOMETA_SOURCE_MAX] = "MW";
570REAL8* fracVec =NULL;
571REAL8* ratioVec = NULL;
572REAL8 norm=0;
573
574int num_mass;
575struct {
576 REAL8 mass1;
577 REAL8 mass2;
578} *mass_data;
579
580int n_times;
581LIGOTimeGPS* inj_times;
582
583struct FakeGalaxy{
584char name[LIGOMETA_SOURCE_MAX];
585REAL8 ra;
586REAL8 dec;
587REAL8 lum;
588REAL8 dist;
589REAL8 fudge;
590struct FakeGalaxy *next; };
591int srcComplete = 0;
592int makeCatalog = 0;
593REAL8 srcCompleteDist;
594
595int num_nr = 0;
596int i = 0;
597SimInspiralTable **nrSimArray = NULL;
598
599/*
600 * *****************************************************
601 * Functions implementing SFR distribution over redshift
602 * *****************************************************
603 */
604
605REAL8 probability_redshift(REAL8 rshift)
606{
607 REAL8 pz;
608
609 pz = -0.000429072589677+(rshift*(-0.036349728568888+(rshift*(0.860892111762314
610 +(rshift*(-0.740935488674010+rshift*(0.265848831356864+rshift*(-0.050041573542298
611 +rshift*(0.005184554232421+rshift*(-0.000281450045300+rshift*0.000006400690921))))))))));
612
613 return pz;
614}
615
616REAL8 luminosity_distance(REAL8 rshift)
617{
618 REAL8 dL;
619
620 dL = -2.89287707063171+(rshift*(4324.33492012756+(rshift*(3249.74193862773
621 +(rshift*(-1246.66339928289+rshift*(335.354613407693+rshift*(-56.1194965448065
622 +rshift*(5.20261234121263+rshift*(-0.203151569744028))))))))));
623
624 return dL;
625}
626
627REAL8 mean_time_step_sfr(REAL8 zmax, REAL8 rate_local)
628{
629 REAL8 logzmax,loglambda,step;
630
631 logzmax=log10(zmax);
632 loglambda = -0.039563*pow(logzmax,6.)-0.15282*pow(logzmax,5.)-0.017596*pow(logzmax,4.)
633 + 0.67193*pow(logzmax,3.)+1.1347*pow(logzmax,2.)-2.3543*logzmax+ 2.0228;
634 step=pow(10.,loglambda)/rate_local;
635
636 return step;
637}
638
639REAL8 drawRedshift(REAL8 zmin, REAL8 zmax, REAL8 pzmax)
640{
641 REAL8 test,z,p;
642 do
643 {
644 test = pzmax * XLALUniformDeviate(randParams);
645 z = (zmax-zmin) * XLALUniformDeviate(randParams)+zmin;
646 p = probability_redshift(z);
647 }
648 while (test>p);
649
650 return z;
651}
652
653REAL8 redshift_mass(REAL8 mass, REAL8 z)
654{
655 REAL8 mz;
656 mz = mass * (1.+z);
657
658 return mz;
659}
660
661
662/*************************************************************
663 * Routines that calculate/adjust SNRs for REAL4 NINJA-1
664 * injections. In principle these are obsolete and could be
665 * deleted.
666 *************************************************************/
667REAL8 snr_in_ifo(const char *ifo, SimInspiralTable *inj)
668{
669 REAL8 this_snr;
670 REAL4TimeVectorSeries *tempStrain=NULL;
671
672 AddNumRelStrainModes( &status, &tempStrain, inj);
673
674 this_snr = calculate_ligo_snr_from_strain( tempStrain, inj, ifo);
675
676 XLALDestroyREAL4VectorSequence (tempStrain->data);
677 tempStrain->data = NULL;
678 LALFree(tempStrain);
679 tempStrain = NULL;
680
681 return this_snr;
682}
683
684REAL8 network_snr(const char *ifo_list, SimInspiralTable *inj)
685{
686 char *tmp;
687 char *ifo;
688 REAL8 snr_total = 0.0;
689 REAL8 this_snr;
690
691 tmp = LALCalloc(1, strlen(ifos) + 1);
692 strcpy(tmp, ifo_list);
693
694 ifo = strtok (tmp,",");
695 while (ifo != NULL)
696 {
697 this_snr = snr_in_ifo(ifo, inj);
698 snr_total += this_snr * this_snr;
699 ifo = strtok (NULL, ",");
700 }
701
702 LALFree(tmp);
703
704 return sqrt(snr_total);
705}
706
707void adjust_snr(SimInspiralTable *inj, REAL8 target_snr, const char *ifo_list)
708{
709 /* Vars for calculating SNRs */
710 REAL8 this_snr;
711 REAL8 UNUSED low_snr, UNUSED high_snr;
712 REAL8 low_dist,high_dist;
713
714 this_snr = network_snr(ifo_list, inj);
715
716 if (this_snr > target_snr)
717 {
718 high_snr = this_snr;
719 high_dist = inj->distance;
720
721 while (this_snr > target_snr)
722 {
723 inj-> distance = inj->distance * 3.0;
724 this_snr = network_snr(ifo_list, inj);
725 }
726 low_snr = this_snr;
727 low_dist = inj->distance;
728 } else {
729 low_snr = this_snr;
730 low_dist = inj->distance;
731
732 while (this_snr < target_snr)
733 {
734 inj->distance = (inj->distance) / 3.0;
735 this_snr = network_snr(ifo_list, inj);
736 }
737 high_snr = this_snr;
738 high_dist = inj->distance;
739 }
740
741 while ( fabs(target_snr - this_snr) > 1.0 )
742 {
743 inj->distance = (high_dist + low_dist) / 2.0;
744 this_snr = network_snr(ifo_list, inj);
745
746 if (this_snr > target_snr)
747 {
748 high_snr = this_snr;
749 high_dist = inj->distance;
750 } else {
751 low_snr = this_snr;
752 low_dist = inj->distance;
753 }
754 }
755}
756
757
758/*************************************************************
759 * Routines that calculate/adjust SNRs for REAL8 NINJA-2
760 * injections, using the initial LIGO/Virgo noise curves from
761 * the noisemodels package. In principle these could be replaced
762 * with the next group, which will default to these noise curves
763 * when alternatives are not provided.
764 *************************************************************/
765REAL8 snr_in_ifo_real8(const char *ifo, SimInspiralTable *inj)
766{
767 REAL8 this_snr;
768 REAL8TimeSeries *strain = NULL;
769
770 strain = XLALNRInjectionStrain(ifo, inj);
771 this_snr = calculate_ligo_snr_from_strain_real8(strain, ifo);
772
773 XLALDestroyREAL8TimeSeries (strain);
774
775 return this_snr;
776}
777
778
779REAL8 network_snr_real8(const char *ifo_list, SimInspiralTable *inj)
780{
781 char *tmp;
782 char *ifo;
783
784 REAL8 snr_total = 0.0;
785 REAL8 this_snr;
786
787 tmp = LALCalloc(1, strlen(ifos) + 1);
788 strcpy(tmp, ifo_list);
789 ifo = strtok (tmp,",");
790
791 while (ifo != NULL)
792 {
793 this_snr = snr_in_ifo_real8(ifo, inj);
794 snr_total += this_snr * this_snr;
795 ifo = strtok (NULL, ",");
796 }
797
798 LALFree(tmp);
799
800 return sqrt(snr_total);
801}
802
803
804void adjust_snr_real8(
805 SimInspiralTable *inj,
806 REAL8 target_snr,
807 const char *ifo_list)
808{
809 /* Vars for calculating SNRs */
810 REAL8 this_snr;
811 REAL8 UNUSED low_snr, UNUSED high_snr;
812 REAL8 low_dist,high_dist;
813
814 this_snr = network_snr_real8(ifo_list, inj);
815
816 if (this_snr > target_snr)
817 {
818 high_snr = this_snr;
819 high_dist = inj->distance;
820
821 while (this_snr > target_snr)
822 {
823 inj-> distance = inj->distance * 3.0;
824 this_snr = network_snr_real8(ifo_list, inj);
825 }
826 low_snr = this_snr;
827 low_dist = inj->distance;
828 }
829 else
830 {
831 low_snr = this_snr;
832 low_dist = inj->distance;
833
834 while (this_snr < target_snr)
835 {
836 inj->distance = (inj->distance) / 3.0;
837 this_snr = network_snr_real8(ifo_list, inj);
838 }
839 high_snr = this_snr;
840 high_dist = inj->distance;
841 }
842
843 while ( fabs(target_snr - this_snr) > 1.0 )
844 {
845 inj->distance = (high_dist + low_dist) / 2.0;
846 this_snr = network_snr_real8(ifo_list, inj);
847
848 if (this_snr > target_snr)
849 {
850 high_snr = this_snr;
851 high_dist = inj->distance;
852 }
853 else
854 {
855 low_snr = this_snr;
856 low_dist = inj->distance;
857 }
858 }
859}
860
861
862/*************************************************************
863 * Routines that calculate/adjust SNRs for REAL8 NINJA-2
864 * injections, using arbitrary LIGO/Virgo noise curves given
865 * in files.
866 *************************************************************/
867REAL8 snr_in_psd_real8(
868 const char *ifo,
869 REAL8FrequencySeries *psd,
870 REAL8 start_freq,
871 SimInspiralTable *inj)
872{
873 REAL8 this_snr;
874 REAL8TimeSeries *strain = NULL;
875
876 strain = XLALNRInjectionStrain(ifo, inj);
877 this_snr = calculate_snr_from_strain_and_psd_real8(strain, psd, start_freq, ifo);
878
879 XLALDestroyREAL8TimeSeries (strain);
880
881 return this_snr;
882}
883
884REAL8 network_snr_with_psds_real8(
885 int num_ifos,
886 const char **ifo_list,
887 REAL8FrequencySeries **psds,
888 REAL8 *start_freqs,
889 SimInspiralTable *inj)
890{
891 REAL8 snr_total = 0.0;
892 REAL8 this_snr;
893
894 for (i=0; i< num_ifos; i++)
895 {
896 this_snr = snr_in_psd_real8(ifo_list[i], psds[i], start_freqs[i], inj);
897 snr_total += this_snr * this_snr;
898 }
899
900 return sqrt(snr_total);
901}
902
903void adjust_snr_with_psds_real8(
904 SimInspiralTable *inj,
905 REAL8 target_snr,
906 int num_ifos,
907 const char **ifo_list,
908 REAL8FrequencySeries **psds,
909 REAL8 *start_freqs)
910{
911 /* Vars for calculating SNRs */
912 REAL8 this_snr;
913
914 this_snr = network_snr_with_psds_real8(num_ifos, ifo_list, psds, start_freqs, inj);
915 inj->distance = inj->distance * (this_snr/target_snr);
916}
917
918
919/*
920 *
921 * code to step forward in the process table
922 *
923 */
924ProcessParamsTable *next_process_param( const char *name, const char *type,
925 const char *fmt, ... )
926{
927 ProcessParamsTable *pp;
928 va_list ap;
929 pp = calloc( 1, sizeof( *pp ) );
930 if ( ! pp )
931 {
932 perror( "next_process_param" );
933 exit( 1 );
934 }
935 strncpy( pp->program, PROGRAM_NAME, LIGOMETA_PROGRAM_MAX );
936 snprintf( pp->param, LIGOMETA_PARAM_MAX, "--%s", name );
937 strncpy( pp->type, type, LIGOMETA_TYPE_MAX - 1 );
938 va_start( ap, fmt );
939 vsnprintf( pp->value, LIGOMETA_VALUE_MAX, fmt, ap );
940 va_end( ap );
941 return pp;
942}
943
944/*
945 *
946 * print-out of the usage
947 *
948 */
949static void print_usage(char *program)
950{
951 fprintf(stderr,
952 "%s [options]\n"\
953 "The following options are recognized. Options not surrounded in []\n"\
954 "are required. Defaults are shown in brackets\n", program );
955 fprintf(stderr,
956 " [--help ] display this message\n"\
957 " [--verbose] print progress information\n"\
958 " [--user-tag] usertag set the usertag \n"\
959 " [--output ] name overwrite the standard file naming convention\n"\
960 " [--write-compress] write a compressed xml file\n\n");\
961 fprintf(stderr,
962 "Waveform details:\n"\
963 " [--seed] randomSeed seed for random number generator (default : 1)\n"\
964 " --f-lower freq lower cut-off frequency.\n"\
965 " --waveform wfm set waveform type to wfm\n"\
966 " --amp-order set PN order in amplitude\n\n");
967 fprintf(stderr,
968 "Time distribution information:\n"\
969 " --gps-start-time start GPS start time for injections\n"\
970 " --gps-end-time end GPS end time for injections\n"\
971 " --ipn-gps-time IPNtime GPS end time for IPN trigger\n"\
972 " --t-distr timeDist set the time step distribution of injections\n"\
973 " fixed: fixed time step\n"\
974 " uniform: uniform distribution\n"\
975 " exponential: exponential distribution for Poisson process\n"\
976 " [--time-step] step space injections by average of step seconds\n"\
977 " (suggestion : 2630 / pi seconds)\n"\
978 " [--time-interval] int distribute injections in an interval, int s\n"\
979 " (default : 0 seconds)\n\n");
980 fprintf(stderr,
981 "Source distribution information:\n"\
982 " --l-distr locDist set the source location distribution,\n"\
983 " locDist must be one of:\n"\
984 " source: use locations from source-file\n"\
985 " exttrig: use external trigger file\n"\
986 " random: uses random locations\n"\
987 " fixed: set fixed location\n"\
988 " ipn: random locations from IPN skypoints\n"\
989 " [--longitude] longitude read longitude if fixed value (degrees)\n"
990 " [--latitude] latitude read latitude if fixed value (degrees)\n"
991 " [--d-distr] distDist use a distribution over physical distance\n"\
992 " source: take distance from galaxy source file\n"\
993 " uniform: uniform distribution in distance\n"\
994 " distancesquared: uniform distribution in distance^2\n"\
995 " log10: uniform distribution in log10(d) \n"\
996 " volume: uniform distribution in volume\n"\
997 " sfr: distribution derived from the SFR\n"\
998 " [--min-distance] DMIN set the minimum (chirp) distance to DMIN kpc\n"\
999 " [--max-distance] DMAX set the maximum (chirp) distance to DMAX kpc\n"\
1000 " min/max distance required if d-distr not 'source'\n"\
1001 " [--source-file] sources read source parameters from sources\n"\
1002 " requires enable/disable milkyway\n"\
1003 " [--sourcecomplete] distance \n"
1004 " complete galaxy catalog out to distance (kPc)\n"\
1005 " [--make-catalog] create a text file of the completed galaxy catalog\n"\
1006 " [--enable-milkyway] lum enables MW injections, set MW luminosity\n"\
1007 " [--disable-milkyway] disables Milky Way injections\n"\
1008 " [--dchirp-distr] use a distribution over chirp distance\n"\
1009 " (normalized to a 1.4,1.4 Msun binary)\n"\
1010 " [--z-distr] use a distribution over redshift\n"\
1011 " currently only 'sfr' is supported\n"\
1012 " [--local-rate] rho set the local coalescence rate for --z-distr sfr\n"\
1013 " (suggestion: 1 per Mpc^3 per Myr)\n"\
1014 " [--min-z] set the minimum redshift: at least 0.2 for sfr\n"\
1015 " [--max-z] set the maximum redshift: at most 1.0 for sfr\n"\
1016 " [--snr-distr] use a distribution over expected (optimal) network SNR\n"\
1017 " uniform: uniform in SNR, log10: uniform in log10(SNR)\n"\
1018 " volume: uniform in 1/SNR^3\n"\
1019 " ( Setting max-snr == min-snr will allow you to choose a fixed SNR )\n"\
1020 " [--ninja-snr] use a NINJA waveform SNR calculation (if not set, use LALSimulation)\n"\
1021 " [--min-snr] SMIN set the minimum network snr\n"\
1022 " [--max-snr] SMAX set the maximum network snr\n"\
1023 " [--min-coinc-snr] sm Set the minimum SNR in two IFOs. Neglected if a single IFO is used\n"\
1024 " [--ligo-psd] filename Ascii, tab-separated file of frequency, value pairs to use for LIGO PSD in snr computation\n"\
1025 " [--ligo-fake-psd] PSD LALsimulation PSD fit to use instead of a file. Allowed values: LALLIGO, LALAdLIGO\n"\
1026 " [--ligo-start-freq] freq Frequency in Hz to use for LIGO snr computation\n"\
1027 " [--virgo-psd] filename Ascii, tab-separated file of frequency, value pairs to use for Virgo PSD in snr computation\n"\
1028 " [--virgo-fake-psd] PSD LALsimulation PSD fit to use instead of a file. Allowed values: LALVirgo, LALAdVirgo\n"\
1029 " [--virgo-start-freq] freq Frequency in Hz to use for Virgo snr computation\n"\
1030 " [--ifos] ifos Comma-separated list of ifos to include in network SNR\n\n"\
1031 " --i-distr INCDIST set the inclination distribution, must be either\n"\
1032 " uniform: distribute uniformly over arccos(i)\n"\
1033 " gaussian: gaussian distributed in (i)\n"\
1034 " fixed: no distribution, fixed values of (i)\n"\
1035 " [--polarization] psi set the polarization angle for all injections (degrees)\n"\
1036 " [--incl-std] inclStd std dev for gaussian inclination dist\n"\
1037 " [--fixed-inc] fixed_inc value for the fixed inclination angle (in degrees) if '--i-distr fixed' is chosen.\n"\
1038 " [--max-inc] max_inc value for the maximum inclination angle (in degrees) if '--i-distr uniform' is chosen. \n"\
1039 " [--coa-phase-distr] cDist set the coalescence phase distribution,\n"\
1040 " cDist must be one of:\n"\
1041 " uniform: use random, uniformly distributed coalescence phase [default]\n"\
1042 " fixed: set fixed coalescence phase\n"\
1043 " [--fixed-coa-phase] phase set the coalescence phase (in degrees) for all injections if --coa-phase-distr=fixed\n"\
1044 " [--ipn-file] ipnskypoints read IPN sky points from file\n"\
1045 " [--exttrig-file] exttrig XML file containing external trigger\n\n");
1046 fprintf(stderr,
1047 "Mass distribution information:\n"\
1048 " --m-distr massDist set the mass distribution of injections\n"\
1049 " must be one of:\n"\
1050 " source: using file containing list of mass pairs\n"\
1051 " nrwaves: using xml file with list of NR waveforms\n"\
1052 " (requires setting max/min total masses)\n"\
1053 " totalMass: uniform distribution in total mass\n"\
1054 " componentMass: uniform in m1 and m2\n"\
1055 " gaussian: gaussian mass distribution\n"\
1056 " log: log distribution in component mass\n"\
1057 " totalMassRatio: uniform distribution in total mass and\n"\
1058 " mass ratio m1 /m2\n"\
1059 " logTotalMassUniformMassRatio: log distribution in total mass\n"\
1060 " and uniform in mass ratio\n"\
1061 " totalMassFraction: uniform distribution in total mass and\n"\
1062 " in m1 /(m1+m2)\n"\
1063 " m1m2SquareGrid: component masses on a square grid\n"\
1064 " fixMasses: fix m1 and m2 to specific values\n"\
1065 " [--ninja2-mass] use the NINJA 2 mass-selection algorithm\n"\
1066 " [--real8-ninja2] when distributing by SNR for NINJA2, assume frames are REAL8\n"\
1067 " [--mass-file] mFile read population mass parameters from mFile\n"\
1068 " [--nr-file] nrFile read mass/spin parameters from xml nrFile\n"\
1069 " [--min-mass1] m1min set the minimum component mass to m1min\n"\
1070 " [--max-mass1] m1max set the maximum component mass to m1max\n"\
1071 " [--min-mass2] m2min set the min component mass2 to m2min\n"\
1072 " [--max-mass2] m2max set the max component mass2 to m2max\n"\
1073 " [--min-mtotal] minTotal sets the minimum total mass to minTotal\n"\
1074 " [--max-mtotal] maxTotal sets the maximum total mass to maxTotal\n"\
1075 " [--fixed-mass1] fixMass1 set mass1 to fixMass1\n"\
1076 " [--fixed-mass2] fixMass2 set mass2 to fixMass2\n"\
1077 " [--mean-mass1] m1mean set the mean value for mass1\n"\
1078 " [--stdev-mass1] m1std set the standard deviation for mass1\n"\
1079 " [--mean-mass2] m2mean set the mean value for mass2\n"\
1080 " [--stdev-mass2] m2std set the standard deviation for mass2\n"\
1081 " [--min-mratio] minr set the minimum mass ratio\n"\
1082 " [--max-mratio] maxr set the maximum mass ratio\n"\
1083 " [--mass1-points] m1pnt set the number of grid points in the m1 direction if '--m-distr=m1m2SquareGrid'\n"\
1084 " [--mass2-points] m2pnt set the number of grid points in the m2 direction if '--m-distr=m1m2SquareGrid'\n\n");
1085 fprintf(stderr,
1086 "Spin distribution information:\n"\
1087 " --disable-spin disables spinning injections\n"\
1088 " --enable-spin enables spinning injections\n"\
1089 " One of these is required.\n"\
1090 " [--spin-gaussian] enable gaussian spin distribution\n"\
1091 " --aligned enforces the spins to be along the direction\n"\
1092 " of orbital angular momentum. Spin z-components are the only non-vanishing (unless '--axis-choice view' convention is chosen)\n"\
1093 " [--axis-choice] choice frame axis choice: 'angmomentum' (default) or 'view' to define convention for spin aligned case\n"\
1094 " [--min-spin1] spin1min Set the minimum spin1 to spin1min (0.0)\n"\
1095 " [--max-spin1] spin1max Set the maximum spin1 to spin1max (0.0)\n"\
1096 " [--mean-spin1] spin1mean Set the mean for |spin1| distribution\n"\
1097 " [--stdev-spin1] spin1std Set the standard deviation for |spin1|\n"\
1098 " [--min-spin2] spin2min Set the minimum spin2 to spin2min (0.0)\n"\
1099 " [--max-spin2] spin2max Set the maximum spin2 to spin2max (0.0)\n"\
1100 " [--mean-spin2] spin2mean Set the mean for |spin2| distribution\n"\
1101 " [--stdev-spin2] spin2std Set the standard deviation for |spin2|\n"\
1102 " [--min-kappa1] kappa1min Set the minimum cos(S1.L_N) to kappa1min (-1.0)\n"\
1103 " [--max-kappa1] kappa1max Set the maximum cos(S1.L_N) to kappa1max (1.0)\n"\
1104 " [--min-abskappa1] abskappa1min \n"\
1105 " Set the minimum absolute value of cos(S1.L_N)\n"\
1106 " to abskappa1min (0.0)\n"\
1107 " [--max-abskappa1] abskappa1max \n"\
1108 " Set the maximum absolute value of cos(S1.L_N) \n"\
1109 " to abskappa1max (1.0)\n\n");
1110 fprintf(stderr,
1111 "Tapering the injection waveform:\n"\
1112 " [--taper-injection] OPT Taper the inspiral template using option OPT\n"\
1113 " (start|end|startend) \n"\
1114 " [--band-pass-injection] sets the tapering method of the injected waveform\n\n");
1115}
1116
1117
1118/*
1119 *
1120 * functions to read source masses
1121 *
1122 */
1123
1124 void
1125read_mass_data( char* filename )
1126{
1127 char line[256];
1128 FILE *fp;
1129 int n = 0;
1130
1131 fp=fopen( filename, "r" );
1132 if ( ! fp )
1133 {
1134 perror( "read_mass_data" );
1135 fprintf( stderr,
1136 "Error while trying to open file %s\n",
1137 filename );
1138 exit( 1 );
1139 }
1140
1141 /* count the number of lines in the file */
1142 num_mass=0;
1143 while ( fgets( line, sizeof( line ), fp ) )
1144 ++num_mass;
1145
1146 /* alloc space for the data */
1147 mass_data = LALCalloc( num_mass, sizeof(*mass_data) );
1148 if ( !mass_data )
1149 {
1150 fprintf( stderr, "Allocation error for mass_data\n" );
1151 exit( 1 );
1152 }
1153
1154 /* 'rewind' the file */
1155 rewind( fp );
1156
1157 /* read the file finally */
1158 while ( fgets( line, sizeof( line ), fp ) )
1159 {
1160 sscanf( line, "%le %le", &mass_data[n].mass1, &mass_data[n].mass2 );
1161 n++;
1162 }
1163
1164 /* close the file */
1165 fclose( fp );
1166}
1167
1168void read_time_data( char* filename)
1169{
1170 char line[256];
1171 FILE *fp;
1172 int n = 0;
1173 INT4 this_time = 0;
1174
1175 fp=fopen( filename, "r" );
1176 if ( ! fp )
1177 {
1178 perror( "read_time_data" );
1179 fprintf( stderr,
1180 "Error while trying to open file %s\n",
1181 filename );
1182 exit( 1 );
1183 }
1184
1185 /* count the number of lines in the file */
1186 n_times=0;
1187 while ( fgets( line, sizeof( line ), fp ) )
1188 ++n_times;
1189
1190 /* alloc space for the data */
1191 inj_times = LALCalloc( n_times, sizeof(*inj_times) );
1192 if ( !inj_times )
1193 {
1194 fprintf( stderr, "Allocation error for inj_times\n" );
1195 exit( 1 );
1196 }
1197
1198 /* 'rewind' the file */
1199 rewind( fp );
1200
1201 /* read the file finally */
1202 while ( fgets( line, sizeof( line ), fp ) )
1203 {
1204 sscanf( line, "%d", &this_time);
1205 if ( this_time < 441417609 )
1206 {
1207 fprintf( stderr, "invalid injection time %d:\n"
1208 "GPS start time is prior to "
1209 "Jan 01, 1994 00:00:00 UTC:\n"
1210 "(%d specified)\n",
1211 n, this_time );
1212 exit( 1 );
1213 }
1214 inj_times[n].gpsSeconds = this_time;
1215 inj_times[n].gpsNanoSeconds = 0;
1216 // printf("%d Time: %d\t%d\n", n, inj_times[n].gpsSeconds, inj_times[n].gpsNanoSeconds);
1217 n++;
1218 }
1219
1220 /* close the file */
1221 fclose( fp );
1222}
1223
1224 void
1225read_nr_data( char* filename )
1226{
1227 SimInspiralTable *nrSimHead = NULL;
1228 SimInspiralTable *thisEvent= NULL;
1229 INT4 j = 0;
1230
1231 nrSimHead = XLALSimInspiralTableFromLIGOLw( filename );
1232 if ( !nrSimHead )
1233 {
1234 fprintf( stderr, "error: unable to read sim_inspiral table from %s\n",
1235 filename );
1236 exit( 1 );
1237 }
1238 for(num_nr = 0, thisEvent = nrSimHead; thisEvent; num_nr++, thisEvent = thisEvent->next);
1239 if ( !num_nr )
1240 {
1241 fprintf( stderr, "error: zero events in sim_inspiral table from %s\n",
1242 filename );
1243 }
1244
1245 /* allocate an array of pointers */
1246 nrSimArray = (SimInspiralTable ** )
1247 LALCalloc( num_nr, sizeof(SimInspiralTable *) );
1248
1249 if ( !nrSimArray )
1250 {
1251 fprintf( stderr, "Allocation error for nr simulations\n" );
1252 exit( 1 );
1253 }
1254
1255 for( j = 0, thisEvent=nrSimHead; j < num_nr;
1256 ++j, thisEvent = thisEvent->next )
1257 {
1258 nrSimArray[j] = thisEvent;
1259 if (j > 0)
1260 {
1261 nrSimArray[j-1]->next = NULL;
1262 }
1263 }
1264}
1265
1266
1267/*
1268 *
1269 * functions to read source distribution
1270 *
1271 */
1272
1273 void
1274read_source_data( char* filename )
1275{
1276 char line[256];
1277 FILE *fp;
1278 int j, k;
1279
1280 fp = fopen (filename, "r" );
1281 if ( ! fp )
1282 {
1283 perror( "read_source_data" );
1284 fprintf( stderr, "Could not find file %s\n", filename );
1285 exit( 1 );
1286 }
1287
1288 /* count the number of entries in this file */
1289 num_source = 0;
1290 while ( fgets( line, sizeof( line ), fp ) )
1291 if ( line[0] == '#' )
1292 continue;
1293 else
1294 ++num_source;
1295
1296 /* rewind the file */
1297 rewind( fp );
1298
1299 /* allocate space */
1300 source_data = LALCalloc( num_source, sizeof( *source_data ) );
1301 if ( ! source_data )
1302 {
1303 fprintf( stderr, "Allocation error for source_data\n" );
1304 exit( 1 );
1305 }
1306
1307 j = 0;
1308 while ( fgets( line, sizeof( line ), fp ) )
1309 if ( line[0] == '#' )
1310 continue;
1311 else
1312 {
1313 char ra_sgn, dec_sgn;
1314 REAL8 ra_h, ra_m, dec_d, dec_m;
1315 int c;
1316
1317 c = sscanf( line, "%s %c%le:%le %c%le:%le %le %le %le",
1318 source_data[j].name, &ra_sgn, &ra_h, &ra_m, &dec_sgn, &dec_d, &dec_m,
1319 &source_data[j].dist, &source_data[j].lum, &source_data[j].fudge );
1320 if ( c != 10 )
1321 {
1322 fprintf( stderr, "error parsing source datafile %s\n", sourceFileName );
1323 exit( 1 );
1324 }
1325
1326 /* by convention, overall sign is carried only on hours/degrees entry */
1327 source_data[j].ra = ( ra_h + ra_m / 60.0 ) * LAL_PI / 12.0;
1328 source_data[j].dec = ( dec_d + dec_m / 60.0 ) * LAL_PI / 180.0;
1329
1330 if ( ra_sgn == '-' )
1331 source_data[j].ra *= -1;
1332 if ( dec_sgn == '-' )
1333 source_data[j].dec *= -1;
1334 ++j;
1335 }
1336
1337 /* close file */
1338 fclose( fp );
1339
1340 /* generate ratio and fraction vectors */
1341
1342 ratioVec = calloc( num_source, sizeof( REAL8 ) );
1343 fracVec = calloc( num_source, sizeof( REAL8 ) );
1344 if ( !ratioVec || !fracVec )
1345 {
1346 fprintf( stderr, "Allocation error for ratioVec/fracVec\n" );
1347 exit( 1 );
1348 }
1349
1350 /* MW luminosity might be zero */
1351 norm = mwLuminosity;
1352
1353 /* calculate the fractions of the different sources */
1354 for ( k = 0; k < num_source; ++k )
1355 norm += ratioVec[k] = source_data[k].lum * source_data[k].fudge;
1356 fracVec[0] = ratioVec[0] / norm;
1357 for ( k = 1; k < num_source; ++k )
1358 fracVec[k] = fracVec[k-1] + ratioVec[k] / norm;
1359}
1360
1361/*
1362 *
1363 * Function to read IPN sky simulations from text file given file - read(file,ra,dec)
1364 *
1365 */
1366
1367 void
1368read_IPN_grid_from_file( char *fname )
1369{
1370 UINT4 j; /* counter */
1371 char line[256]; /* string holders */
1372 FILE *data; /* file object */
1373
1374 /* read file */
1375 data = fopen(fname, "r");
1376 if ( ! data )
1377 {
1378 fprintf( stderr, "Could not find file %s\n", fname );
1379 exit( 1 );
1380 }
1381
1382 /* find number of lines */
1383 numSkyPoints = 0;
1384 while ( fgets( line, sizeof( line ), data ) )
1385 ++numSkyPoints;
1386
1387 /* seek to start of file again */
1388 fseek(data, 0, SEEK_SET);
1389
1390 /* assign memory for sky points */
1391 skyPoints = LALCalloc(numSkyPoints, sizeof(*skyPoints));
1392 if ( ! skyPoints )
1393 {
1394 fprintf( stderr, "Allocation error for skyPoints\n" );
1395 exit( 1 );
1396 }
1397
1398 j = 0;
1399 while ( fgets( line, sizeof( line ), data ) )
1400 {
1401 REAL8 ra, dec;
1402 int c;
1403
1404 c = sscanf( line, "%le %le", &ra, &dec );
1405 if ( c != 2 )
1406 {
1407 fprintf( stderr, "error parsing IPN sky points datafile %s\n", IPNSkyPositionsFile );
1408 exit( 1 );
1409 }
1410
1411 /* convert to radians */
1412 skyPoints[j].ra = ra * ( LAL_PI / 180.0 ); /* from degrees (IPN file) to radians */
1413 skyPoints[j].dec = dec * ( LAL_PI / 180.0 );
1414 ++j;
1415 }
1416
1417 /* close file */
1418 fclose( data );
1419}
1420
1421/*
1422 *
1423 * Function to complete galaxy catalog
1424 *
1425 */
1426
1427 void
1428sourceComplete(void)
1429{
1430 /* Catalog Completion Constants */
1431 REAL8 Mbstar = -20.45;
1432 /* Mbstar = magnitude at which the number of galaxies begins to fall off exponentially, */
1433 /* corrected for reddening (to agree with the lum density of 0.0198) */
1434 REAL8 phistar = 0.0081/0.92; /* normalization constant */
1435 REAL8 alpha = -0.9; /* determines slope at faint end of luminosity function */
1436 REAL8 initDistance = 0.0; /* minimum Distance for galaxy catalog*/
1437 REAL8 DeltaD = 100.0; /* Distance step for stepping through galaxy catalog (kPc) */
1438 REAL8 maxDistance = srcCompleteDist; /* Distance to which you want to correct the catalog (kPc)*/
1439 REAL8 M_min = -12.0; /* minimum blue light magnitude */
1440 REAL8 M_max = -25; /* maximum blue light magnitude */
1441 REAL8 edgestep = 0.1; /* magnitude bin size */
1442
1443 /* Vectors */
1444 REAL8Vector *phibins = NULL; /* Magnitude bins for calculating Schechter function */
1445 REAL8Vector *Distance = NULL; /* Distances from initDistance to maxDistance in steps of DeltaD */
1446 REAL8Vector *phi = NULL; /* Schechter magnitude function */
1447 REAL8Vector *phiN = NULL; /* Number of expected galaxies in each magnitude bin */
1448 REAL8Vector *N = NULL; /* Actual number of galaxies in each magnitude bin */
1449 REAL8Vector *pN = NULL; /* Running tally of the fake galaxies added to the catalog */
1450 REAL8Vector *Corrections = NULL; /* Number of galaxies to be added in each magnitude bin */
1451
1452 /* Other Variables */
1453 int edgenum = (int) ceil((M_min-M_max)/edgestep); /* Number of magnitude bins */
1454 const char *galaxyname = "Fake"; /* Beginning of name for all added (non-real) galaxies */
1455 int distnum = (maxDistance-initDistance)/DeltaD; /* Number of elements in Distance vector */
1456 int k_at_25Mpc = floor((25000-initDistance)/DeltaD); /* Initial index for Distance vector - no galaxies added before 25Mpc */
1457 int j,k,q; /* Indices for loops */
1458 REAL8 mag; /* Converted blue light luminosity of each galaxy */
1459 int mag_index; /* Index of each galaxy when binning by magnitude */
1460 FILE *fp; /* File for output of corrected galaxy catalog */
1461 REAL8 pow1 = 0.0; /* Used to calculate Schechter function */
1462 REAL8 pow2 = 0.0; /* Used to calculate Schechter function */
1463
1464 REAL8 UNUSED shellLum = 0.0;
1465
1466 /* Parameters for generating random sky positions */
1467 SimInspiralTable *randPositionTable;
1468 static RandomParams* randPositions = NULL;
1469 int rand_skylocation_seed = 3456;
1470
1471 /* Set up linked list for added galaxies*/
1472 struct FakeGalaxy *myFakeGalaxy;
1473 struct FakeGalaxy *head; /*=myFakeGalaxy;*/
1474 struct FakeGalaxy *saved_next;
1475
1476 /* Create the Vectors */
1477 phibins = XLALCreateREAL8Vector(edgenum);
1478 Distance = XLALCreateREAL8Vector(distnum+1);
1479 phi = XLALCreateREAL8Vector(edgenum);
1480 phiN = XLALCreateREAL8Vector(edgenum);
1481 N = XLALCreateREAL8Vector(edgenum);
1482 pN = XLALCreateREAL8Vector(edgenum);
1483 Corrections = XLALCreateREAL8Vector(edgenum);
1484
1485 /* Initialize sky location parameters and FakeGalaxy linked list */
1486 randPositionTable = calloc(1, sizeof(SimInspiralTable));
1487 LALCreateRandomParams( &status, &randPositions, rand_skylocation_seed);
1488 galaxynum = 0;
1489 myFakeGalaxy = (struct FakeGalaxy*) calloc(1, sizeof(struct FakeGalaxy));
1490 head = myFakeGalaxy;
1491
1492 /* Initialize the vectors */
1493 for (j=0; j<edgenum; j++)
1494 {
1495 phibins->data[j] = M_max+j*edgestep;
1496 phiN->data[j] = 0;
1497 N->data[j] = 0;
1498 pN->data[j] = 0;
1499 Corrections->data[j] = 0;
1500
1501 /* Calculate the theoretical blue light magnitude in each magnitude bin */
1502 pow1 = -1*pow(10, (-0.4*(phibins->data[j]-Mbstar)));
1503 pow2 = pow(10, (-0.4*(phibins->data[j]-Mbstar)));
1504 phi->data[j] = 0.92*phistar*exp(pow1)*pow(pow2, alpha+1);
1505 }
1506
1507 /* Initialize the Distance array */
1508 for (j=0; j<=distnum; j++)
1509 {
1510 Distance->data[j] = initDistance+j*DeltaD;
1511 }
1512
1513 /* Iterate through Distance vector and bin galaxies according to magnitude at each distance */
1514 for (k = k_at_25Mpc; k<distnum; k++)
1515 {
1516
1517 /* Reset N to zero before you count the galaxies with distances less than the current Distance */
1518 for (q = 0; q<edgenum;q++)
1519 {
1520 N->data[q]=0;
1521 }
1522
1523 /* Count the number of galaxies in the spherical volume with radius Distance->data[k+1] and bin them in magnitude */
1524 for( q = 0; q<num_source; q++)
1525 {
1526 if ( (source_data[q].dist<=Distance->data[k+1]) )
1527 {
1528 /* Convert galaxy luminosity to blue light magnitude */
1529 mag = -2.5*(log10(source_data[q].lum)+7.808);
1530 /* Calculate which magnitude bin it falls in */
1531 mag_index = (int) floor((mag-M_max)/edgestep);
1532 /* Create a histogram array of the number of galaxies in each magnitude bin */
1533 if (mag_index >= 0 && mag_index<edgenum)
1534 {
1535 N->data[mag_index] += 1.0;
1536 }
1537 else printf("WARNING GALAXY DOESNT FIT IN BIN\n");
1538 }
1539 }
1540
1541 /* Add galaxies to the catalog based on the difference between the expected number of galaxies and the number of galaxies in the catalog */
1542 for (j = 0; j<edgenum; j++)
1543 {
1544 /* Number of galaxies expected in the spherical volume with radius Distance->data[k+1] */
1545 phiN->data[j] =edgestep*phi->data[j]*(4.0/3.0)*LAL_PI*(pow(Distance->data[k+1]/1000.0,3));
1546 /* Difference between the counted number of galaxies and the expected number of galaxies */
1547 Corrections->data[j] = phiN->data[j] - N->data[j] - pN->data[j];
1548 /* If there are galaxies missing, add them */
1549 if (Corrections->data[j]>0.0)
1550 {
1551 for (q=0;q<floor(Corrections->data[j]);q++)
1552 {
1553 randPositionTable = XLALRandomInspiralSkyLocation( randPositionTable, randPositions);
1554 myFakeGalaxy->dist = Distance->data[k+1];
1555 myFakeGalaxy->ra = randPositionTable->longitude;
1556 myFakeGalaxy->dec = randPositionTable->latitude;
1557 myFakeGalaxy->fudge = 1;
1558 sprintf(myFakeGalaxy->name, "%s%d", galaxyname, galaxynum);
1559 myFakeGalaxy->lum = pow(10.0, (phibins->data[j]/(-2.5)-7.808));
1560 myFakeGalaxy->next = (struct FakeGalaxy*) calloc(1,sizeof(struct FakeGalaxy));
1561 myFakeGalaxy = myFakeGalaxy->next;
1562 galaxynum++;
1563 pN->data[j] += 1.0;
1564 }
1565 }
1566 }
1567 }
1568
1569 /* Combine source_data (original catalog) and FakeGalaxies into one array */
1570 temparray = calloc((num_source+galaxynum), sizeof(*source_data));
1571 if ( !temparray )
1572 { fprintf( stderr, "Allocation error for temparray\n" );
1573 exit( 1 );
1574 }
1575
1576 for (j=0;j<num_source;j++)
1577 {
1578 temparray[j].dist = source_data[j].dist;
1579 temparray[j].lum = source_data[j].lum;
1580 sprintf(temparray[j].name, "%s", source_data[j].name);
1581 temparray[j].ra = source_data[j].ra;
1582 temparray[j].dec = source_data[j].dec;
1583 temparray[j].fudge = source_data[j].fudge;
1584 }
1585 myFakeGalaxy = head;
1586 for (j=num_source;j<(num_source+galaxynum);j++)
1587 {
1588 temparray[j].dist = myFakeGalaxy->dist;
1589 temparray[j].lum = myFakeGalaxy->lum;
1590 sprintf(temparray[j].name, "%s", myFakeGalaxy->name);
1591 temparray[j].ra = myFakeGalaxy->ra;
1592 temparray[j].dec = myFakeGalaxy->dec;
1593 temparray[j].fudge = myFakeGalaxy->fudge;
1594 myFakeGalaxy = myFakeGalaxy->next;
1595 }
1596 myFakeGalaxy->next = NULL;
1597
1598 /* Point old_source_data at source_data */
1599 old_source_data = source_data;
1600
1601 /* Point source_data at the new array*/
1602 source_data = temparray;
1603 shellLum = 0;
1604
1605 if (makeCatalog == 1)
1606 {
1607 /* Write the corrected catalog to a file */
1608 fp = fopen("correctedcatalog.txt", "w+");
1609 for (j=0; j<(num_source+galaxynum);j++) {
1610 fprintf(fp, "%s %g %g %g %g %g \n", source_data[j].name, source_data[j].ra,
1611 source_data[j].dec, source_data[j].dist, source_data[j].lum, source_data[j].fudge );
1612 }
1613 fclose(fp);
1614 }
1615
1616 /* Recalculate some variables from read_source_data that will have changed due to the addition of fake galaxies */
1617 ratioVec = (REAL8*) calloc( (num_source+galaxynum), sizeof( REAL8 ) );
1618 fracVec = (REAL8*) calloc( (num_source+galaxynum), sizeof( REAL8 ) );
1619 if ( !ratioVec || !fracVec )
1620 {
1621 fprintf( stderr, "Allocation error for ratioVec/fracVec\n" );
1622 exit( 1 );
1623 }
1624
1625 /* MW luminosity might be zero */
1626 norm = mwLuminosity;
1627
1628 /* calculate the fractions of the different sources */
1629 for ( i = 0; i <(num_source+galaxynum); ++i )
1630 norm += ratioVec[i] = source_data[i].lum * source_data[i].fudge;
1631 fracVec[0] = ratioVec[0] / norm;
1632 for ( i = 1; i <(num_source+galaxynum); ++i )
1633 fracVec[i] = fracVec[i-1] + ratioVec[i] / norm;
1634
1635 /* Free some stuff */
1636 myFakeGalaxy = head;
1637 for (j=0; j<galaxynum; j++) {
1638 saved_next = myFakeGalaxy->next;
1639 free(myFakeGalaxy);
1640 myFakeGalaxy = saved_next;
1641 }
1642 LALFree(old_source_data);
1643 LALFree( skyPoints );
1644 LALDestroyRandomParams(&status, &randPositions);
1645
1646 XLALDestroyREAL8Vector(phibins);
1647 XLALDestroyREAL8Vector(Corrections);
1648 XLALDestroyREAL8Vector(Distance);
1649 XLALDestroyREAL8Vector(phi);
1650 XLALDestroyREAL8Vector(phiN);
1651 XLALDestroyREAL8Vector(N);
1652 XLALDestroyREAL8Vector(pN);
1653}
1654
1655/*
1656 *
1657 * functions to draw masses from source distribution
1658 *
1659 */
1660
1661 void
1662drawMassFromSource( SimInspiralTable* table )
1663{
1664 REAL4 m1, m2, eta;
1665 int mass_index=0;
1666
1667 /* choose masses from the mass-list */
1668 mass_index = (int)( num_mass * XLALUniformDeviate( randParams ) );
1669 m1 = redshift_mass(mass_data[mass_index].mass1,redshift);
1670 m2 = redshift_mass(mass_data[mass_index].mass2,redshift);
1671
1672 eta=m1 * m2 / ( ( m1 + m2 ) * ( m1 + m2 ) );
1673 table->mass1 = m1;
1674 table->mass2 = m2;
1675 table->eta = eta;
1676 table->mchirp = pow( eta, 0.6) * (m1 + m2);
1677}
1678
1679
1680/*
1681 *
1682 * functions to draw masses and spins from NR distribution
1683 *
1684 */
1685
1686 void
1687drawMassSpinFromNR( SimInspiralTable* table )
1688{
1689 int mass_index=0;
1690
1691 /* choose masses from the mass-list */
1692 mass_index = (int)( num_nr * XLALUniformDeviate( randParams ) );
1693 XLALRandomNRInjectTotalMass( table, randParams, minMtotal, maxMtotal,
1694 nrSimArray[mass_index]);
1695}
1696
1697
1698 void
1699drawMassSpinFromNRNinja2( SimInspiralTable* inj )
1700{
1701 /* For ninja2 we first select a mass, then find */
1702 /* a waveform that can be injected at that mass */
1703
1704 int j,k;
1705 REAL8 startFreq, startFreqHz, massTotal;
1706 int indx,tmp,*indicies;
1707
1708 /* Permute the indicies in a random order */
1709 /* This lets us check each available waveform */
1710 /* once and lets us know when no option works */
1711 indicies = (int *) LALCalloc( num_nr, sizeof(int) );
1712
1713 for ( j = 0; j < num_nr; j++ )
1714 indicies[j] = j;
1715
1716 for ( j = 0; j < num_nr; j++ )
1717 {
1718 indx = (int) ( (num_nr-j) * XLALUniformDeviate( randParams ) ) + j;
1719 tmp = indicies[j];
1720 indicies[j] = indicies[indx];
1721 indicies[indx] = tmp;
1722 }
1723
1724 massTotal = (maxMtotal - minMtotal) * XLALUniformDeviate( randParams ) + minMtotal;
1725
1726 for ( j = 0; j < num_nr; j++ )
1727 {
1728 k = indicies[j];
1729 if (nrSimArray[k]->f_lower > 0.0000001)
1730 startFreq = nrSimArray[k]->f_lower;
1731 else
1732 startFreq = start_freq_from_frame_url(nrSimArray[k]->numrel_data);
1733 /*startFreqHz = startFreq / (LAL_TWOPI * massTotal * LAL_MTSUN_SI);*/
1734 startFreqHz = startFreq / (massTotal);
1735 /* if this startFreqHz makes us happy, inject it */
1736 if (startFreqHz <= inj->f_lower)
1737 {
1738 /* This is a copy of XLALRandomNRInjectTotalMass without */
1739 /* the random mass selection. TODO: refactor that method */
1740 inj->eta = nrSimArray[k]->eta;
1741 inj->mchirp = massTotal * pow(inj->eta, 3.0/5.0);
1742
1743 /* set mass1 and mass2 */
1744 inj->mass1 = (massTotal / 2.0) * (1 + pow( (1 - 4 * inj->eta), 0.5) );
1745 inj->mass2 = (massTotal / 2.0) * (1 - pow( (1 - 4 * inj->eta), 0.5) );
1746
1747 /* copy over the spin parameters */
1748 inj->spin1x = nrSimArray[k]->spin1x;
1749 inj->spin1y = nrSimArray[k]->spin1y;
1750 inj->spin1z = nrSimArray[k]->spin1z;
1751 inj->spin2x = nrSimArray[k]->spin2x;
1752 inj->spin2y = nrSimArray[k]->spin2y;
1753 inj->spin2z = nrSimArray[k]->spin2z;
1754
1755 /* copy over the numrel information */
1756 inj->numrel_mode_min = nrSimArray[k]->numrel_mode_min;
1757 inj->numrel_mode_max = nrSimArray[k]->numrel_mode_max;
1758 snprintf(inj->numrel_data, LIGOMETA_STRING_MAX, "%s",
1759 nrSimArray[k]->numrel_data);
1760
1761 XLALFree(indicies);
1762 return;
1763 }
1764 }
1765
1766 /* If we hit the end of the list, oops */
1767 XLALFree(indicies);
1768 /* should throw an error here... */
1769 fprintf(stderr,"No waveform could be injected at MTotal=%f Msun\n", massTotal);
1770}
1771
1772/*
1773 *
1774 * functions to draw sky location from source distribution
1775 *
1776 */
1777
1778 void
1779drawFromSource(
1780 REAL8 *rightAscension,
1781 REAL8 *declination,
1782 REAL8 *distance,
1783 CHAR name[LIGOMETA_SOURCE_MAX] )
1784{
1785 REAL4 u;
1786 int j;
1787
1788 u=XLALUniformDeviate( randParams );
1789
1790 /* draw from the source table */
1791 for ( j = 0; j < num_source; ++j )
1792 {
1793 if ( u < fracVec[j] )
1794 {
1795 /* put the parameters */
1796 *rightAscension = source_data[j].ra;
1797 *declination = source_data[j].dec;
1798 *distance = source_data[j].dist/1000.0;
1799 memcpy( name, source_data[j].name,
1800 sizeof(CHAR) * LIGOMETA_SOURCE_MAX );
1801 return;
1802 }
1803 }
1804
1805 /* now then, draw from MilkyWay
1806 * WARNING: This sets location AND distance */
1807 XLALRandomInspiralMilkywayLocation( rightAscension, declination, distance,
1808 randParams );
1809 memcpy( name, MW_name, sizeof(CHAR) * 30 );
1810}
1811
1812/*
1813 *
1814 * function to draw IPN sky location from IPN simulation points
1815 *
1816 */
1817
1818 void
1819drawFromIPNsim(
1820 REAL8 *rightAscension,
1821 REAL8 *declination )
1822{
1823 REAL4 u;
1824 INT4 j;
1825
1826 u=XLALUniformDeviate( randParams );
1827 j=( int ) (u*numSkyPoints);
1828
1829 /* draw from the IPN source table */
1830 if ( j < numSkyPoints )
1831 {
1832 /* put the parameters */
1833 *rightAscension = skyPoints[j].ra;
1834 *declination = skyPoints[j].dec;
1835 return;
1836 }
1837}
1838
1839/*
1840 *
1841 * function to draw sky location from exttrig source file
1842 *
1843 */
1844
1845 void
1846drawLocationFromExttrig( SimInspiralTable* table )
1847{
1848 LIGOTimeGPS timeGRB; /* real time of the GRB */
1849 REAL4 ra_rad, de_rad;
1850 REAL8 gmst1, gmst2;
1851
1852 /* convert the position (stored as degree) to radians first */
1853 ra_rad = exttrigHead->event_ra * LAL_PI_180;
1854 de_rad = exttrigHead->event_dec * LAL_PI_180;
1855
1856 /* populate the time structures */
1857 timeGRB.gpsSeconds = exttrigHead->start_time;
1858 timeGRB.gpsNanoSeconds = exttrigHead->start_time_ns;
1859
1860 gmst1 = XLALGreenwichMeanSiderealTime(&timeGRB);
1861 gmst2 = XLALGreenwichMeanSiderealTime(&table->geocent_end_time);
1862
1863 /* populate the table */
1864 table->longitude = ra_rad- gmst1 + gmst2;
1865 table->latitude = de_rad;
1866}
1867
1868
1869/*
1870 *
1871 * generate all parameters (sky position and angles) for a random inspiral
1872 *
1873 */
1874
1875int main( int argc, char *argv[] )
1876{
1877 LIGOTimeGPS gpsStartTime = {-1,0};
1878 LIGOTimeGPS gpsEndTime = {-1,0};
1879 LIGOTimeGPS IPNgpsTime = {-1,0};
1880 LIGOTimeGPS currentGpsTime;
1881 long gpsDuration;
1882
1883 REAL8 meanTimeStep = -1;
1884 REAL8 timeInterval = 0;
1885 REAL4 fLower = -1;
1886 UINT4 useChirpDist = 0;
1887 REAL4 minMass10, maxMass10, minMass20, maxMass20, minMtotal0, maxMtotal0,
1888 meanMass10, meanMass20, massStdev10, massStdev20; /* masses at z=0 */
1889 REAL8 pzmax=0; /* maximal value of the probability distribution of the redshift */
1890 INT4 ncount;
1891 size_t ninj;
1892 int rand_seed = 1;
1893
1894 /* waveform */
1895 CHAR waveform[LIGOMETA_WAVEFORM_MAX];
1896 CHAR axisChoiceString[]="angmomentum";
1897 CHAR dummy[256];
1898 INT4 amp_order = -1;
1899 /* xml output data */
1900 CHAR fname[256];
1901 CHAR *userTag = NULL;
1902 ProcessTable *proctable;
1903 ProcessParamsTable *procparams;
1904 SimInspiralTable *injections;
1905 SimRingdownTable *ringparams;
1906 ProcessParamsTable *this_proc_param;
1907 LIGOLwXMLStream *xmlfp;
1908
1909 REAL8 drawnDistance = 0.0;
1910 REAL8 drawnRightAscension = 0.0;
1911 REAL8 drawnDeclination = 0.0;
1912 CHAR drawnSourceName[LIGOMETA_SOURCE_MAX];
1913 REAL8 IPNgmst1 = 0.0;
1914 REAL8 IPNgmst2 = 0.0;
1915
1916 REAL8 targetSNR;
1917
1918 CHAR *ligoPsdFileName = NULL;
1919 REAL8 ligoStartFreq = -1;
1920 CHAR *virgoPsdFileName = NULL;
1921 REAL8 virgoStartFreq = -1;
1922 CHAR *ligoFakePsd=NULL;
1923 CHAR *virgoFakePsd=NULL;
1924 REAL8FrequencySeries *ligoPsd = NULL;
1925 REAL8FrequencySeries *virgoPsd = NULL;
1926 XLAL_INIT_MEM(status);
1927
1928 /* LALgetopt arguments */
1929 struct LALoption long_options[] =
1930 {
1931 {"help", no_argument, 0, 'h'},
1932 {"verbose", no_argument, &vrbflg, 1 },
1933 {"source-file", required_argument, 0, 'f'},
1934 {"mass-file", required_argument, 0, 'm'},
1935 {"nr-file", required_argument, 0, 'c'},
1936 {"exttrig-file", required_argument, 0, 'E'},
1937 {"f-lower", required_argument, 0, 'F'},
1938 {"gps-start-time", required_argument, 0, 'a'},
1939 {"gps-end-time", required_argument, 0, 'b'},
1940 {"ipn-gps-time", required_argument, 0, '"'},
1941 {"t-distr", required_argument, 0, '('},
1942 {"time-step", required_argument, 0, 't'},
1943 {"time-interval", required_argument, 0, 'i'},
1944 {"time-file", required_argument, 0, 1035},
1945 {"seed", required_argument, 0, 's'},
1946 {"waveform", required_argument, 0, 'w'},
1947 {"amp-order", required_argument, 0, 'q'},
1948 {"user-tag", required_argument, 0, 'Z'},
1949 {"userTag", required_argument, 0, 'Z'},
1950 {"m-distr", required_argument, 0, 'd'},
1951 {"min-mass1", required_argument, 0, 'j'},
1952 {"max-mass1", required_argument, 0, 'k'},
1953 {"min-mass2", required_argument, 0, 'J'},
1954 {"max-mass2", required_argument, 0, 'K'},
1955 {"min-mtotal", required_argument, 0, 'A'},
1956 {"max-mtotal", required_argument, 0, 'L'},
1957 {"fixed-mass1", required_argument, 0, ']'},
1958 {"fixed-mass2", required_argument, 0, '['},
1959 {"mean-mass1", required_argument, 0, 'n'},
1960 {"mean-mass2", required_argument, 0, 'N'},
1961 {"ninja2-mass", no_argument, &ninjaMass, 1},
1962 {"real8-ninja2", no_argument, &real8Ninja2, 1},
1963 {"mass1-points", required_argument, 0, ':'},
1964 {"mass2-points", required_argument, 0, ';'},
1965 {"stdev-mass1", required_argument, 0, 'o'},
1966 {"stdev-mass2", required_argument, 0, 'O'},
1967 {"min-mratio", required_argument, 0, 'x'},
1968 {"max-mratio", required_argument, 0, 'y'},
1969 {"d-distr", required_argument, 0, 'e'},
1970 {"min-distance", required_argument, 0, 'p'},
1971 {"max-distance", required_argument, 0, 'r'},
1972 {"dchirp-distr", required_argument, 0, ','},
1973 {"z-distr", required_argument, 0, '5'},
1974 {"min-z", required_argument, 0, '6'},
1975 {"max-z", required_argument, 0, '7'},
1976 {"local-rate", required_argument, 0, ')'},
1977 {"snr-distr", required_argument, 0, '1'},
1978 {"min-snr", required_argument, 0, '2'},
1979 {"max-snr", required_argument, 0, '3'},
1980 {"min-coinc-snr", required_argument, 0, 1707},
1981 {"ifos", required_argument, 0, '4'},
1982 {"ninja-snr", no_argument, &ninjaSNR, 1},
1983 {"ligo-psd", required_argument, 0, 500},
1984 {"ligo-fake-psd", required_argument, 0, 501},
1985 {"ligo-start-freq", required_argument, 0, 502},
1986 {"virgo-psd", required_argument, 0, 600},
1987 {"virgo-fake-psd", required_argument, 0, 601},
1988 {"virgo-start-freq", required_argument, 0, 602},
1989 {"l-distr", required_argument, 0, 'l'},
1990 {"longitude", required_argument, 0, 'v'},
1991 {"latitude", required_argument, 0, 'z'},
1992 {"i-distr", required_argument, 0, 'I'},
1993 {"incl-std", required_argument, 0, 'B'},
1994 {"fixed-inc", required_argument, 0, 'C'},
1995 {"max-inc", required_argument, 0, 1001},
1996 {"polarization", required_argument, 0, 'S'},
1997 {"coa-phase-distr", required_argument, 0, 1007},
1998 {"fixed-coa-phase", required_argument, 0, 1008},
1999 {"sourcecomplete", required_argument, 0, 'H'},
2000 {"make-catalog", no_argument, 0, '.'},
2001 {"enable-milkyway", required_argument, 0, 'M'},
2002 {"disable-milkyway", no_argument, 0, 'D'},
2003 {"min-spin1", required_argument, 0, 'g'},
2004 {"min-kappa1", required_argument, 0, 'Q'},
2005 {"max-kappa1", required_argument, 0, 'R'},
2006 {"min-abskappa1", required_argument, 0, 'X'},
2007 {"max-abskappa1", required_argument, 0, 'Y'},
2008 {"max-spin1", required_argument, 0, 'G'},
2009 {"min-spin2", required_argument, 0, 'u'},
2010 {"max-spin2", required_argument, 0, 'U'},
2011 {"output", required_argument, 0, 'P'},
2012 {"version", no_argument, 0, 'V'},
2013 {"enable-spin", no_argument, 0, 'T'},
2014 {"disable-spin", no_argument, 0, 'W'},
2015 {"aligned", no_argument, 0, '@'},
2016 {"axis-choice", required_argument, 0, 1009},
2017 {"write-compress", no_argument, &outCompress, 1},
2018 {"taper-injection", required_argument, 0, '*'},
2019 {"band-pass-injection", no_argument, 0, '}'},
2020 {"ipn-file", required_argument, 0, '^'},
2021 {"spin-gaussian", no_argument, 0, 1002},
2022 {"stdev-spin1", required_argument, 0, 1003},
2023 {"stdev-spin2", required_argument, 0, 1004},
2024 {"mean-spin1", required_argument, 0, 1005},
2025 {"mean-spin2", required_argument, 0, 1006},
2026 {0, 0, 0, 0}
2027 };
2028 int c;
2029
2030 /* set up initial debugging values */
2031 lal_errhandler = LAL_ERR_EXIT;
2032
2033 /* create the process and process params tables */
2034 proctable = (ProcessTable *) calloc( 1, sizeof(ProcessTable) );
2035 XLALGPSTimeNow(&(proctable->start_time));
2036 XLALPopulateProcessTable(proctable, PROGRAM_NAME, lalAppsVCSIdentInfo.vcsId,
2037 lalAppsVCSIdentInfo.vcsStatus, lalAppsVCSIdentInfo.vcsDate, 0);
2038 snprintf( proctable->comment, LIGOMETA_COMMENT_MAX, " " );
2039 this_proc_param = procparams = (ProcessParamsTable *)
2040 calloc( 1, sizeof(ProcessParamsTable) );
2041
2042 /* clear the waveform field */
2043 memset( waveform, 0, LIGOMETA_WAVEFORM_MAX * sizeof(CHAR) );
2044
2045 /* parse the arguments */
2046 while ( 1 )
2047 {
2048 /* LALgetopt_long stores long option here */
2049 int option_index = 0;
2050 long int gpsinput;
2051 size_t LALoptarg_len;
2052
2053 c = LALgetopt_long_only( argc, argv,
2054 "hf:m:a:b:t:s:w:i:M:*", long_options, &option_index );
2055
2056 /* detect the end of the options */
2057 if ( c == - 1 )
2058 {
2059 break;
2060 }
2061
2062 switch ( c )
2063 {
2064 case 0:
2065 /* if this option set a flag, do nothing else now */
2066 if ( long_options[option_index].flag != 0 )
2067 {
2068 break;
2069 }
2070 else
2071 {
2072 fprintf( stderr, "error parsing option %s with argument %s\n",
2073 long_options[option_index].name, LALoptarg );
2074 exit( 1 );
2075 }
2076 break;
2077
2078 case 'f':
2079 LALoptarg_len = strlen( LALoptarg ) + 1;
2080 sourceFileName = calloc( 1, LALoptarg_len * sizeof(char) );
2081 memcpy( sourceFileName, LALoptarg, LALoptarg_len * sizeof(char) );
2082 this_proc_param = this_proc_param->next =
2083 next_process_param( long_options[option_index].name, "string",
2084 "%s", LALoptarg );
2085 break;
2086
2087 case 'm':
2088 LALoptarg_len = strlen( LALoptarg ) + 1;
2089 massFileName = calloc( 1, LALoptarg_len * sizeof(char) );
2090 memcpy( massFileName, LALoptarg, LALoptarg_len * sizeof(char) );
2091 this_proc_param = this_proc_param->next =
2092 next_process_param( long_options[option_index].name, "string",
2093 "%s", LALoptarg );
2094 break;
2095
2096 case 'c':
2097 LALoptarg_len = strlen( LALoptarg ) + 1;
2098 nrFileName = calloc( 1, LALoptarg_len * sizeof(char) );
2099 memcpy( nrFileName, LALoptarg, LALoptarg_len * sizeof(char) );
2100 this_proc_param = this_proc_param->next =
2101 next_process_param( long_options[option_index].name, "string",
2102 "%s", LALoptarg );
2103 break;
2104
2105 case 1035:
2106 LALoptarg_len = strlen( LALoptarg ) + 1;
2107 injtimesFileName = calloc( 1, LALoptarg_len * sizeof(char) );
2108 memcpy( injtimesFileName, LALoptarg, LALoptarg_len * sizeof(char) );
2109 this_proc_param = this_proc_param->next =
2110 next_process_param( long_options[option_index].name, "string",
2111 "%s", LALoptarg );
2112 break;
2113
2114 case 'E':
2115 LALoptarg_len = strlen( LALoptarg ) + 1;
2116 exttrigFileName = calloc( 1, LALoptarg_len * sizeof(char) );
2117 memcpy( exttrigFileName, LALoptarg, LALoptarg_len * sizeof(char) );
2118 this_proc_param = this_proc_param->next =
2119 next_process_param( long_options[option_index].name, "string",
2120 "%s", LALoptarg );
2121 break;
2122
2123 case 'F':
2124 fLower = atof( LALoptarg );
2125 this_proc_param = this_proc_param->next =
2126 next_process_param( long_options[option_index].name, "float",
2127 "%f", fLower );
2128 break;
2129
2130 case 'a':
2131 gpsinput = atol( LALoptarg );
2132 if ( gpsinput < 441417609 )
2133 {
2134 fprintf( stderr, "invalid argument to --%s:\n"
2135 "GPS start time is prior to "
2136 "Jan 01, 1994 00:00:00 UTC:\n"
2137 "(%ld specified)\n",
2138 long_options[option_index].name, gpsinput );
2139 exit( 1 );
2140 }
2141 gpsStartTime.gpsSeconds = gpsinput;
2142 gpsStartTime.gpsNanoSeconds = 0;
2143 this_proc_param = this_proc_param->next =
2144 next_process_param( long_options[option_index].name, "int",
2145 "%ld", gpsinput );
2146 break;
2147
2148 case 'b':
2149 gpsinput = atol( LALoptarg );
2150 if ( gpsinput < 441417609 )
2151 {
2152 fprintf( stderr, "invalid argument to --%s:\n"
2153 "GPS start time is prior to "
2154 "Jan 01, 1994 00:00:00 UTC:\n"
2155 "(%ld specified)\n",
2156 long_options[option_index].name, gpsinput );
2157 exit( 1 );
2158 }
2159 gpsEndTime.gpsSeconds = gpsinput;
2160 gpsEndTime.gpsNanoSeconds = 0;
2161 this_proc_param = this_proc_param->next =
2162 next_process_param( long_options[option_index].name, "int",
2163 "%ld", gpsinput );
2164 break;
2165
2166 case '"':
2167 gpsinput = atol( LALoptarg );
2168 if ( gpsinput < 441417609 )
2169 {
2170 fprintf( stderr, "invalid argument to --%s:\n"
2171 "GPS start time is prior to "
2172 "Jan 01, 1994 00:00:00 UTC:\n"
2173 "(%ld specified)\n",
2174 long_options[option_index].name, gpsinput );
2175 exit( 1 );
2176 }
2177 IPNgpsTime.gpsSeconds = gpsinput;
2178 this_proc_param = this_proc_param->next =
2179 next_process_param( long_options[option_index].name, "int",
2180 "%ld", gpsinput );
2181 break;
2182
2183 case 's':
2184 rand_seed = atoi( LALoptarg );
2185 this_proc_param = this_proc_param->next =
2186 next_process_param( long_options[option_index].name, "int",
2187 "%d", rand_seed );
2188 break;
2189
2190 case '(':
2191 LALoptarg_len = strlen( LALoptarg ) + 1;
2192 memcpy( dummy, LALoptarg, LALoptarg_len );
2193
2194 if (!strcmp(dummy, "fixed"))
2195 {
2196 tDistr=LALINSPIRAL_FIXED_TIME_DIST;
2197 }
2198 else if (!strcmp(dummy, "uniform"))
2199 {
2200 tDistr=LALINSPIRAL_UNIFORM_TIME_DIST;
2201 }
2202 else if (!strcmp(dummy, "exponential"))
2203 {
2204 tDistr=LALINSPIRAL_EXPONENTIAL_TIME_DIST;
2205 }
2206 else if (!strcmp(dummy, "file"))
2207 {
2208 tDistr=LALINSPIRAL_FILE_TIME_DIST;
2209 }
2210 else
2211 {
2212 tDistr=LALINSPIRAL_UNKNOWN_TIME_DIST;
2213 fprintf( stderr, "invalid argument to --%s:\n"
2214 "unknown time distribution: %s must be one of\n"
2215 "fixed, uniform or exponential\n",
2216 long_options[option_index].name, LALoptarg );
2217 exit( 1 );
2218 }
2219 break;
2220
2221 case ')':
2222 localRate = atof( LALoptarg );
2223 this_proc_param = this_proc_param->next =
2224 next_process_param( long_options[option_index].name, "float",
2225 "%le", localRate );
2226 if ( !(localRate > 0.) )
2227 {
2228 fprintf( stderr, "invalid argument to --%s:\n"
2229 "local coalescence rate must be positive"
2230 "(%f specified)\n",
2231 long_options[option_index].name, localRate );
2232 exit( 1 );
2233 }
2234 break;
2235
2236 case 't':
2237 meanTimeStep = atof( LALoptarg );
2238 this_proc_param = this_proc_param->next =
2239 next_process_param( long_options[option_index].name, "float",
2240 "%le", meanTimeStep );
2241 break;
2242
2243 case 'i':
2244 timeInterval = atof( LALoptarg );
2245 this_proc_param = this_proc_param->next =
2246 next_process_param( long_options[option_index].name, "float",
2247 "%le", timeInterval );
2248 break;
2249
2250 case 'w':
2251 snprintf( waveform, LIGOMETA_WAVEFORM_MAX, "%s", LALoptarg );
2252 this_proc_param = this_proc_param->next =
2253 next_process_param( long_options[option_index].name, "string",
2254 "%s", LALoptarg );
2255 break;
2256
2257 case 'q':
2258 amp_order = atof( LALoptarg );
2259 this_proc_param = this_proc_param->next =
2260 next_process_param( long_options[option_index].name, "int",
2261 "%ld", amp_order );
2262 break;
2263
2264 case 'M':
2265 /* set the luminosity of the Milky Way */
2266 mwLuminosity = atof( LALoptarg );
2267 if ( mwLuminosity < 0 )
2268 {
2269 fprintf( stderr, "invalid argument to --%s:\n"
2270 "Milky Way luminosity must be positive"
2271 "(%f specified)\n",
2272 long_options[option_index].name, mwLuminosity );
2273 exit( 1 );
2274 }
2275
2276 this_proc_param = this_proc_param->next =
2277 next_process_param( long_options[option_index].name, "float",
2278 "%le", mwLuminosity );
2279 break;
2280
2281 case 'D':
2282 /* set the luminosity of the Milky Way */
2283 this_proc_param = this_proc_param->next =
2284 next_process_param( long_options[option_index].name, "string",
2285 "" );
2286 mwLuminosity = 0;
2287 break;
2288
2289 case 'Z':
2290 /* create storage for the usertag */
2291 LALoptarg_len = strlen( LALoptarg ) + 1;
2292 userTag = (CHAR *) calloc( LALoptarg_len, sizeof(CHAR) );
2293 memcpy( userTag, LALoptarg, LALoptarg_len );
2294
2295 this_proc_param = this_proc_param->next = (ProcessParamsTable *)
2296 calloc( 1, sizeof(ProcessParamsTable) );
2297 snprintf( this_proc_param->program, LIGOMETA_PROGRAM_MAX, "%s",
2298 PROGRAM_NAME );
2299 snprintf( this_proc_param->param, LIGOMETA_PARAM_MAX, "--userTag" );
2300 snprintf( this_proc_param->type, LIGOMETA_TYPE_MAX, "string" );
2301 snprintf( this_proc_param->value, LIGOMETA_VALUE_MAX, "%s",
2302 LALoptarg );
2303 break;
2304
2305 case 'd':
2306 LALoptarg_len = strlen( LALoptarg ) + 1;
2307 memcpy( dummy, LALoptarg, LALoptarg_len );
2308 this_proc_param = this_proc_param->next = (ProcessParamsTable *)
2309 calloc( 1, sizeof(ProcessParamsTable) );
2310 snprintf( this_proc_param->program, LIGOMETA_PROGRAM_MAX, "%s",
2311 PROGRAM_NAME );
2312 snprintf( this_proc_param->param, LIGOMETA_PARAM_MAX, "--m-distr" );
2313 snprintf( this_proc_param->type, LIGOMETA_TYPE_MAX, "string" );
2314 snprintf( this_proc_param->value, LIGOMETA_VALUE_MAX, "%s",
2315 LALoptarg );
2316
2317 if (!strcmp(dummy, "source"))
2318 {
2319 mDistr=massFromSourceFile;
2320 }
2321 else if (!strcmp(dummy, "nrwaves"))
2322 {
2323 mDistr=massFromNRFile;
2324 }
2325 else if (!strcmp(dummy, "totalMass"))
2326 {
2327 mDistr=uniformTotalMass;
2328 }
2329 else if (!strcmp(dummy, "componentMass"))
2330 {
2331 mDistr=uniformComponentMass;
2332 }
2333 else if (!strcmp(dummy, "gaussian"))
2334 {
2335 mDistr=gaussianMassDist;
2336 }
2337 else if (!strcmp(dummy, "log"))
2338 {
2339 mDistr=logComponentMass;
2340 }
2341 else if (!strcmp(dummy, "totalMassRatio"))
2342 {
2343 mDistr=uniformTotalMassRatio;
2344 }
2345 else if (!strcmp(dummy, "logTotalMassUniformMassRatio"))
2346 {
2347 mDistr=logMassUniformTotalMassRatio;
2348 }
2349 else if (!strcmp(dummy, "m1m2SquareGrid"))
2350 {
2351 mDistr=m1m2SquareGrid;
2352 }
2353 else if (!strcmp(dummy, "fixMasses"))
2354 {
2355 mDistr=fixMasses;
2356 }
2357 else if (!strcmp(dummy, "totalMassFraction"))
2358 {
2359 mDistr=uniformTotalMassFraction;
2360 }
2361 else
2362 {
2363 fprintf( stderr, "invalid argument to --%s:\n"
2364 "unknown mass distribution: %s must be one of\n"
2365 "(source, nrwaves, totalMass, componentMass, gaussian, log,\n"
2366 "totalMassRatio, totalMassFraction, logTotalMassUniformMassRatio,\n"
2367 "m1m2SquareGrid, fixMasses)\n",
2368 long_options[option_index].name, LALoptarg );
2369 exit( 1 );
2370 }
2371 break;
2372
2373 case 'j':
2374 minMass1 = atof( LALoptarg );
2375 this_proc_param = this_proc_param->next =
2376 next_process_param( long_options[option_index].name,
2377 "float", "%le", minMass1 );
2378 break;
2379
2380 case 'k':
2381 maxMass1 = atof( LALoptarg );
2382 this_proc_param = this_proc_param->next =
2383 next_process_param( long_options[option_index].name,
2384 "float", "%le", maxMass1 );
2385 break;
2386
2387 case 'J':
2388 minMass2 = atof( LALoptarg );
2389 this_proc_param = this_proc_param->next =
2390 next_process_param( long_options[option_index].name,
2391 "float", "%le", minMass2 );
2392 break;
2393
2394 case 'K':
2395 maxMass2 = atof( LALoptarg );
2396 this_proc_param = this_proc_param->next =
2397 next_process_param( long_options[option_index].name,
2398 "float", "%le", maxMass2 );
2399 break;
2400
2401 case 'A':
2402 minMtotal = atof( LALoptarg );
2403 this_proc_param = this_proc_param->next =
2404 next_process_param( long_options[option_index].name,
2405 "float", "%le", minMtotal );
2406 break;
2407
2408 case 'L':
2409 maxMtotal = atof( LALoptarg );
2410 this_proc_param = this_proc_param->next =
2411 next_process_param( long_options[option_index].name,
2412 "float", "%le", maxMtotal );
2413 break;
2414
2415 case 'n':
2416 meanMass1 = atof( LALoptarg );
2417 this_proc_param = this_proc_param->next =
2418 next_process_param( long_options[option_index].name,
2419 "float", "%le", meanMass1 );
2420 break;
2421
2422 case 'N':
2423 meanMass2 = atof( LALoptarg );
2424 this_proc_param = this_proc_param->next =
2425 next_process_param( long_options[option_index].name,
2426 "float", "%le", meanMass2 );
2427 break;
2428
2429 case 'o':
2430 massStdev1 = atof( LALoptarg );
2431 this_proc_param = this_proc_param->next =
2432 next_process_param( long_options[option_index].name,
2433 "float", "%le", massStdev1 );
2434 break;
2435
2436 case 'O':
2437 massStdev2 = atof( LALoptarg );
2438 this_proc_param = this_proc_param->next =
2439 next_process_param( long_options[option_index].name,
2440 "float", "%le", massStdev2 );
2441 break;
2442
2443 case 'x':
2444 minMassRatio = atof( LALoptarg );
2445 this_proc_param = this_proc_param->next =
2446 next_process_param( long_options[option_index].name,
2447 "float", "%le", minMassRatio );
2448 break;
2449
2450 case 'y':
2451 maxMassRatio = atof( LALoptarg );
2452 this_proc_param = this_proc_param->next =
2453 next_process_param( long_options[option_index].name,
2454 "float", "%le", maxMassRatio );
2455 break;
2456
2457 case ':':
2458 pntMass1 = atof( LALoptarg );
2459 this_proc_param = this_proc_param->next =
2460 next_process_param( long_options[option_index].name,
2461 "int", "%d", pntMass1 );
2462 break;
2463
2464 case ';':
2465 pntMass2 = atof( LALoptarg );
2466 this_proc_param = this_proc_param->next =
2467 next_process_param( long_options[option_index].name,
2468 "int", "%d", pntMass2 );
2469 break;
2470
2471 case ']':
2472 fixedMass1 = atof( LALoptarg );
2473 this_proc_param = this_proc_param->next =
2474 next_process_param( long_options[option_index].name,
2475 "float", "%f", fixedMass1 );
2476 break;
2477
2478 case '[':
2479 fixedMass2 = atof( LALoptarg );
2480 this_proc_param = this_proc_param->next =
2481 next_process_param( long_options[option_index].name,
2482 "float", "%f", fixedMass2 );
2483 break;
2484
2485 case 'e':
2486 LALoptarg_len = strlen( LALoptarg ) + 1;
2487 memcpy( dummy, LALoptarg, LALoptarg_len );
2488 this_proc_param = this_proc_param->next = (ProcessParamsTable *)
2489 calloc( 1, sizeof(ProcessParamsTable) );
2490 snprintf( this_proc_param->program, LIGOMETA_PROGRAM_MAX, "%s",
2491 PROGRAM_NAME );
2492 snprintf( this_proc_param->param,LIGOMETA_PARAM_MAX,"--d-distr" );
2493 snprintf( this_proc_param->type, LIGOMETA_TYPE_MAX, "string" );
2494 snprintf( this_proc_param->value,LIGOMETA_VALUE_MAX,"%s", LALoptarg );
2495 haveLoudness += 1; /* counter to check for clashing options */
2496
2497 if (!strcmp(dummy, "source"))
2498 {
2499 dDistr=distFromSourceFile;
2500 }
2501 else if (!strcmp(dummy, "uniform"))
2502 {
2503 dDistr=uniformDistance;
2504 }
2505 else if (!strcmp(dummy, "distancesquared"))
2506 {
2507 dDistr=uniformDistanceSquared;
2508 }
2509 else if (!strcmp(dummy, "log10"))
2510 {
2511 dDistr=uniformLogDistance;
2512 }
2513 else if (!strcmp(dummy, "volume"))
2514 {
2515 dDistr=uniformVolume;
2516 }
2517 else
2518 {
2519 fprintf( stderr, "invalid argument to --%s:\n"
2520 "unknown distance distribution: "
2521 "%s, must be one of (uniform, distancesquared, volume, log10, source)\n",
2522 long_options[option_index].name, LALoptarg );
2523 exit( 1 );
2524 }
2525 break;
2526
2527 case ',':
2528 LALoptarg_len = strlen( LALoptarg ) + 1;
2529 memcpy( dummy, LALoptarg, LALoptarg_len );
2530 this_proc_param = this_proc_param->next = (ProcessParamsTable *)
2531 calloc( 1, sizeof(ProcessParamsTable) );
2532 snprintf( this_proc_param->program, LIGOMETA_PROGRAM_MAX, "%s",
2533 PROGRAM_NAME );
2534 snprintf( this_proc_param->param,LIGOMETA_PARAM_MAX,"--dchirp-distr" );
2535 snprintf( this_proc_param->type, LIGOMETA_TYPE_MAX, "string" );
2536 snprintf( this_proc_param->value,LIGOMETA_VALUE_MAX,"%s", LALoptarg );
2537 haveLoudness += 1; /* counter to check for clashing options */
2538 useChirpDist = 1;
2539
2540 if (!strcmp(dummy, "uniform"))
2541 {
2542 dDistr=uniformDistance;
2543 }
2544 else if (!strcmp(dummy, "distancesquared"))
2545 {
2546 dDistr=uniformDistanceSquared;
2547 }
2548 else if (!strcmp(dummy, "log10"))
2549 {
2550 dDistr=uniformLogDistance;
2551 }
2552 else if (!strcmp(dummy, "volume"))
2553 {
2554 dDistr=uniformVolume;
2555 }
2556 else
2557 {
2558 fprintf( stderr, "invalid argument to --%s:\n"
2559 "unknown distribution: "
2560 "%s, must be one of (uniform, distancesquared, volume, log10)\n",
2561 long_options[option_index].name, LALoptarg );
2562 exit( 1 );
2563 }
2564 break;
2565
2566 case 'p':
2567 /* minimum distance from earth */
2568 minD = (REAL4) atof( LALoptarg );
2569 if ( minD <= 0 )
2570 {
2571 fprintf( stderr, "invalid argument to --%s:\n"
2572 "minimum distance must be > 0: "
2573 "(%f kpc specified)\n",
2574 long_options[option_index].name, minD );
2575 exit( 1 );
2576 }
2577 this_proc_param = this_proc_param->next =
2578 next_process_param( long_options[option_index].name,
2579 "float", "%e", minD );
2580 break;
2581
2582 case 'r':
2583 /* max distance from earth */
2584 maxD = (REAL4) atof( LALoptarg );
2585 if ( maxD <= 0 )
2586 {
2587 fprintf( stderr, "invalid argument to --%s:\n"
2588 "maximum distance must be greater than 0: "
2589 "(%f kpc specified)\n",
2590 long_options[option_index].name, maxD );
2591 exit( 1 );
2592 }
2593 this_proc_param = this_proc_param->next =
2594 next_process_param( long_options[option_index].name,
2595 "float", "%e", maxD );
2596 break;
2597
2598 case '5':
2599 LALoptarg_len = strlen( LALoptarg ) + 1;
2600 memcpy( dummy, LALoptarg, LALoptarg_len );
2601 this_proc_param = this_proc_param->next = (ProcessParamsTable *)
2602 calloc( 1, sizeof(ProcessParamsTable) );
2603 snprintf( this_proc_param->program, LIGOMETA_PROGRAM_MAX, "%s",
2604 PROGRAM_NAME );
2605 snprintf( this_proc_param->param,LIGOMETA_PARAM_MAX,"--z-distr" );
2606 snprintf( this_proc_param->type, LIGOMETA_TYPE_MAX, "string" );
2607 snprintf( this_proc_param->value,LIGOMETA_VALUE_MAX,"%s", LALoptarg );
2608 haveLoudness += 1; /* counter to check for clashing options */
2609
2610 if (!strcmp(dummy, "sfr"))
2611 {
2612 dDistr = starFormationRate;
2613 }
2614 else
2615 {
2616 fprintf( stderr, "invalid argument to --%s:\n"
2617 "unknown redshift distribution: "
2618 "%s, must be sfr (other distributions may be implemented in future)\n",
2619 long_options[option_index].name, LALoptarg );
2620 exit( 1 );
2621 }
2622 break;
2623
2624 case '6':
2625 minZ = atof( LALoptarg );
2626 this_proc_param = this_proc_param->next =
2627 next_process_param( long_options[option_index].name,
2628 "float", "%le", minZ );
2629 if ( minZ < 0 )
2630 {
2631 fprintf(stderr,"invalid argument to --%s:\n"
2632 "%s must not be less than 0.\n",
2633 long_options[option_index].name, LALoptarg );
2634 exit( 1 );
2635 }
2636 break;
2637
2638 case '7':
2639 maxZ = atof( LALoptarg );
2640 this_proc_param = this_proc_param->next =
2641 next_process_param( long_options[option_index].name,
2642 "float", "%le", maxZ );
2643 if ( maxZ < 0 )
2644 {
2645 fprintf(stderr,"invalid argument to --%s:\n"
2646 "%s must not be less than 0.\n",
2647 long_options[option_index].name, LALoptarg );
2648 exit( 1 );
2649 }
2650 break;
2651
2652 case '1':
2653 LALoptarg_len = strlen( LALoptarg ) + 1;
2654 memcpy( dummy, LALoptarg, LALoptarg_len );
2655 this_proc_param = this_proc_param->next = (ProcessParamsTable *)
2656 calloc( 1, sizeof(ProcessParamsTable) );
2657 snprintf( this_proc_param->program, LIGOMETA_PROGRAM_MAX, "%s",
2658 PROGRAM_NAME );
2659 snprintf( this_proc_param->param,LIGOMETA_PARAM_MAX,"--snr-distr" );
2660 snprintf( this_proc_param->type, LIGOMETA_TYPE_MAX, "string" );
2661 snprintf( this_proc_param->value,LIGOMETA_VALUE_MAX,"%s", LALoptarg );
2662 haveLoudness += 1; /* counter to check for clashing options */
2663
2664 if (!strcmp(dummy, "uniform"))
2665 {
2666 dDistr=uniformSnr;
2667 }
2668 else if (!strcmp(dummy, "log10"))
2669 {
2670 dDistr=uniformLogSnr;
2671 }
2672 else if (!strcmp(dummy, "volume"))
2673 {
2674 dDistr=uniformVolumeSnr;
2675 }
2676 else
2677 {
2678 fprintf( stderr, "invalid argument to --%s:\n"
2679 "unknown SNR distribution: "
2680 "%s, must be uniform, log10, or volume \n",
2681 long_options[option_index].name, LALoptarg );
2682 exit( 1 );
2683 }
2684 break;
2685
2686 case '2':
2687 minSNR = atof( LALoptarg );
2688 this_proc_param = this_proc_param->next =
2689 next_process_param( long_options[option_index].name,
2690 "float", "%le", minSNR );
2691 if ( minSNR < 2 )
2692 {
2693 fprintf(stderr,"invalid argument to --%s:\n"
2694 "%s must be greater than 2\n",
2695 long_options[option_index].name, LALoptarg );
2696 exit( 1 );
2697 }
2698 break;
2699
2700 case '3':
2701 maxSNR = atof( LALoptarg );
2702 this_proc_param = this_proc_param->next =
2703 next_process_param( long_options[option_index].name,
2704 "float", "%le", maxSNR );
2705 if ( maxSNR < 2 )
2706 {
2707 fprintf(stderr,"invalid argument to --%s:\n"
2708 "%s must be greater than 2\n",
2709 long_options[option_index].name, LALoptarg );
2710 exit( 1 );
2711 }
2712 break;
2713
2714 case '4':
2715 LALoptarg_len = strlen( LALoptarg ) + 1;
2716 ifos = calloc( 1, LALoptarg_len * sizeof(char) );
2717 memcpy( ifos, LALoptarg, LALoptarg_len * sizeof(char) );
2718 this_proc_param = this_proc_param->next =
2719 next_process_param( long_options[option_index].name, "string",
2720 "%s", LALoptarg );
2721 break;
2722
2723 case 'l':
2724 LALoptarg_len = strlen( LALoptarg ) + 1;
2725 memcpy( dummy, LALoptarg, LALoptarg_len );
2726 this_proc_param = this_proc_param->next = (ProcessParamsTable *)
2727 calloc( 1, sizeof(ProcessParamsTable) );
2728 snprintf( this_proc_param->program, LIGOMETA_PROGRAM_MAX, "%s",
2729 PROGRAM_NAME );
2730 snprintf( this_proc_param->param, LIGOMETA_PARAM_MAX, "--l-distr" );
2731 snprintf( this_proc_param->type, LIGOMETA_TYPE_MAX, "string" );
2732 snprintf( this_proc_param->value, LIGOMETA_VALUE_MAX, "%s",
2733 LALoptarg );
2734
2735 if (!strcmp(dummy, "source"))
2736 {
2737 lDistr=locationFromSourceFile;
2738 }
2739 else if (!strcmp(dummy, "exttrig"))
2740 {
2741 lDistr=locationFromExttrigFile;
2742 }
2743 else if (!strcmp(dummy, "random"))
2744 {
2745 lDistr=uniformSkyLocation;
2746 }
2747 else if (!strcmp(dummy, "fixed"))
2748 {
2749 lDistr=fixedSkyLocation;
2750 }
2751 else if (!strcmp(dummy, "ipn"))
2752 {
2753 lDistr=locationFromIPNFile;
2754 }
2755 else
2756 {
2757 fprintf( stderr, "invalid argument to --%s:\n"
2758 "unknown location distribution: "
2759 "%s must be one of (source, random)\n",
2760 long_options[option_index].name, LALoptarg );
2761 exit( 1 );
2762 }
2763
2764 break;
2765
2766 case 'H':
2767 /* Turn on galaxy catalog completion function */
2768 srcComplete = 1;
2769 srcCompleteDist = (REAL8) atof( LALoptarg );
2770 this_proc_param = this_proc_param->next =
2771 next_process_param( long_options[option_index].name,
2772 "string", "%s", LALoptarg );
2773 break;
2774
2775 case '.':
2776 /* Create a text file of completed catalog */
2777 makeCatalog = 1;
2778 break;
2779
2780 case 'v':
2781 /* fixed location (longitude) */
2782 longitude = atof( LALoptarg )*LAL_PI_180 ;
2783 if (longitude <= ( LAL_PI + epsAngle ) && \
2784 longitude >= ( -LAL_PI - epsAngle ))
2785 {
2786 this_proc_param = this_proc_param->next =
2787 next_process_param( long_options[option_index].name,
2788 "float", "%e", longitude );
2789 }
2790 else
2791 {
2792 fprintf(stderr,"invalid argument to --%s:\n"
2793 "%s must be between -180. and 180. degrees\n",
2794 long_options[option_index].name, LALoptarg );
2795 exit( 1 );
2796 }
2797 break;
2798
2799 case 'z':
2800 /* fixed location (latitude) */
2801 latitude = (REAL4) atof( LALoptarg )*LAL_PI_180;
2802 if (latitude <= ( LAL_PI/2. + epsAngle ) && \
2803 latitude >= ( -LAL_PI/2. - epsAngle ))
2804 {
2805 this_proc_param = this_proc_param->next =
2806 next_process_param( long_options[option_index].name,
2807 "float", "%e", latitude );
2808 }
2809 else
2810 {
2811 fprintf(stderr,"invalid argument to --%s:\n"
2812 "%s must be between -90. and 90. degrees\n",
2813 long_options[option_index].name, LALoptarg );
2814 exit( 1 );
2815 }
2816 break;
2817
2818 case 'I':
2819 LALoptarg_len = strlen( LALoptarg ) + 1;
2820 memcpy( dummy, LALoptarg, LALoptarg_len );
2821 this_proc_param = this_proc_param->next = (ProcessParamsTable *)
2822 calloc( 1, sizeof(ProcessParamsTable) );
2823 snprintf( this_proc_param->program, LIGOMETA_PROGRAM_MAX, "%s",
2824 PROGRAM_NAME );
2825 snprintf( this_proc_param->param, LIGOMETA_PARAM_MAX, "--i-distr" );
2826 snprintf( this_proc_param->type, LIGOMETA_TYPE_MAX, "string" );
2827 snprintf( this_proc_param->value, LIGOMETA_VALUE_MAX, "%s",
2828 LALoptarg );
2829
2830 if (!strcmp(dummy, "uniform"))
2831 {
2832 iDistr=uniformInclDist;
2833 }
2834 else if (!strcmp(dummy, "gaussian"))
2835 {
2836 iDistr=gaussianInclDist;
2837 }
2838 else if (!strcmp(dummy, "fixed"))
2839 {
2840 iDistr=fixedInclDist;
2841 }
2842 else
2843 {
2844 fprintf( stderr, "invalid argument to --%s:\n"
2845 "unknown inclination distribution: "
2846 "%s must be one of (uniform, gaussian, fixed)\n",
2847 long_options[option_index].name, LALoptarg );
2848 exit( 1 );
2849 }
2850 break;
2851
2852 case 'B':
2853 /* gaussian width for inclination */
2854 inclStd = (REAL4) atof( LALoptarg );
2855 if ( inclStd <= 0 )
2856 {
2857 fprintf( stderr, "invalid argument to --%s:\n"
2858 "inclination gaussian width must be greater than 0: "
2859 "(%f specified)\n",
2860 long_options[option_index].name, inclStd );
2861 exit( 1 );
2862 }
2863 this_proc_param = this_proc_param->next =
2864 next_process_param( long_options[option_index].name,
2865 "float", "%e", inclStd );
2866 break;
2867
2868 case 'C':
2869 /* fixed angle of inclination */
2870 fixed_inc = (REAL4) atof( LALoptarg )/180.*LAL_PI;
2871 this_proc_param = this_proc_param->next =
2872 next_process_param( long_options[option_index].name,
2873 "float", "%e", fixed_inc );
2874 break;
2875
2876 case 1001:
2877 /* maximum angle of inclination */
2878 max_inc = (REAL4) atof( LALoptarg )/180.*LAL_PI;
2879 if ( (atof(LALoptarg) < 0.) || (atof(LALoptarg) >= 180.) ) {
2880 fprintf( stderr, "invalid argument to --%s:\n"
2881 "maximum inclination angle must be between 0 and 180 degrees:"
2882 "(%s specified)\n",
2883 long_options[option_index].name, LALoptarg );
2884 exit( 1 );
2885 }
2886 this_proc_param = this_proc_param->next =
2887 next_process_param( long_options[option_index].name,
2888 "float", "%e", max_inc );
2889 break;
2890
2891 case 1007:
2892 /* coalescence phase distribution */
2893 if ( strcmp( LALoptarg, "uniform" ) == 0)
2894 coaPhaseFixed = 0;
2895 else if ( strcmp( LALoptarg, "fixed" ) == 0)
2896 coaPhaseFixed = 1;
2897 else {
2898 fprintf( stderr, "invalid argument to --%s:\n"
2899 "must either uniform or fixed (%s specified)\n",
2900 long_options[option_index].name, LALoptarg );
2901 exit( 1 );
2902 }
2903 this_proc_param = this_proc_param->next = (ProcessParamsTable *)
2904 calloc( 1, sizeof(ProcessParamsTable) );
2905 snprintf( this_proc_param->program, LIGOMETA_PROGRAM_MAX, "%s",
2906 PROGRAM_NAME );
2907 snprintf( this_proc_param->param, LIGOMETA_PARAM_MAX, "--coa-phase-distr" );
2908 snprintf( this_proc_param->type, LIGOMETA_TYPE_MAX, "string" );
2909 snprintf( this_proc_param->value, LIGOMETA_VALUE_MAX, "%s",
2910 LALoptarg );
2911 break;
2912
2913 case 1008:
2914 /* fixed coalescence phase */
2915 fixedCoaPhase = (REAL4) atof( LALoptarg );
2916 if ( (fixedCoaPhase < 0.) || (fixedCoaPhase >= 360.) ) {
2917 fprintf( stderr, "invalid argument to --%s:\n"
2918 "fixed coalescence phase must be between 0 and 360 degrees:"
2919 "(%s specified)\n",
2920 long_options[option_index].name, LALoptarg );
2921 exit( 1 );
2922 }
2923 this_proc_param = this_proc_param->next =
2924 next_process_param( long_options[option_index].name,
2925 "float", "%e", fixedCoaPhase );
2926 fixedCoaPhase *= LAL_PI / 180.;
2927 break;
2928
2929 case 'S':
2930 /* set the polarization angle */
2931 psi = (REAL4) atof( LALoptarg )/180.*LAL_PI;
2932 if ( (atof(LALoptarg) < 0.) || (atof(LALoptarg) >= 360.) ) {
2933 fprintf( stderr, "invalid argument to --%s:\n"
2934 "polarization angle must be between 0 and 360 degrees: "
2935 "(%s specified)\n",
2936 long_options[option_index].name, LALoptarg );
2937 exit( 1 );
2938 }
2939 this_proc_param = this_proc_param->next =
2940 next_process_param( long_options[option_index].name,
2941 "float", "%e", psi );
2942 break;
2943
2944 case 'P':
2945 LALoptarg_len = strlen( LALoptarg ) + 1;
2946 outputFileName = calloc( 1, LALoptarg_len * sizeof(char) );
2947 memcpy( outputFileName, LALoptarg, LALoptarg_len * sizeof(char) );
2948 this_proc_param = this_proc_param->next =
2949 next_process_param( long_options[option_index].name,
2950 "string", "%s", LALoptarg );
2951 break;
2952
2953 case 500: /* LIGO psd file */
2954 LALoptarg_len = strlen( LALoptarg ) + 1;
2955 ligoPsdFileName = calloc( 1, LALoptarg_len * sizeof(char) );
2956 memcpy( ligoPsdFileName, LALoptarg, LALoptarg_len * sizeof(char) );
2957 this_proc_param = this_proc_param->next =
2958 next_process_param( long_options[option_index].name,
2959 "string", "%s", LALoptarg );
2960 break;
2961
2962 case 501: /* LIGO fake LALSim PSD */
2963 LALoptarg_len = strlen( LALoptarg ) + 1;
2964 ligoFakePsd = calloc( 1, LALoptarg_len * sizeof(char) );
2965 memcpy( ligoFakePsd, LALoptarg, LALoptarg_len * sizeof(char) );
2966 this_proc_param = this_proc_param->next =
2967 next_process_param( long_options[option_index].name,
2968 "string", "%s", LALoptarg );
2969 break;
2970
2971 case 502: /* LIGO start frequency */
2972 ligoStartFreq = (REAL8) atof( LALoptarg );
2973 this_proc_param = this_proc_param->next =
2974 next_process_param( long_options[option_index].name,
2975 "float", "%f", ligoStartFreq );
2976 break;
2977
2978 case 600: /* Virgo psd file */
2979 LALoptarg_len = strlen( LALoptarg ) + 1;
2980 virgoPsdFileName = calloc( 1, LALoptarg_len * sizeof(char) );
2981 memcpy( virgoPsdFileName, LALoptarg, LALoptarg_len * sizeof(char) );
2982 this_proc_param = this_proc_param->next =
2983 next_process_param( long_options[option_index].name,
2984 "string", "%s", LALoptarg );
2985 break;
2986
2987 case 601: /* Virgo fake LALSim PSD */
2988 LALoptarg_len = strlen( LALoptarg ) + 1;
2989 virgoFakePsd = calloc( 1, LALoptarg_len * sizeof(char) );
2990 memcpy( virgoFakePsd, LALoptarg, LALoptarg_len * sizeof(char) );
2991 this_proc_param = this_proc_param->next =
2992 next_process_param( long_options[option_index].name,
2993 "string", "%s", LALoptarg );
2994 break;
2995
2996 case 602: /* Virgo start frequency */
2997 virgoStartFreq = (REAL8) atof( LALoptarg );
2998 this_proc_param = this_proc_param->next =
2999 next_process_param( long_options[option_index].name,
3000 "float", "%f", virgoStartFreq );
3001 break;
3002
3003 case 1707: /* Set min coincident SNR in two IFOs */
3004 single_IFO_SNR_threshold=(REAL8) atof(LALoptarg);
3005 this_proc_param = this_proc_param->next =
3006 next_process_param( long_options[option_index].name,
3007 "float", "%e", single_IFO_SNR_threshold );
3008 break;
3009
3010 case 'g':
3011 minSpin1 = atof( LALoptarg );
3012 this_proc_param = this_proc_param->next =
3013 next_process_param( long_options[option_index].name,
3014 "float", "%le", minSpin1 );
3015 break;
3016
3017 case 'G':
3018 maxSpin1 = atof( LALoptarg );
3019 this_proc_param = this_proc_param->next =
3020 next_process_param( long_options[option_index].name,
3021 "float", "%le", maxSpin1 );
3022 break;
3023
3024 case 'Q':
3025 minKappa1 = atof( LALoptarg );
3026 this_proc_param = this_proc_param->next =
3027 next_process_param( long_options[option_index].name,
3028 "float", "%le", minKappa1 );
3029 break;
3030
3031 case 'R':
3032 maxKappa1 = atof( LALoptarg );
3033 this_proc_param = this_proc_param->next =
3034 next_process_param( long_options[option_index].name,
3035 "float", "%le", maxKappa1 );
3036 break;
3037
3038 case 'X':
3039 minabsKappa1 = atof( LALoptarg );
3040 this_proc_param = this_proc_param->next =
3041 next_process_param( long_options[option_index].name,
3042 "float", "%le", minabsKappa1 );
3043 break;
3044
3045 case 'Y':
3046 maxabsKappa1 = atof( LALoptarg );
3047 this_proc_param = this_proc_param->next =
3048 next_process_param( long_options[option_index].name,
3049 "float", "%le", maxabsKappa1 );
3050 break;
3051
3052 case 'u':
3053 minSpin2 = atof( LALoptarg );
3054 this_proc_param = this_proc_param->next =
3055 next_process_param( long_options[option_index].name,
3056 "float", "%le", minSpin2 );
3057 break;
3058
3059 case 'U':
3060 maxSpin2 = atof( LALoptarg );
3061 this_proc_param = this_proc_param->next =
3062 next_process_param( long_options[option_index].name,
3063 "float", "%le", maxSpin2 );
3064 break;
3065
3066 case 'V':
3067 /* print version information and exit */
3068 fprintf( stdout, "LIGO/LSC inspiral injection engine\n");
3069 XLALOutputVCSInfo(stderr, lalAppsVCSInfoList, 0, "%% ");
3070 exit( 0 );
3071 break;
3072
3073 case 'T':
3074 /* enable spinning injections */
3075 this_proc_param = this_proc_param->next =
3076 next_process_param( long_options[option_index].name, "string",
3077 "" );
3078 spinInjections = 1;
3079 break;
3080
3081 case 'W':
3082 /* disable spinning injections */
3083 this_proc_param = this_proc_param->next =
3084 next_process_param( long_options[option_index].name, "string",
3085 "" );
3086 spinInjections = 0;
3087 break;
3088
3089 case '@':
3090 /* enforce aligned spins */
3091 this_proc_param = this_proc_param->next =
3092 next_process_param( long_options[option_index].name, "string",
3093 "" );
3094 spinAligned = 1;
3095 break;
3096
3097 case 1009:
3098 /* frame axis choice */
3099 snprintf( axisChoiceString, AXIS_MAX, "%s", LALoptarg );
3100 this_proc_param = this_proc_param->next =
3101 next_process_param( long_options[option_index].name, "string",
3102 "%s", LALoptarg );
3103 break;
3104
3105 case '}':
3106 /* enable band-passing */
3107 this_proc_param = this_proc_param->next =
3108 next_process_param( long_options[option_index].name, "string",
3109 "" );
3110 bandPassInj = 1;
3111 break;
3112
3113 case '*':
3114 /* Set injection tapering */
3115 if ( ! strcmp( "start", LALoptarg ) )
3116 {
3117 taperInj = LAL_SIM_INSPIRAL_TAPER_START;
3118 }
3119 else if ( ! strcmp( "end", LALoptarg ) )
3120 {
3121 taperInj = LAL_SIM_INSPIRAL_TAPER_END;
3122 }
3123 else if ( ! strcmp( "startend", LALoptarg ) )
3124 {
3125 taperInj = LAL_SIM_INSPIRAL_TAPER_STARTEND;
3126 }
3127 else
3128 {
3129 fprintf( stderr, "invalid argument to --%s:\n"
3130 "unknown option specified: %s\n"
3131 "(Must be one of start|end|startend)\n",
3132 long_options[option_index].name, LALoptarg );
3133 }
3134 this_proc_param = this_proc_param->next =
3135 next_process_param( long_options[option_index].name,
3136 "string", LALoptarg );
3137 break;
3138
3139 case 'h':
3140 print_usage(argv[0]);
3141 exit( 0 );
3142 break;
3143
3144 case '?':
3145 print_usage(argv[0]);
3146 exit( 1 );
3147 break;
3148
3149 case '^':
3150 LALoptarg_len = strlen( LALoptarg ) + 1;
3151 IPNSkyPositionsFile = calloc( 1, LALoptarg_len * sizeof(char) );
3152 memcpy( IPNSkyPositionsFile, LALoptarg, LALoptarg_len * sizeof(char) );
3153 this_proc_param = this_proc_param->next =
3154 next_process_param( long_options[option_index].name, "string",
3155 "%s", LALoptarg );
3156 break;
3157
3158 case 1002:
3159 spinDistr = gaussianSpinDist;
3160 this_proc_param = this_proc_param->next =
3161 next_process_param( long_options[option_index].name, "string", "" );
3162 break;
3163
3164 case 1003:
3165 Spin1Std = atof( LALoptarg );
3166 this_proc_param = this_proc_param->next =
3167 next_process_param( long_options[option_index].name,
3168 "float", "%le", Spin1Std );
3169 break;
3170
3171 case 1004:
3172 Spin2Std = atof( LALoptarg );
3173 this_proc_param = this_proc_param->next =
3174 next_process_param( long_options[option_index].name,
3175 "float", "%le", Spin2Std );
3176 break;
3177 case 1005:
3178 meanSpin1 = atof( LALoptarg );
3179 this_proc_param = this_proc_param->next =
3180 next_process_param( long_options[option_index].name,
3181 "float", "%le", meanSpin1 );
3182 break;
3183 case 1006:
3184 meanSpin2 = atof( LALoptarg );
3185 this_proc_param = this_proc_param->next =
3186 next_process_param( long_options[option_index].name,
3187 "float", "%le", meanSpin2 );
3188 break;
3189
3190
3191 default:
3192 fprintf( stderr, "unknown error while parsing options\n" );
3193 print_usage(argv[0]);
3194 exit( 1 );
3195 }
3196 }
3197
3198 /* must set MW flag */
3199 if ( mwLuminosity < 0 && dDistr == distFromSourceFile )
3200 {
3201 fprintf( stderr,
3202 "Must specify either --enable-milkyway LUM or --disable-milkyway\n"\
3203 " when using --d-distr=source\n" );
3204 exit( 1 );
3205 }
3206
3207 if ( (gpsStartTime.gpsSeconds==-1 || gpsEndTime.gpsSeconds==-1) && tDistr != LALINSPIRAL_FILE_TIME_DIST)
3208 {
3209 fprintf( stderr,
3210 "Must specify both --gps-start-time and --gps-end-time.\n");
3211 exit( 1 );
3212 }
3213
3214 gpsDuration=gpsEndTime.gpsSeconds-gpsStartTime.gpsSeconds;
3215
3216 if ( (dDistr == unknownLoudnessDist) || (haveLoudness != 1) )
3217 {
3218 fprintf(stderr,"Must specify exactly one distribution out of\n"\
3219 "--d-distr, --dchirp-distr, --z-distr or --snr-distr.\n");
3220 exit( 1 );
3221 }
3222
3223 if ( lDistr == unknownLocationDist )
3224 {
3225 fprintf(stderr,"Must specify a location distribution (--l-distr).\n");
3226 exit( 1 );
3227 }
3228
3229 if ( lDistr == fixedSkyLocation && longitude == 181. )
3230 {
3231 fprintf(stderr,
3232 "Must specify both --longitude and --latitude when using \n"\
3233 "--l-distr=fixed\n");
3234 exit( 1 );
3235 }
3236
3237 if ( lDistr == fixedSkyLocation && latitude == 91. )
3238 {
3239 fprintf(stderr,
3240 "Must specify both --longitude and --latitude when using \n"\
3241 "--l-distr=fixed\n");
3242 exit( 1 );
3243 }
3244
3245 if ( mDistr == unknownMassDist )
3246 {
3247 fprintf(stderr,"Must specify a mass distribution (--m-distr).\n");
3248 exit( 1 );
3249 }
3250
3251 if ( iDistr == unknownInclDist )
3252 {
3253 fprintf(stderr,"Must specify an inclination distribution (--i-distr).\n");
3254 exit( 1 );
3255 }
3256
3257 /* if using source file, check that file and MW choice selected */
3258 if ( dDistr==distFromSourceFile || lDistr==locationFromSourceFile )
3259 {
3260 if ( ! sourceFileName )
3261 {
3262 fprintf( stderr,
3263 "Must specify --source-file when using --d-distr or --l-distr source \n" );
3264 exit( 1 );
3265 }
3266
3267 if ( ( dDistr == distFromSourceFile ) && ( minD>0.0 || maxD>0.0 ) )
3268 {
3269 fprintf( stderr,
3270 "Cannot specify --min-distance or --max-distance\n"\
3271 "if --d-distr=source\n");
3272 exit( 1 );
3273 }
3274
3275 /* read the source distribution here */
3276 read_source_data( sourceFileName );
3277
3278 /* complete the galaxy catalog */
3279 if (srcComplete == 1)
3280 {
3281 sourceComplete();
3282 }
3283 }
3284
3285 /* if using IPN sky points file, check that file exists and read it */
3286 if ( lDistr==locationFromIPNFile )
3287 {
3288 if ( ! IPNSkyPositionsFile )
3289 {
3290 fprintf( stderr,
3291 "Must specify --ipn-file when using IPN sky points distribution\n" );
3292 exit( 1 );
3293 }
3294
3295 /* read the source distribution here */
3296 read_IPN_grid_from_file( IPNSkyPositionsFile );
3297 }
3298
3299 /* check compatibility of distance/loudness options */
3300 if ( ( dDistr == uniformDistance || dDistr == uniformDistanceSquared ||
3301 dDistr == uniformLogDistance || dDistr == uniformVolume ) &&
3302 ( minD<=0.0 || maxD<=0.0 ) )
3303 {
3304 fprintf( stderr,
3305 "Positive minimum and maximum distances must be specified\n");
3306 exit( 1 );
3307 }
3308 if ( dDistr == uniformDistance || dDistr == uniformDistanceSquared ||
3309 dDistr == uniformLogDistance || dDistr == uniformVolume ||
3310 dDistr == distFromSourceFile )
3311 {
3312 if ( minZ>0.0 || maxZ>0.0 || localRate>0.0 || ninjaSNR ||
3313 minSNR>0.0 || maxSNR>0.0 || ifos!=NULL )
3314 {
3315 fprintf( stderr,
3316 "One or more options on redshift or SNR are incompatible\n"\
3317 "with --d-distr or --dchirp-distr !\n");
3318 exit( 1 );
3319 }
3320 }
3321 if ( dDistr == starFormationRate )
3322 {
3323 if ( minD>0.0 || maxD>0.0 || ninjaSNR || minSNR>0.0 || maxSNR>0.0 || ifos!=NULL )
3324 {
3325 fprintf( stderr,
3326 "One or more options on distance or SNR are incompatible\n"\
3327 "with --z-distr !\n");
3328 exit( 1 );
3329 }
3330 if ( minZ<0.2 || maxZ>1.0 )
3331 {
3332 fprintf( stderr,
3333 "Redshift can only take values between 0.2 and 1 for --z-distr=sfr\n");
3334 exit( 1 );
3335 }
3336 if ( localRate<=0. )
3337 {
3338 fprintf( stderr,
3339 "Local coalescence rate must be positive for --z-distr=sfr\n");
3340 exit( 1 );
3341 }
3342 if ( meanTimeStep>=0. )
3343 {
3344 fprintf( stderr, "Time step cannot be specified for --z-distr=sfr\n"\
3345 "(it is calculated from local coalescence rate)\n");
3346 exit( 1 );
3347 }
3348 }
3349 if ( ( dDistr == uniformSnr || dDistr == uniformLogSnr ||
3350 dDistr == uniformVolumeSnr ) &&
3351 ( minD>0.0 || maxD>0.0 || minZ>0.0 || maxZ>0.0 || localRate>0.0 ) )
3352 {
3353 fprintf( stderr,
3354 "One or more options on distance or redshift are incompatible\n"\
3355 "with --snr-distr !\n");
3356 exit( 1 );
3357 }
3358
3359 /* SNR distribution options */
3360
3361 if ( dDistr == uniformSnr || dDistr == uniformLogSnr ||
3362 dDistr == uniformVolumeSnr )
3363 {
3364 /* basic input checks */
3365 if ( minSNR == -1 || maxSNR == -1 || ifos == NULL )
3366 {
3367 fprintf( stderr,
3368 "Must provide all of --min-snr, --max-snr and --ifos to distribute by SNR\n" );
3369 exit( 1 );
3370 }
3371 if (single_IFO_SNR_threshold > maxSNR)
3372 {
3373 fprintf(stderr,
3374 "Minimum coincident SNR should not be larger than maximum network SNR. Exiting...\n");
3375 exit(1);
3376 }
3377
3378 /* Check that each ifo has its PSD file or fakePSD */
3379 char *tmp, *ifo;
3380
3381 /* Get the number and names of IFOs */
3382 tmp = LALCalloc(1, strlen(ifos) + 1);
3383 strcpy(tmp, ifos);
3384 ifo = strtok(tmp,",");
3385
3386 while (ifo != NULL)
3387 {
3388 numifos += 1;
3389 ifo = strtok(NULL, ",");
3390 }
3391
3392 ifonames=realloc(ifonames,(numifos+2)*sizeof(CHAR **));
3393 strcpy(tmp, ifos);
3394 ifo = strtok(tmp,",");
3395 ifonames[0]=malloc(strlen(ifo)+1);
3396 sprintf(ifonames[0],"%s",ifo);
3397
3398 i=1;
3399 while (ifo != NULL)
3400 {
3401 ifo = strtok(NULL, ",");
3402 if (ifo!=NULL)
3403 {
3404 ifonames[i]=malloc(strlen(ifo)+1);
3405 sprintf(ifonames[i],"%s",ifo);
3406 }
3407 i++;
3408 }
3409
3410 ifonames[numifos]=NULL;
3411 i=0;
3412 while (ifonames[i]!= NULL)
3413 {
3414 if (!strcmp(ifonames[i],"H1") || !strcmp(ifonames[i],"L1"))
3415 {
3416 /* Check either PSD file or fakePSD are given */
3417 if(!( ligoFakePsd || ligoPsdFileName ))
3418 {
3419 fprintf( stderr,
3420 "Must provide PSD file or the name of analytic PSD for LIGO if --snr-distr is given and H1 or L1 are in --ifos. \n" );
3421 exit( 1 );
3422 }
3423 /* Check we didn't give both fake PSD and filename*/
3424 if ( ligoFakePsd && ligoPsdFileName )
3425 {
3426 fprintf( stderr,"Must provide only one between --ligo-psd and --ligo-fake-psd \n" );
3427 exit( 1 );
3428 }
3429 /* Check flow for SNR calculation was given */
3430 if (ligoStartFreq < 0)
3431 {
3432 fprintf( stderr, "Must specify --ligo-start-freq together with --ligo-psd or --ligo-fake-psd.\n");
3433 exit( 1 );
3434 }
3435 /* Ninja only work with PSD file, check user provided it */
3436 if (ninjaSNR && !(ligoPsdFileName))
3437 {
3438 fprintf( stderr, "Ninja injections do not support SNR calculation with simulated PSD.\n"
3439 "Please provide PSD file for LIGO (with --ligo-psd filename.dat). Exiting...\n");
3440 exit( 1 );
3441 }
3442 }
3443 else if (!strcmp(ifonames[i],"V1"))
3444 {
3445 /* Check either PSD file or fakePSD are given */
3446 if (!(virgoFakePsd || virgoPsdFileName ))
3447 {
3448 fprintf( stderr,
3449 "Must provide PSD file or the name of analytic PSD for Virgo if --snr-distr is given and V1 is in --ifos. \n" );
3450 exit( 1 );
3451 }
3452 /* Check we didn't give both fake PSD and filename*/
3453 if ( virgoFakePsd && virgoPsdFileName )
3454 {
3455 fprintf( stderr,"Must provide only one between --virgo-psd and --virgo-fake-psd \n" );
3456 exit( 1 );
3457 }
3458 /* Check flow for SNR calculation was given */
3459 if (virgoStartFreq < 0)
3460 {
3461 fprintf( stderr,"Must specify --virgo-start-freq with --virgo-psd or --virgo-fake-psd.\n");
3462 exit( 1 );
3463 }
3464 /* Ninja only work with PSD file, check user provided it */
3465 if (ninjaSNR && !(virgoPsdFileName))
3466 {
3467 fprintf( stderr, "Ninja injections do not support SNR calculation with simulated PSD.\n"
3468 "Please provide PSD file for Virgo (with --virgo-psd filename.dat). Exiting...\n");
3469 exit( 1 );
3470 }
3471 }
3472 i++;
3473 }
3474
3475 if (tmp) LALFree(tmp);
3476 if (ifo) LALFree(ifo);
3477
3478 if ( maxSNR < minSNR )
3479 {
3480 fprintf( stderr, "max SNR must be greater than min SNR\n");
3481 exit( 1 );
3482 }
3483 if (single_IFO_SNR_threshold<0.0)
3484 {
3485 fprintf( stderr,
3486 "The single IFO SNR threshold must be positive. Exiting...\n" );
3487 exit( 1 );
3488 }
3489
3490 /* Check custom PSDs */
3491 if (ligoPsdFileName) {
3492 if (XLALPsdFromFile(&ligoPsd, ligoPsdFileName) != XLAL_SUCCESS)
3493 {
3494 fprintf(stderr, "Unable to load PSD file %s.\n", ligoPsdFileName);
3495 exit( 1 );
3496 }
3497 /* We're done with the filename */
3498 free(ligoPsdFileName);
3499 }
3500
3501 if (virgoPsdFileName) {
3502 if (XLALPsdFromFile(&virgoPsd, virgoPsdFileName) != XLAL_SUCCESS)
3503 {
3504 fprintf(stderr, "Unable to load PSD file %s.\n", virgoPsdFileName);
3505 exit( 1 );
3506 }
3507 /* We're done with the filename */
3508 free(virgoPsdFileName);
3509 }
3510 }
3511
3512 /* check if the source file is specified for distance but NOT for
3513 location */
3514 if ( dDistr==distFromSourceFile && lDistr!=locationFromSourceFile )
3515 {
3516 fprintf( stderr,
3517 "WARNING: source file specified for distance "
3518 "but NOT for location. This might give strange distributions\n" );
3519 }
3520
3521 /* check if the location file is specified for location but NOT for
3522 * distances: GRB case */
3523 if ( dDistr!=distFromSourceFile && lDistr==locationFromSourceFile &&
3524 mwLuminosity>0.0 )
3525 {
3526 fprintf( stderr,
3527 "WARNING: source file specified for locations "
3528 "but NOT for distances, while Milky Way injections "
3529 "are allowed. This might give strange distributions\n" );
3530 }
3531
3532 /* read in the data from the external trigger file */
3533 if ( lDistr == locationFromExttrigFile && !exttrigFileName )
3534 {
3535 fprintf( stderr,
3536 "If --l-distr exttrig is specified, must specify "
3537 "external trigger XML file using --exttrig-file.\n");
3538 exit( 1 );
3539 }
3540 if ( lDistr == locationFromExttrigFile && exttrigFileName )
3541 {
3542 numExtTriggers=LALExtTriggerTableFromLIGOLw( &exttrigHead, exttrigFileName,
3543 0, 1);
3544 fprintf(stderr,
3545 "Number of triggers read from the external trigger file: %d\n",
3546 numExtTriggers);
3547
3548 if (numExtTriggers>1)
3549 {
3550 fprintf(stderr,
3551 "WARNING: Only 1 external trigger expected in the file '%s'",
3552 exttrigFileName );
3553 }
3554 if (numExtTriggers==0)
3555 {
3556 fprintf(stderr,
3557 "ERROR: No external trigger found in file '%s'",
3558 exttrigFileName );
3559 exit(1);
3560 }
3561 }
3562
3563 /* check inclination distribution */
3564 if ( ( iDistr == gaussianInclDist ) && ( inclStd < 0.0 ) )
3565 {
3566 fprintf( stderr,
3567 "Must specify width for gaussian inclination distribution; \n"
3568 "use --incl-std.\n" );
3569 exit( 1 );
3570 }
3571 if ( ( iDistr == fixedInclDist ) && ( fixed_inc < 0. ) )
3572 {
3573 fprintf( stderr,
3574 "Must specify an inclination if you want it fixed; \n"
3575 "use --fixed-inc.\n" );
3576 exit( 1 );
3577 }
3578
3579 /* require --f-lower be explicit */
3580 if ( fLower <= 0.0 )
3581 {
3582 fprintf( stderr, "--f-lower must be specified and non-zero\n" );
3583 exit( 1 );
3584 }
3585
3586 /* check files have been specified correctly for mass distributions */
3587 if ( !massFileName && mDistr==massFromSourceFile )
3588 {
3589 fprintf( stderr,
3590 "Must specify either a file contining the masses (--mass-file) \n"
3591 "or choose another mass-distribution (--m-distr).\n" );
3592 exit( 1 );
3593 }
3594 if ( !nrFileName && mDistr==massFromNRFile )
3595 {
3596 fprintf( stderr,
3597 "Must specify either a file contining the masses (--nr-file) \n"
3598 "or choose another mass-distribution (--m-distr).\n" );
3599 exit( 1 );
3600 }
3601
3602 /* read the masses from the mass file here, check for junk options */
3603 if ( massFileName && mDistr==massFromSourceFile )
3604 {
3605 if ( minMass1>0.0 || minMass2>0.0 || maxMass1>0.0 || maxMass2>0.0 ||
3606 minMtotal>0.0 || maxMtotal>0.0 || meanMass1>0.0 || meanMass2>0.0 ||
3607 massStdev1>0.0 || massStdev2>0.0 || minMassRatio>0.0 || maxMassRatio>0.0 ||
3608 pntMass1>1 || pntMass2>1 || fixedMass1>0.0 || fixedMass2>0.0 )
3609 {
3610 fprintf( stderr,
3611 "One or more mass distribution options are incompatible with \n"
3612 "using a file containing masses (--mass-file).\n" );
3613 exit( 1 );
3614 }
3615 else
3616 {
3617 read_mass_data( massFileName );
3618 }
3619 }
3620
3621 /* NR option requires min- & max-mtotal, all other options are junk */
3622 if ( nrFileName && mDistr==massFromNRFile )
3623 {
3624 if ( minMtotal<=0.0 || maxMtotal<=0.0 )
3625 {
3626 fprintf( stderr,
3627 "Must specify positive min-mtotal and max-mtotal when using \n"
3628 "a file containing NR injection masses (--nr-file).\n" );
3629 exit( 1 );
3630 }
3631 else if ( minMass1>0.0 || minMass2>0.0 || maxMass1>0.0 || maxMass2>0.0 ||
3632 meanMass1>0.0 || meanMass2>0.0 ||
3633 massStdev1>0.0 || massStdev2>0.0 || minMassRatio>0.0 || maxMassRatio>0.0 ||
3634 pntMass1>1 || pntMass2>1 || fixedMass1>0.0 || fixedMass2>0.0 )
3635 {
3636 fprintf( stderr,
3637 "One or more mass distribution options are incompatible with \n"
3638 "using a file containing masses (--nr-file).\n" );
3639 exit( 1 );
3640 }
3641 else
3642 {
3643 read_nr_data ( nrFileName );
3644 }
3645 }
3646
3647 /* inverse logic check */
3648 if ( ( massFileName && mDistr!=massFromSourceFile ) ||
3649 ( nrFileName && mDistr!=massFromNRFile ) )
3650 {
3651 fprintf( stderr,
3652 "Cannot specify a source or NR injection mass file for your choice \n"
3653 "of --m-distr.\n" );
3654 exit( 1 );
3655 }
3656
3657 /* check options for component-based distributions */
3658 if ( mDistr==uniformTotalMass || mDistr==uniformComponentMass ||
3659 mDistr==logComponentMass || mDistr==gaussianMassDist ||
3660 mDistr==m1m2SquareGrid )
3661 {
3662 /* first check: require component mass ranges */
3663 if ( minMass1<=0.0 || minMass2<=0.0 || maxMass1<=0.0 || maxMass2<=0.0 )
3664 {
3665 fprintf( stderr,
3666 "Must specify positive minimum and maximum component masses for \n"
3667 "your choice of --m-distr.\n" );
3668 exit( 1 );
3669 }
3670 /* second check: exclude junk options */
3671 if ( minMassRatio>=0.0 || maxMassRatio>=0.0 ||
3672 fixedMass1>=0.0 || fixedMass2>=0.0 )
3673 {
3674 fprintf( stderr,
3675 "Cannot specify --min-mratio, --max-mratio, or fixed m1,m2 for \n"
3676 "your choice of --m-distr.\n" );
3677 exit( 1 );
3678 }
3679 /* if max-mtotal and min-mtotal values were not specified, assign them */
3680 if ( minMtotal<0.0 )
3681 { minMtotal = minMass1 + minMass2; }
3682 if ( maxMtotal<0.0 )
3683 { maxMtotal = maxMass1 + maxMass2; }
3684 /* third check: proper values of min-mtotal and max-mtotal */
3685 if ( maxMtotal<(minMass1 + minMass2) || minMtotal>(maxMass1 + maxMass2) )
3686 {
3687 fprintf( stderr,
3688 "Maximum (minimum) total mass must be larger (smaller) than \n"
3689 "minMass1+minMass2 (maxMass1+maxMass2). Check your arithmetic\n");
3690 exit( 1 );
3691 }
3692 }
3693
3694 /* check options for mtotal/q-based distributions */
3695 if ( mDistr==uniformTotalMassRatio || mDistr==logMassUniformTotalMassRatio ||
3696 mDistr==uniformTotalMassFraction )
3697 {
3698 /* first check: required options */
3699 if ( minMassRatio<=0.0 || maxMassRatio<=0.0 ||
3700 minMtotal<=0.0 || maxMtotal<=0.0 )
3701 {
3702 fprintf( stderr,
3703 "Must specify positive min-mratio and max-mratio and min/max mtotal for \n"
3704 "your choice of --m-distr.\n" );
3705 exit( 1 );
3706 }
3707 /* second check: exclude junk options */
3708 if ( minMass1>=0.0 || minMass2>=0.0 || maxMass1>=0.0 || maxMass2>=0.0 ||
3709 meanMass1>=0.0 || meanMass2>=0.0 || massStdev1>=0.0 || massStdev2>=0.0 ||
3710 pntMass1>1 || pntMass2>1 || fixedMass1>=0.0 || fixedMass2>=0.0 )
3711 {
3712 fprintf( stderr,
3713 "Cannot specify options related to component masses for your choice \n"
3714 "of --m-distr.\n" );
3715 exit( 1 );
3716 }
3717 }
3718
3719 /* check for gaussian mass distribution parameters */
3720 if ( mDistr==gaussianMassDist )
3721 {
3722 if ( meanMass1 <= 0.0 || massStdev1 <= 0.0 ||
3723 meanMass2 <= 0.0 || massStdev2 <= 0.0 )
3724 {
3725 fprintf( stderr,
3726 "Must specify positive --mean-mass1/2 and --stdev-mass1/2 if choosing \n"
3727 " --m-distr gaussian\n" );
3728 exit( 1 );
3729 }
3730 if ( minMass1==maxMass1 || minMass2==maxMass2 )
3731 {
3732 fprintf( stderr,
3733 "Must specify a nonzero range of mass1 and mass2 if choosing \n"
3734 " --m-distr gaussian\n" );
3735 exit( 1 );
3736 }
3737 }
3738 /* inverse logic check for junk options */
3739 if ( ( meanMass1>=0.0 || meanMass2>=0.0 || massStdev1>=0.0 || massStdev2>=0.0 )
3740 && ( mDistr!=gaussianMassDist ) )
3741 {
3742 fprintf( stderr,
3743 "Cannot specify --mean-mass1/2 or --stdev-mass1/2 unless choosing \n"
3744 " --m-distr gaussian\n" );
3745 exit( 1 );
3746 }
3747
3748 /* checks for m1m2 mass grid options */
3749 if ( mDistr==m1m2SquareGrid )
3750 {
3751 if ( pntMass1<2 || pntMass2<2 )
3752 {
3753 fprintf( stderr,
3754 "Values of --mass1-points and --mass2-points must be specified \n"
3755 "and >= 2 if choosing --m-distr m1m2SquareGrid\n" );
3756 exit( 1 );
3757 }
3758 else
3759 {
3760 deltaMass1 = ( maxMass1 - minMass1 ) / (REAL4) ( pntMass1 -1 );
3761 deltaMass2 = ( maxMass2 - minMass2 ) / (REAL4) ( pntMass2 -1 );
3762 }
3763 }
3764 /* inverse logic check */
3765 if ( ( pntMass1>1 || pntMass2>1 ) && ( mDistr!=m1m2SquareGrid ) )
3766 {
3767 fprintf( stderr,
3768 "Cannot specify --mass1-points or mass2-points unless choosing \n"
3769 " --m-distr m1m2SquareGrid\n" );
3770 exit( 1 );
3771 }
3772
3773 /* checks for fixed-mass injections */
3774 if ( mDistr==fixMasses )
3775 {
3776 if ( fixedMass1<=0.0 || fixedMass2<=0.0 )
3777 {
3778 fprintf( stderr, "--fixed-mass1 and --fixed-mass2 must be specified "
3779 "and positive if choosing --m-distr fixMasses\n" );
3780 exit( 1 );
3781 }
3782 /* exclude junk options */
3783 if ( minMass1>=0.0 || minMass2>=0.0 || maxMass1>=0.0 || maxMass2>=0.0 ||
3784 minMtotal>=0.0 || maxMtotal>=0.0 || meanMass1>=0.0 || meanMass2>=0.0 ||
3785 massStdev1>=0.0 || massStdev2>=0.0 || minMassRatio>=0.0 ||
3786 maxMassRatio>=0.0 || pntMass1>1 || pntMass2>1 )
3787 {
3788 fprintf( stderr,
3789 "One or more mass options are incompatible with --m-distr fixMasses, \n"
3790 "only --fixed-mass1 and --fixed-mass2 may be specified\n" );
3791 exit( 1 );
3792 }
3793 }
3794
3795 /* check if waveform is specified */
3796 if ( !*waveform )
3797 {
3798 fprintf( stderr, "No waveform specified (--waveform).\n" );
3799 exit( 1 );
3800 }
3801
3802 if ( spinInjections==-1 && mDistr != massFromNRFile )
3803 {
3804 fprintf( stderr,
3805 "Must specify --disable-spin or --enable-spin\n"
3806 "unless doing NR injections\n" );
3807 exit( 1 );
3808 }
3809
3810 if ( spinInjections==0 && spinAligned==1 )
3811 {
3812 fprintf( stderr,
3813 "Must enable spin to obtain aligned spin injections.\n" );
3814 exit( 1 );
3815 }
3816
3817 if ( spinInjections==1 && spinAligned==-1 &&
3818 ( !strncmp(waveform, "IMRPhenomB", 10) || !strncmp(waveform, "IMRPhenomC", 10) ) )
3819 {
3820 fprintf( stderr,
3821 "Spinning IMRPhenomB or -C injections must have the --aligned option.\n" );
3822 exit( 1 );
3823 }
3824
3825 if ( spinInjections==1 )
3826 {
3827 if ( spinDistr == unknownSpinDist ) /* Not currently used */
3828 {
3829 fprintf(stderr,"Must specify a spin magnitude distribution (--spin-distr).\n");
3830 exit( 1 );
3831 }
3832
3833 /* check that spin stddev is positive */
3834 if ( spinDistr==gaussianSpinDist && (Spin1Std <= 0.0 || Spin2Std <= 0.0))
3835 {
3836 fprintf( stderr,
3837 "Must specify positive |spin| standard deviations when using"
3838 " --spin-gaussian\n" );
3839 exit( 1 );
3840 }
3841
3842 /* check that spins are in range 0 - 1 */
3843 if (minSpin1 < 0. || minSpin2 < 0. || maxSpin1 > 1. || maxSpin2 >1.)
3844 {
3845 fprintf( stderr,
3846 "Spins can only take values between 0 and 1.\n" );
3847 exit( 1 );
3848 }
3849
3850 /* check max and mins are the correct way around */
3851 if (minSpin1 > maxSpin1 || minSpin2 > maxSpin2 )
3852 {
3853 fprintf( stderr,
3854 "Minimal spins must be less than maximal spins.\n" );
3855 exit( 1 );
3856 }
3857
3858 /* check that spin means are within a reasonable range */
3859 if (spinDistr==gaussianSpinDist && ( minSpin1 - meanSpin1 > 2.0*Spin1Std || meanSpin1 - maxSpin1 > 2.0*Spin2Std ))
3860 {
3861 fprintf(stderr,"Mean of |spin1| distribution is way out of range.\n");
3862 exit( 1 );
3863 }
3864 if (spinDistr==gaussianSpinDist && ( minSpin2 - meanSpin2 > 2.0*Spin2Std || meanSpin2 - maxSpin2 > 2.0*Spin2Std ))
3865 {
3866 fprintf(stderr,"Mean of |spin2| distribution is way out of range.\n");
3867 exit( 1 );
3868 }
3869
3870 /* check that selection criteria for kappa are unique */
3871 if ( (minKappa1 > -1.0 || maxKappa1 < 1.0) &&
3872 (minabsKappa1 > 0.0 || maxabsKappa1 < 1.0) )
3873 {
3874 fprintf( stderr,
3875 "Either the options [--min-kappa1,--max-kappa1] or\n"
3876 "[--min-abskappa1,--max-abskappa1] can be specified\n" );
3877 exit( 1 );
3878 }
3879
3880 /* check that kappa is in range */
3881 if (minKappa1 < -1.0 || maxKappa1 > 1.0)
3882 {
3883 fprintf( stderr,
3884 "Kappa can only take values between -1 and +1\n" );
3885 exit( 1 );
3886 }
3887
3888 /* check that kappa min-max are set correctly */
3889 if (minKappa1 > maxKappa1)
3890 {
3891 fprintf( stderr,
3892 "Minimal kappa must be less than maximal kappa\n" );
3893 exit( 1 );
3894 }
3895
3896 /* check that abskappa is in range */
3897 if (minabsKappa1 < 0.0 || maxabsKappa1 > 1.0)
3898 {
3899 fprintf( stderr,
3900 "The absolute value of kappa can only take values between 0 and +1\n" );
3901 exit( 1 );
3902 }
3903
3904 /* check that kappa min-max are set correctly */
3905 if (minabsKappa1 > maxabsKappa1)
3906 {
3907 fprintf( stderr,
3908 "Minimal kappa must be less than maximal kappa\n" );
3909 exit( 1 );
3910 }
3911 }
3912
3913 if ( dDistr == starFormationRate )
3914 {
3915 /* recalculate mean time step from the SFR */
3916 meanTimeStep = mean_time_step_sfr(maxZ, localRate);
3917 }
3918
3919 if (meanTimeStep<=0 && tDistr != LALINSPIRAL_FILE_TIME_DIST)
3920 {
3921 fprintf( stderr,
3922 "Minimum time step value must be larger than zero\n" );
3923 exit( 1 );
3924 }
3925
3926 if (!injtimesFileName && tDistr == LALINSPIRAL_FILE_TIME_DIST)
3927 {
3928 fprintf(stderr, "No filename for injection GPStimes is given. Use --time-file.\n");
3929 }
3930
3931 if ( injtimesFileName && tDistr != LALINSPIRAL_FILE_TIME_DIST )
3932 {
3933 fprintf( stderr,
3934 "Cannot specify an injection times file for your choice of --t-distr.\n" );
3935 exit( 1 );
3936 }
3937
3938 if (timeInterval > 0. && (tDistr == LALINSPIRAL_EXPONENTIAL_TIME_DIST || tDistr == LALINSPIRAL_FILE_TIME_DIST) )
3939 {
3940 fprintf( stderr,
3941 "time interval must be zero\n" );
3942 exit( 1 );
3943 }
3944
3945 if ( injtimesFileName && tDistr == LALINSPIRAL_FILE_TIME_DIST)
3946 {
3947 if (meanTimeStep > 0.)
3948 {
3949 fprintf(stderr, "Minimum time step value must be larger than zero\n" );
3950 exit(1);
3951 }
3952 // printf("Reading injection times from file %s\n", injtimesFileName);
3953 read_time_data(injtimesFileName);
3954 }
3955
3956 if ( userTag && outCompress )
3957 {
3958 snprintf( fname, sizeof(fname), "HL-INJECTIONS_%d_%s-%d-%ld.xml.gz",
3959 rand_seed, userTag, gpsStartTime.gpsSeconds, gpsDuration );
3960 }
3961 else if ( userTag && !outCompress )
3962 {
3963 snprintf( fname, sizeof(fname), "HL-INJECTIONS_%d_%s-%d-%ld.xml",
3964 rand_seed, userTag, gpsStartTime.gpsSeconds, gpsDuration );
3965 }
3966 else if ( !userTag && outCompress )
3967 {
3968 snprintf( fname, sizeof(fname), "HL-INJECTIONS_%d-%d-%ld.xml.gz",
3969 rand_seed, gpsStartTime.gpsSeconds, gpsDuration );
3970 }
3971 else
3972 {
3973 snprintf( fname, sizeof(fname), "HL-INJECTIONS_%d-%d-%ld.xml",
3974 rand_seed, gpsStartTime.gpsSeconds, gpsDuration );
3975 }
3976 if ( outputFileName )
3977 {
3978 snprintf( fname, sizeof(fname), "%s",
3979 outputFileName);
3980 }
3981
3982 /* increment the random seed by the GPS start time:*/
3983 rand_seed += gpsStartTime.gpsSeconds;
3984
3985 /* set up the LAL random number generator */
3986 LALCreateRandomParams( &status, &randParams, rand_seed );
3987
3988 this_proc_param = procparams;
3989 procparams = procparams->next;
3990 free( this_proc_param );
3991
3992 /* create the first injection */
3993 simTable = injections = (SimInspiralTable *)
3994 calloc( 1, sizeof(SimInspiralTable) );
3995
3996 simRingTable = ringparams = (SimRingdownTable *)
3997 calloc( 1, sizeof(SimRingdownTable) );
3998
3999 /* set redshift to zero */
4000 redshift=0.;
4001
4002 /* set mass distribution parameters to their value at z = 0 */
4003 minMass10 = minMass1;
4004 maxMass10 = maxMass1;
4005 minMass20 = minMass2;
4006 maxMass20 = maxMass2;
4007 minMtotal0 = minMtotal;
4008 maxMtotal0 = maxMtotal;
4009 meanMass10 = meanMass1;
4010 meanMass20 = meanMass2;
4011 massStdev10 = massStdev1;
4012 massStdev20 = massStdev2;
4013
4014 /* calculate the maximal value of the probability distribution of the redshift */
4015 if (dDistr == starFormationRate)
4016 {
4017 pzmax = probability_redshift(maxZ);
4018 }
4019
4020 /* loop over parameter generation until end time is reached */
4021 ninj = 0;
4022 ncount = 0;
4023 currentGpsTime = gpsStartTime;
4024 if (tDistr == LALINSPIRAL_FILE_TIME_DIST){
4025 currentGpsTime.gpsSeconds = inj_times[0].gpsSeconds;
4026 currentGpsTime.gpsNanoSeconds = inj_times[0].gpsNanoSeconds;
4027 }
4028
4029 while ( 1 )
4030 {
4031 /* increase counter */
4032 ninj++;
4033
4034 /* store time in table */
4035 simTable=XLALRandomInspiralTime( simTable, randParams,
4036 currentGpsTime, timeInterval );
4037
4038 /* populate waveform and other parameters */
4039 memcpy( simTable->waveform, waveform,
4040 sizeof(CHAR) * LIGOMETA_WAVEFORM_MAX );
4041 simTable->f_lower = fLower;
4042 simTable->amp_order = amp_order;
4043
4044 /* draw redshift and apply to mass parameters */
4045 if (dDistr==starFormationRate)
4046 {
4047 redshift = drawRedshift(minZ,maxZ,pzmax);
4048
4049 minMass1 = redshift_mass(minMass10, redshift);
4050 maxMass1 = redshift_mass(maxMass10, redshift);
4051 meanMass1 = redshift_mass(meanMass10, redshift);
4052 massStdev1 = redshift_mass(massStdev10, redshift);
4053 minMass2 = redshift_mass(minMass20, redshift);
4054 maxMass2 = redshift_mass(maxMass20, redshift);
4055 meanMass2 = redshift_mass(meanMass20, redshift);
4056 massStdev2 = redshift_mass(massStdev20, redshift);
4057 minMtotal = redshift_mass(minMtotal0, redshift);
4058 maxMtotal = redshift_mass(maxMtotal0, redshift);
4059 }
4060
4061 /* populate masses */
4062 if ( mDistr==massFromSourceFile )
4063 {
4064 drawMassFromSource( simTable );
4065 }
4066 else if ( mDistr==massFromNRFile )
4067 {
4068 if (ninjaMass)
4069 drawMassSpinFromNRNinja2( simTable );
4070 else
4071 drawMassSpinFromNR( simTable );
4072 }
4073 else if ( mDistr==gaussianMassDist )
4074 {
4075 simTable=XLALGaussianInspiralMasses( simTable, randParams,
4076 minMass1, maxMass1,
4077 meanMass1, massStdev1,
4078 minMass2, maxMass2,
4079 meanMass2, massStdev2);
4080 }
4081 else if ( mDistr==uniformTotalMassRatio )
4082 {
4083 simTable=XLALRandomInspiralTotalMassRatio(simTable, randParams,
4084 mDistr, minMtotal, maxMtotal, minMassRatio, maxMassRatio );
4085 }
4086 else if ( mDistr==logMassUniformTotalMassRatio )
4087 {
4088 simTable=XLALRandomInspiralTotalMassRatio(simTable, randParams,
4089 mDistr, minMtotal, maxMtotal, minMassRatio, maxMassRatio );
4090 }
4091 else if ( mDistr==m1m2SquareGrid )
4092 {
4093 simTable=XLALm1m2SquareGridInspiralMasses( simTable, minMass1, minMass2,
4094 minMtotal, maxMtotal, deltaMass1, deltaMass2, pntMass1, pntMass2,
4095 ninj, &ncount);
4096 }
4097 else if ( mDistr==fixMasses )
4098 {
4099 simTable=XLALFixedInspiralMasses( simTable, fixedMass1, fixedMass2);
4100 }
4101 else if ( mDistr==uniformTotalMassFraction )
4102 {
4103 simTable=XLALRandomInspiralTotalMassFraction(simTable, randParams,
4104 mDistr, minMtotal, maxMtotal, minMassRatio, maxMassRatio );
4105 }
4106 else {
4107 simTable=XLALRandomInspiralMasses( simTable, randParams, mDistr,
4108 minMass1, maxMass1,
4109 minMass2, maxMass2,
4110 minMtotal, maxMtotal);
4111 }
4112
4113 /* draw location and distances */
4114 drawFromSource( &drawnRightAscension, &drawnDeclination, &drawnDistance,
4115 drawnSourceName );
4116 drawFromIPNsim( &drawnRightAscension, &drawnDeclination );
4117
4118 /* populate distances */
4119 if ( dDistr == distFromSourceFile )
4120 {
4121 if ( maxD > 0 )
4122 {
4123 while ( drawnDistance > maxD/1000.0 )
4124 {
4125 drawFromSource( &drawnRightAscension, &drawnDeclination,
4126 &drawnDistance, drawnSourceName );
4127 }
4128 }
4129 simTable->distance = drawnDistance;
4130 }
4131 else if ( dDistr == starFormationRate )
4132 {
4133 /* fit of luminosity distance between z=0-1, in Mpc for h0=0.7, omega_m=0.3, omega_v=0.7*/
4134 simTable->distance = luminosity_distance(redshift);
4135 }
4136 else if (dDistr== uniformDistance || dDistr== uniformDistanceSquared || dDistr== uniformLogDistance || dDistr==uniformVolume)
4137 {
4138 simTable=XLALRandomInspiralDistance(simTable, randParams,
4139 dDistr, minD/1000.0, maxD/1000.0);
4140 }
4141 else if (dDistr==uniformSnr || dDistr==uniformLogSnr || dDistr==uniformVolumeSnr)
4142 {
4143 /* Set distance to just any value, e.g. 100, which will be used to scale the SNR */
4144 simTable->distance=100.0;
4145 }
4146 /* Possible errors (i.e. no dDistr given) should have already been caught */
4147 /* check just to be sure in case someone adds new LoudnessDistribution */
4148 else
4149 {
4150 fprintf(stderr,"Error while generating the distance of the event.\n");
4151 exit(1);
4152 }
4153 /* Scale by chirp mass if desired, relative to a 1.4,1.4 object */
4154 if ( useChirpDist )
4155 {
4156 REAL4 scaleFac;
4157 scaleFac = simTable->mchirp/(2.8*pow(0.25,0.6));
4158 simTable->distance = simTable->distance*pow(scaleFac,5./6.);
4159 }
4160
4161 /* populate location */
4162 if ( lDistr == locationFromSourceFile )
4163 {
4164 simTable->longitude = drawnRightAscension;
4165 simTable->latitude = drawnDeclination;
4166 memcpy( simTable->source, drawnSourceName,
4167 sizeof(CHAR) * LIGOMETA_SOURCE_MAX );
4168 }
4169 else if ( lDistr == locationFromExttrigFile )
4170 {
4171 drawLocationFromExttrig( simTable );
4172 }
4173 else if ( lDistr == fixedSkyLocation)
4174 {
4175 simTable->longitude = longitude;
4176 simTable->latitude = latitude;
4177 }
4178 else if ( lDistr == uniformSkyLocation )
4179 {
4180 simTable=XLALRandomInspiralSkyLocation(simTable, randParams);
4181 }
4182 else if ( lDistr == locationFromIPNFile )
4183 {
4184 IPNgmst1 = XLALGreenwichMeanSiderealTime(&IPNgpsTime);
4185 IPNgmst2 = XLALGreenwichMeanSiderealTime(&simTable->geocent_end_time);
4186 simTable->longitude = drawnRightAscension - IPNgmst1 + IPNgmst2;
4187 simTable->latitude = drawnDeclination;
4188 }
4189 else
4190 {
4191 fprintf( stderr,
4192 "Unknown location distribution specified. Possible choices: "
4193 "source, exttrig, random or fixed\n" );
4194 exit( 1 );
4195 }
4196
4197 /* populate polarization, inclination, and coa_phase */
4198 do
4199 {
4200 simTable=XLALRandomInspiralOrientation(simTable, randParams,
4201 iDistr, inclStd);
4202 } while ( (fabs(cos(simTable->inclination))<cos(max_inc)) );
4203
4204 /* override inclination */
4205 if ( iDistr == fixedInclDist )
4206 {
4207 simTable->inclination = fixed_inc;
4208 }
4209
4210 /* override polarization angle */
4211 if ( psi != -1.0 )
4212 {
4213 simTable->polarization = psi;
4214 }
4215
4216 /* override coalescence phase */
4217 if ( coaPhaseFixed )
4218 {
4219 simTable->coa_phase = fixedCoaPhase;
4220 }
4221
4222 /* populate spins, if required */
4223 if (spinInjections==1) {
4224 if ( spinAligned==1 ) {
4225 if ( !strncmp(axisChoiceString,"angmomentum", 11) )
4226 alignInj = alongzAxis;
4227 else {
4228 if ( !strncmp(axisChoiceString,"view", 4) )
4229 alignInj = inxzPlane;
4230 else {
4231 fprintf( stderr, "Unknown axis-choice specification: 'angmomentum' and 'view' allowed, %s given.\n",axisChoiceString);
4232 exit( 1 );
4233 }
4234 }
4235 }
4236 else
4237 alignInj = notAligned;
4238
4239 simTable = XLALRandomInspiralSpins( simTable, randParams,
4240 minSpin1, maxSpin1,
4241 minSpin2, maxSpin2,
4242 minKappa1, maxKappa1,
4243 minabsKappa1, maxabsKappa1,
4244 alignInj, spinDistr,
4245 meanSpin1, Spin1Std,
4246 meanSpin2, Spin2Std );
4247 }
4248
4249 /* adjust SNR to desired distribution using NINJA calculation */
4250 if ( ifos != NULL && ninjaSNR )
4251 {
4252 if (dDistr==uniformSnr){
4253 targetSNR=draw_uniform_snr(minSNR,maxSNR);
4254 }
4255 else if(dDistr==uniformLogSnr){
4256 targetSNR=draw_log10_snr(minSNR,maxSNR);
4257 }
4258 else if (dDistr==uniformVolumeSnr){
4259 targetSNR=draw_volume_snr(minSNR,maxSNR);
4260 }
4261 else{
4262 fprintf(stderr,"Allowed values for --snr-distr are uniform, log10 and volume. Exiting...\n");
4263 exit(1);
4264 }
4265
4266 if (! real8Ninja2)
4267 {
4268 adjust_snr(simTable, targetSNR, ifos);
4269 }
4270 else
4271 {
4272 REAL8 *start_freqs;
4273 const char **ifo_list;
4274 REAL8FrequencySeries **psds;
4275 int count, num_ifos = 0;
4276 char *tmp, *ifo;
4277
4278 tmp = LALCalloc(1, strlen(ifos) + 1);
4279 strcpy(tmp, ifos);
4280 ifo = strtok (tmp,",");
4281
4282 while (ifo != NULL)
4283 {
4284 num_ifos += 1;
4285 ifo = strtok (NULL, ",");
4286 }
4287
4288 start_freqs = (REAL8 *) LALCalloc(num_ifos, sizeof(REAL8));
4289 ifo_list = (const char **) LALCalloc(num_ifos, sizeof(char *));
4290 psds = (REAL8FrequencySeries **) LALCalloc(num_ifos, sizeof(REAL8FrequencySeries *));
4291
4292 strcpy(tmp, ifos);
4293 ifo = strtok (tmp,",");
4294 count = 0;
4295
4296 while (ifo != NULL)
4297 {
4298 ifo_list[count] = ifo;
4299
4300 if (ifo_list[count][0] == 'V')
4301 {
4302 start_freqs[count] = virgoStartFreq;
4303 psds[count] = virgoPsd;
4304 }
4305 else
4306 {
4307 start_freqs[count] = ligoStartFreq;
4308 psds[count] = ligoPsd;
4309 }
4310 count++;
4311 ifo = strtok (NULL, ",");
4312 }
4313
4314 adjust_snr_with_psds_real8(simTable, targetSNR, num_ifos, ifo_list, psds, start_freqs);
4315
4316 LALFree(start_freqs);
4317 LALFree(ifo_list);
4318 LALFree(psds);
4319 LALFree(tmp);
4320 }
4321 }
4322
4323 /* adjust SNR to desired distribution using LALSimulation WF Generator */
4324 if (ifos!=NULL && !ninjaSNR)
4325 {
4326 char *ifo;
4327 REAL8 *start_freqs;
4328 REAL8FrequencySeries **psds;
4329 i=1;
4330
4331 /*reset counter */
4332 ifo = ifonames[0];
4333 i = 0;
4334 /* Create variables for PSDs and starting frequencies */
4335 start_freqs = (REAL8 *) LALCalloc(numifos+1, sizeof(REAL8));
4336 psds = (REAL8FrequencySeries **) LALCalloc(numifos+1, sizeof(REAL8FrequencySeries *));
4337
4338 /* Hardcoded values of srate and segment length. If changed here they must also be changed in inspiralutils.c/calculate_lalsim_snr */
4339 REAL8 srate = 4096.0;
4340 /* Increase srate for EOB WFs */
4341 const char* WF=simTable->waveform;
4342 if (strstr(WF,"EOB"))
4343 srate = 8192.0;
4344 /* We may want to increase the segment length when starting at low frequencies */
4345 REAL8 segment = 64.0;
4346 size_t seglen = (size_t) segment*srate;
4347
4348 /* Fill psds and start_freqs */
4349 /* If the user did not provide files for the PSDs, use XLALSimNoisePSD to fill in ligoPsd and virgoPsd */
4350 while(ifo !=NULL)
4351 {
4352 if(!strcmp("V1",ifo))
4353 {
4354 start_freqs[i]=virgoStartFreq;
4355 if (!virgoPsd)
4356 {
4357 virgoPsd=XLALCreateREAL8FrequencySeries("VPSD", &(simTable->geocent_end_time), 0, 1.0/segment, &lalHertzUnit, seglen/2+1);
4358 get_FakePsdFromString(virgoPsd,virgoFakePsd, virgoStartFreq);
4359 }
4360 if (!virgoPsd) fprintf(stderr,"Failed to produce Virgo PSD series. Exiting...\n");
4361 psds[i]=virgoPsd;
4362 }
4363 else if (!strcmp("L1",ifo) || !strcmp("H1",ifo))
4364 {
4365 start_freqs[i]=ligoStartFreq;
4366 if (!ligoPsd)
4367 {
4368 ligoPsd=XLALCreateREAL8FrequencySeries("LPSD", &(simTable->geocent_end_time), 0, 1.0/segment, &lalHertzUnit, seglen/2+1);
4369 get_FakePsdFromString(ligoPsd,ligoFakePsd,ligoStartFreq);
4370 }
4371 if (!ligoPsd) fprintf(stderr,"Failed to produce LIGO PSD series. Exiting...\n");
4372 psds[i]=ligoPsd;
4373 }
4374 else
4375 {
4376 fprintf(stderr,"Unknown IFO. Allowed IFOs are H1,L1 and V1. Exiting...\n");
4377 exit(-1);
4378 }
4379 i++;
4380 ifo=ifonames[i];
4381 }
4382
4383 /* If exactly one detector is specified, turn the single IFO snr check off */
4384 if (numifos<2)
4385 {
4386 fprintf(stdout,"Warning: You are using less than 2 IFOs. Disabling the single IFO SNR threshold check...\n");
4387 single_IFO_SNR_threshold=0.0;
4388 }
4389
4390 /* This function draws a proposed SNR and rescales the distance accordingly */
4391 scale_lalsim_distance(simTable, ifonames, psds, start_freqs, dDistr);
4392
4393 /* Clean */
4394 if (psds) LALFree(psds);
4395 if (start_freqs) LALFree(start_freqs);
4396 }
4397
4398 /* populate the site specific information: end times and effective distances */
4399 LALPopulateSimInspiralSiteInfo( &status, simTable );
4400
4401 /* populate the taper options */
4402 {
4403 switch (taperInj)
4404 {
4405 case LAL_SIM_INSPIRAL_TAPER_NONE:
4406 snprintf( simTable->taper, LIGOMETA_INSPIRALTAPER_MAX,
4407 "%s", "TAPER_NONE");
4408 break;
4409 case LAL_SIM_INSPIRAL_TAPER_START:
4410 snprintf( simTable->taper, LIGOMETA_INSPIRALTAPER_MAX,
4411 "%s", "TAPER_START");
4412 break;
4413 case LAL_SIM_INSPIRAL_TAPER_END:
4414 snprintf( simTable->taper, LIGOMETA_INSPIRALTAPER_MAX,
4415 "%s", "TAPER_END");
4416 break;
4417 case LAL_SIM_INSPIRAL_TAPER_STARTEND:
4418 snprintf( simTable->taper, LIGOMETA_INSPIRALTAPER_MAX,
4419 "%s", "TAPER_STARTEND");
4420 break;
4421 default: /* Never reach here */
4422 fprintf( stderr, "unknown error while populating sim_inspiral taper options\n" );
4423 exit(1);
4424 }
4425
4426 }
4427
4428 /* populate the bandpass options */
4429 simTable->bandpass = bandPassInj;
4430
4431
4432 /* increment current time, avoiding roundoff error;
4433 check if end of loop is reached */
4434 if (tDistr == LALINSPIRAL_EXPONENTIAL_TIME_DIST)
4435 {
4436 XLALGPSAdd( &currentGpsTime, -(REAL8)meanTimeStep * log( XLALUniformDeviate(randParams) ) );
4437 }
4438 else if (tDistr == LALINSPIRAL_FILE_TIME_DIST)
4439 {
4440 if (ninj >= (size_t) n_times)
4441 break;
4442 currentGpsTime.gpsSeconds = inj_times[ninj].gpsSeconds;
4443 currentGpsTime.gpsNanoSeconds = inj_times[ninj].gpsNanoSeconds;
4444 }
4445 else
4446 {
4447 currentGpsTime = gpsStartTime;
4448 XLALGPSAdd( &currentGpsTime, ninj * meanTimeStep );
4449 }
4450 if ( XLALGPSCmp( &currentGpsTime, &gpsEndTime ) >= 0 && tDistr!=LALINSPIRAL_FILE_TIME_DIST )
4451 break;
4452
4453 /* allocate and go to next SimInspiralTable */
4454 simTable = simTable->next = (SimInspiralTable *)
4455 calloc( 1, sizeof(SimInspiralTable) );
4456 simRingTable = simRingTable->next = (SimRingdownTable *)
4457 calloc( 1, sizeof(SimRingdownTable) );
4458
4459 }
4460
4461 /* destroy the structure containing the random params */
4462 LAL_CALL( LALDestroyRandomParams( &status, &randParams ), &status);
4463
4464 /* If we read from an external trigger file, free our external trigger.
4465 exttrigHead is guaranteed to have no children to free. */
4466 if ( exttrigHead != NULL ) {
4467 LALFree(exttrigHead);
4468 }
4469
4470 /* destroy the NR data */
4471 if ( num_nr )
4472 {
4473 for( i = 0; i < num_nr; i++ )
4474 {
4475 LALFree( nrSimArray[i] );
4476 }
4477 LALFree( nrSimArray );
4478 }
4479
4480 xmlfp = XLALOpenLIGOLwXMLFile( fname );
4481 if (!xmlfp) XLAL_ERROR(XLAL_EIO);
4482
4483 XLALGPSTimeNow(&(proctable->end_time));
4484
4485 int retcode = XLALWriteLIGOLwXMLProcessTable(xmlfp, proctable);
4486 if (retcode != XLAL_SUCCESS)
4487 {
4488 XLAL_ERROR(retcode);
4489 }
4490
4491 if ( procparams )
4492 {
4493 retcode = XLALWriteLIGOLwXMLProcessParamsTable(xmlfp, procparams);
4494 if (retcode != XLAL_SUCCESS)
4495 {
4496 XLAL_ERROR(retcode);
4497 }
4498 }
4499
4500 XLALSimInspiralAssignIDs ( injections, 0, 0 );
4501 if ( injections )
4502 {
4503 retcode = XLALWriteLIGOLwXMLSimInspiralTable(xmlfp, injections);
4504 if ( retcode != XLAL_SUCCESS )
4505 {
4506 XLAL_ERROR(retcode);
4507 }
4508 }
4509
4510 retcode = XLALCloseLIGOLwXMLFile ( xmlfp );
4511 if (retcode != XLAL_SUCCESS)
4512 {
4513 XLAL_ERROR(retcode);
4514 }
4515
4516
4517 if (source_data)
4518 LALFree(source_data);
4519 if (mass_data)
4520 LALFree(mass_data);
4521 if (skyPoints)
4522 LALFree(skyPoints);
4523
4524 if ( ligoPsd )
4525 XLALDestroyREAL8FrequencySeries( ligoPsd );
4526 if ( virgoPsd )
4527 XLALDestroyREAL8FrequencySeries( virgoPsd );
4528 if (ifonames) LALFree(ifonames);
4529
4530 LALCheckMemoryLeaks();
4531 return 0;
4532}
4533
4534static void scale_lalsim_distance(SimInspiralTable *inj,
4535 char **IFOnames,
4536 REAL8FrequencySeries **psds,
4537 REAL8 *start_freqs,
4538 LoudnessDistribution snrDistr)
4539{
4540 REAL8 proposedSNR=0.0;
4541 REAL8 local_min=0.0;
4542 REAL8 net_snr=0.0;
4543 REAL8 *SNRs=NULL;
4544 UINT4 above_threshold=0;
4545 REAL8 ratio=1.0;
4546 UINT4 j=0;
4547 UINT4 num_ifos=0;
4548 /* If not already done, set distance to 100Mpc, just to have something while calculating the actual SNR */
4549 if (inj->distance<=0)
4550 inj->distance=100.0;
4551
4552 if (IFOnames ==NULL)
4553 {
4554 fprintf(stderr,"scale_lalsim_distance() called with IFOnames=NULL. Exiting...\n");
4555 exit(1);
4556 }
4557 char * ifo=IFOnames[0];
4558
4559 /* Get the number of IFOs from IFOnames*/
4560 while(ifo !=NULL)
4561 {
4562 num_ifos++;
4563 ifo=IFOnames[num_ifos];
4564 }
4565
4566 SNRs=calloc(num_ifos+1 ,sizeof(REAL8));
4567 /* Calculate the single IFO and network SNR for the dummy distance of 100Mpc */
4568 for (j=0;j<num_ifos;j++)
4569 {
4570 SNRs[j]=calculate_lalsim_snr(inj,IFOnames[j],psds[j],start_freqs[j]);
4571 net_snr+=SNRs[j]*SNRs[j];
4572 }
4573 net_snr=sqrt(net_snr);
4574
4575 local_min=minSNR;
4576 /* Draw a proposed new SNR. Check that two or more IFOs are */
4577 /* above threshold (if given and if num_ifos>=2) */
4578 do
4579 {
4580 above_threshold=num_ifos;
4581 /* Generate a new SNR from given distribution */
4582 if (snrDistr==uniformSnr)
4583 {
4584 proposedSNR=draw_uniform_snr(local_min,maxSNR);
4585 }
4586 else if(snrDistr==uniformLogSnr)
4587 {
4588 proposedSNR=draw_log10_snr(local_min,maxSNR);
4589 }
4590 else if (snrDistr==uniformVolumeSnr)
4591 {
4592 proposedSNR=draw_volume_snr(local_min,maxSNR);
4593 }
4594 else
4595 {
4596 fprintf(stderr,"Allowed values for snr-distr are uniform, uniformLogSnr and uniformVolumeSnr. Exiting...\n");
4597 exit(1);
4598 }
4599
4600 if (vrbflg) { printf("proposed SNR %lf. Proposed new dist %lf \n",
4601 proposedSNR,inj->distance*net_snr/proposedSNR); }
4602
4603 ratio=net_snr/proposedSNR;
4604
4605 /* Check that single ifo SNRs above threshold in at least two IFOs */
4606 for (j=0;j<num_ifos;j++)
4607 {
4608 if (SNRs[j]<single_IFO_SNR_threshold*ratio)
4609 above_threshold--;
4610 }
4611 /* Set the min to the proposed SNR, so that next drawing for */
4612 /* this event (if necessary) will give higher SNR */
4613 local_min=proposedSNR;
4614
4615 /* We hit the upper bound of the Network SNR. It is simply not possible to */
4616 /* have >2 IFOs with single IFO above threshold without getting the network */
4617 /* SNR above maxSNR. */
4618 /* Use the last proposed value (~maxSNR) and continue */
4619 if (maxSNR-proposedSNR<0.1 && single_IFO_SNR_threshold>0.0)
4620 {
4621 fprintf(stdout,"WARNING: Could not get two or more IFOs having SNR>%.1f without\n"
4622 "making the network SNR larger that its maximum value %.1f. Setting SNR to %lf.\n",
4623 single_IFO_SNR_threshold,maxSNR,proposedSNR);
4624
4625 /* set above_threshold to 3 to go out */
4626 above_threshold=3;
4627 }
4628 else if (vrbflg)
4629 {
4630 if (above_threshold<2 && num_ifos>=2)
4631 fprintf(stdout,"WARNING: Proposed SNR does not get two or more IFOs having SNR>%.1f. Re-drawing... \n",single_IFO_SNR_threshold);
4632 }
4633
4634 } while(!(above_threshold>=2) && num_ifos>=2);
4635
4636 inj->distance=inj->distance*ratio;
4637
4638 /* We may want the resulting SNR printed to a file. If so, uncomment those lines
4639 char SnrName[200];
4640 sprintf(SnrName,"SNRfile.txt");
4641 FILE * snrout = fopen(SnrName,"a");
4642 fprintf(snrout,"%lf\n",proposedSNR);
4643 fclose(snrout);
4644 */
4645
4646 if (SNRs) free(SNRs);
4647
4648}
4649
4650static REAL8 draw_volume_snr(REAL8 minsnr,REAL8 maxsnr)
4651{
4652 REAL8 proposedSNR=0.0;
4653 proposedSNR=1.0/(maxsnr*maxsnr*maxsnr) +
4654 (1.0/(minsnr*minsnr*minsnr)-1.0/(maxsnr*maxsnr*maxsnr))*XLALUniformDeviate(randParams);
4655 proposedSNR=1.0/cbrt(proposedSNR);
4656 return proposedSNR;
4657}
4658
4659static REAL8 draw_uniform_snr(REAL8 minsnr,REAL8 maxsnr)
4660{
4661 REAL8 proposedSNR=0.0;
4662 proposedSNR=minsnr+ (maxsnr-minsnr)*XLALUniformDeviate(randParams);
4663 return proposedSNR;
4664}
4665
4666static REAL8 draw_log10_snr(REAL8 minsnr,REAL8 maxsnr)
4667{
4668 REAL8 proposedlogSNR=0.0;
4669 REAL8 logminsnr=log10(minsnr);
4670 REAL8 logmaxsnr=log10(maxsnr);
4671 proposedlogSNR=logminsnr+ (logmaxsnr-logminsnr)*XLALUniformDeviate(randParams);
4672 return pow(10.0,proposedlogSNR);
4673}
program
const char * program
lalAppsVCSInfoList
const LALVCSInfoList lalAppsVCSInfoList
NULL-terminated list of VCS and build information for LALApps and its dependencies
lalAppsVCSIdentInfo
const LALVCSInfo lalAppsVCSIdentInfo
Identable VCS and build information for LALApps.
lal_errhandler
lal_errhandler_t lal_errhandler
LAL_ERR_EXIT
int LAL_ERR_EXIT(LALStatus *stat, const char *func, const char *file, const int line, volatile const char *id)
LAL_CALL
#define LAL_CALL(function, statusptr)
LALExtTriggerTableFromLIGOLw
int LALExtTriggerTableFromLIGOLw(ExtTriggerTable **eventHead, CHAR *fileName, INT4 startEvent, INT4 stopEvent)
Definition: LALExtTriggerTableFromLIGOLw.c:57
j
int j
k
int k
LALCheckMemoryLeaks
void LALCheckMemoryLeaks(void)
LALCalloc
#define LALCalloc(m, n)
LALFree
#define LALFree(p)
LALgetopt_long_only
int LALgetopt_long_only(int argc, char *const *argv, const char *options, const struct LALoption *long_options, int *opt_index)
LALoptarg
char * LALoptarg
no_argument
#define no_argument
required_argument
#define required_argument
XLALCloseLIGOLwXMLFile
int XLALCloseLIGOLwXMLFile(LIGOLwXMLStream *xml)
XLALOpenLIGOLwXMLFile
LIGOLwXMLStream * XLALOpenLIGOLwXMLFile(const char *path)
XLALWriteLIGOLwXMLProcessTable
int XLALWriteLIGOLwXMLProcessTable(LIGOLwXMLStream *, const ProcessTable *)
XLALWriteLIGOLwXMLProcessParamsTable
int XLALWriteLIGOLwXMLProcessParamsTable(LIGOLwXMLStream *, const ProcessParamsTable *)
XLALWriteLIGOLwXMLSimInspiralTable
int XLALWriteLIGOLwXMLSimInspiralTable(LIGOLwXMLStream *, const SimInspiralTable *)
XLALSimInspiralTableFromLIGOLw
SimInspiralTable * XLALSimInspiralTableFromLIGOLw(const char *fileName)
LALPopulateSimInspiralSiteInfo
void LALPopulateSimInspiralSiteInfo(LALStatus *status, SimInspiralTable *output)
LIGOMETA_STRING_MAX
#define LIGOMETA_STRING_MAX
LIGOMETA_INSPIRALTAPER_MAX
#define LIGOMETA_INSPIRALTAPER_MAX
LIGOMETA_COMMENT_MAX
#define LIGOMETA_COMMENT_MAX
LIGOMETA_TYPE_MAX
#define LIGOMETA_TYPE_MAX
LIGOMETA_WAVEFORM_MAX
#define LIGOMETA_WAVEFORM_MAX
LIGOMETA_PROGRAM_MAX
#define LIGOMETA_PROGRAM_MAX
LIGOMETA_SOURCE_MAX
#define LIGOMETA_SOURCE_MAX
LIGOMETA_VALUE_MAX
#define LIGOMETA_VALUE_MAX
LIGOMETA_PARAM_MAX
#define LIGOMETA_PARAM_MAX
XLALPopulateProcessTable
int XLALPopulateProcessTable(ProcessTable *ptable, const char *program_name, const char *cvs_revision, const char *cvs_source, const char *cvs_date, long process_id)
XLALSimInspiralAssignIDs
long XLALSimInspiralAssignIDs(SimInspiralTable *head, long process_id, long simulation_id)
dummy
INT4 dummy
fprintf
#define fprintf
e
double e
fLower
REAL4 fLower
Definition: blindinj.c:127
pp
const double pp
u
const double u
data
sigmaKerr data[0]
XLALDestroyREAL4VectorSequence
void XLALDestroyREAL4VectorSequence(REAL4VectorSequence *vecseq)
XLALGPSTimeNow
LIGOTimeGPS * XLALGPSTimeNow(LIGOTimeGPS *gpstime)
XLALCreateREAL8FrequencySeries
REAL8FrequencySeries * XLALCreateREAL8FrequencySeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaF, const LALUnit *sampleUnits, size_t length)
XLALDestroyREAL8FrequencySeries
void XLALDestroyREAL8FrequencySeries(REAL8FrequencySeries *series)
XLALRandomInspiralMilkywayLocation
void XLALRandomInspiralMilkywayLocation(REAL8 *rightAscension, REAL8 *declination, REAL8 *distance, RandomParams *randParams)
SpinDistribution
SpinDistribution
XLALFixedInspiralMasses
SimInspiralTable * XLALFixedInspiralMasses(SimInspiralTable *inj, REAL4 mass1Fix, REAL4 mass2Fix)
XLALRandomInspiralSkyLocation
SimInspiralTable * XLALRandomInspiralSkyLocation(SimInspiralTable *inj, RandomParams *randParams)
XLALRandomInspiralSpins
SimInspiralTable * XLALRandomInspiralSpins(SimInspiralTable *inj, RandomParams *randParams, REAL4 spin1Min, REAL4 spin1Max, REAL4 spin2Min, REAL4 spin2Max, REAL4 kappa1Min, REAL4 kappa1Max, REAL4 abskappa1Min, REAL4 abskappa1Max, AlignmentType alignInj, SpinDistribution distribution, REAL4 spin1Mean, REAL4 spin1Std, REAL4 spin2Mean, REAL4 spin2Std)
LoudnessDistribution
LoudnessDistribution
MassDistribution
MassDistribution
XLALRandomInspiralDistance
SimInspiralTable * XLALRandomInspiralDistance(SimInspiralTable *inj, RandomParams *randParams, LoudnessDistribution dDist, REAL4 distMin, REAL4 distMax)
XLALGaussianInspiralMasses
SimInspiralTable * XLALGaussianInspiralMasses(SimInspiralTable *inj, RandomParams *randParams, REAL4 mass1Min, REAL4 mass1Max, REAL4 mass1Mean, REAL4 mass1Std, REAL4 mass2Min, REAL4 mass2Max, REAL4 mass2Mean, REAL4 mass2Std)
lalinspiral_time_distribution
lalinspiral_time_distribution
XLALRandomInspiralTime
SimInspiralTable * XLALRandomInspiralTime(SimInspiralTable *inj, RandomParams *randParams, LIGOTimeGPS startTime, REAL4 timeWindow)
XLALRandomInspiralTotalMassFraction
SimInspiralTable * XLALRandomInspiralTotalMassFraction(SimInspiralTable *inj, RandomParams *randParams, MassDistribution mDist, REAL4 minTotalMass, REAL4 maxTotalMass, REAL4 minMassRatio, REAL4 maxMassRatio)
InclDistribution
InclDistribution
AlignmentType
AlignmentType
XLALm1m2SquareGridInspiralMasses
SimInspiralTable * XLALm1m2SquareGridInspiralMasses(SimInspiralTable *inj, REAL4 mass1Min, REAL4 mass2Min, REAL4 minTotalMass, REAL4 maxTotalMass, REAL4 mass1Delta, REAL4 mass2Delta, INT4 mass1Pnt, INT4 mass2Pnt, INT4 injNum, INT4 *count)
XLALRandomNRInjectTotalMass
SimInspiralTable * XLALRandomNRInjectTotalMass(SimInspiralTable *inj, RandomParams *randParams, REAL4 minTotalMass, REAL4 maxTotalMass, SimInspiralTable *nrInjParams)
SkyLocationDistribution
SkyLocationDistribution
XLALRandomInspiralMasses
SimInspiralTable * XLALRandomInspiralMasses(SimInspiralTable *inj, RandomParams *randParams, MassDistribution mDistr, REAL4 mass1Min, REAL4 mass1Max, REAL4 mass2Min, REAL4 mass2Max, REAL4 minTotalMass, REAL4 maxTotalMass)
XLALRandomInspiralOrientation
SimInspiralTable * XLALRandomInspiralOrientation(SimInspiralTable *inj, RandomParams *randParams, InclDistribution iDist, REAL4 inclinationPeak)
XLALRandomInspiralTotalMassRatio
SimInspiralTable * XLALRandomInspiralTotalMassRatio(SimInspiralTable *inj, RandomParams *randParams, MassDistribution mDist, REAL4 minTotalMass, REAL4 maxTotalMass, REAL4 minMassRatio, REAL4 maxMassRatio)
gaussianSpinDist
gaussianSpinDist
unknownSpinDist
unknownSpinDist
uniformSpinDist
uniformSpinDist
starFormationRate
starFormationRate
distFromSourceFile
distFromSourceFile
unknownLoudnessDist
unknownLoudnessDist
uniformSnr
uniformSnr
uniformLogDistance
uniformLogDistance
uniformDistance
uniformDistance
uniformVolume
uniformVolume
uniformLogSnr
uniformLogSnr
uniformVolumeSnr
uniformVolumeSnr
uniformDistanceSquared
uniformDistanceSquared
logComponentMass
logComponentMass
m1m2SquareGrid
m1m2SquareGrid
unknownMassDist
unknownMassDist
uniformTotalMass
uniformTotalMass
uniformTotalMassRatio
uniformTotalMassRatio
uniformTotalMassFraction
uniformTotalMassFraction
massFromNRFile
massFromNRFile
massFromSourceFile
massFromSourceFile
fixMasses
fixMasses
uniformComponentMass
uniformComponentMass
gaussianMassDist
gaussianMassDist
logMassUniformTotalMassRatio
logMassUniformTotalMassRatio
LALINSPIRAL_UNIFORM_TIME_DIST
LALINSPIRAL_UNIFORM_TIME_DIST
LALINSPIRAL_FILE_TIME_DIST
LALINSPIRAL_FILE_TIME_DIST
LALINSPIRAL_UNKNOWN_TIME_DIST
LALINSPIRAL_UNKNOWN_TIME_DIST
LALINSPIRAL_EXPONENTIAL_TIME_DIST
LALINSPIRAL_EXPONENTIAL_TIME_DIST
LALINSPIRAL_FIXED_TIME_DIST
LALINSPIRAL_FIXED_TIME_DIST
fixedInclDist
fixedInclDist
uniformInclDist
uniformInclDist
unknownInclDist
unknownInclDist
gaussianInclDist
gaussianInclDist
inxzPlane
inxzPlane
alongzAxis
alongzAxis
notAligned
notAligned
locationFromExttrigFile
locationFromExttrigFile
locationFromSourceFile
locationFromSourceFile
uniformSkyLocation
uniformSkyLocation
fixedSkyLocation
fixedSkyLocation
locationFromIPNFile
locationFromIPNFile
unknownLocationDist
unknownLocationDist
LAL_PI_180
#define LAL_PI_180
LAL_PI
#define LAL_PI
XLAL_INIT_MEM
#define XLAL_INIT_MEM(x)
REAL8
double REAL8
CHAR
char CHAR
UINT4
uint32_t UINT4
INT4
int32_t INT4
REAL4
float REAL4
XLALFree
void XLALFree(void *p)
LALSimInspiralApplyTaper
LALSimInspiralApplyTaper
LAL_SIM_INSPIRAL_TAPER_START
LAL_SIM_INSPIRAL_TAPER_START
LAL_SIM_INSPIRAL_TAPER_STARTEND
LAL_SIM_INSPIRAL_TAPER_STARTEND
LAL_SIM_INSPIRAL_TAPER_END
LAL_SIM_INSPIRAL_TAPER_END
LAL_SIM_INSPIRAL_TAPER_NONE
LAL_SIM_INSPIRAL_TAPER_NONE
XLALOutputVCSInfo
int XLALOutputVCSInfo(FILE *fp, const LALVCSInfoList vcs_list, const int verbose, const char *prefix)
q
static const INT4 q
LALCreateRandomParams
void LALCreateRandomParams(LALStatus *status, RandomParams **params, INT4 seed)
XLALUniformDeviate
REAL4 XLALUniformDeviate(RandomParams *params)
LALDestroyRandomParams
void LALDestroyRandomParams(LALStatus *status, RandomParams **params)
XLALDestroyREAL8TimeSeries
void XLALDestroyREAL8TimeSeries(REAL8TimeSeries *series)
lalHertzUnit
const LALUnit lalHertzUnit
XLALCreateREAL8Vector
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
XLALDestroyREAL8Vector
void XLALDestroyREAL8Vector(REAL8Vector *vector)
XLAL_ERROR
#define XLAL_ERROR(...)
XLAL_SUCCESS
XLAL_SUCCESS
XLAL_EIO
XLAL_EIO
XLALGreenwichMeanSiderealTime
REAL8 XLALGreenwichMeanSiderealTime(const LIGOTimeGPS *gpstime)
XLALGPSAdd
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
XLALGPSCmp
int XLALGPSCmp(const LIGOTimeGPS *t0, const LIGOTimeGPS *t1)
ifo
char * ifo
Definition: gwf2xml.c:57
count
long long count
gpsStartTime
LIGOTimeGPS gpsStartTime
Definition: inspfrinj.c:314
gpsEndTime
LIGOTimeGPS gpsEndTime
Definition: inspfrinj.c:316
userTag
CHAR * userTag
Definition: inspfrinj.c:343
injections
SimInspiralTable * injections
Definition: inspfrinj.c:339
sourceComplete
void sourceComplete(void)
Definition: inspinj.c:1428
fixed_inc
REAL4 fixed_inc
Definition: inspinj.c:515
spinDistr
SpinDistribution spinDistr
Definition: inspinj.c:472
dDistr
LoudnessDistribution dDistr
Definition: inspinj.c:468
main
int main(int argc, char *argv[])
Definition: inspinj.c:1875
minMass1
REAL4 minMass1
Definition: inspinj.c:501
lum
REAL8 lum
Definition: inspinj.c:565
ninjaSNR
INT4 ninjaSNR
Definition: inspinj.c:488
mwLuminosity
REAL4 mwLuminosity
Definition: inspinj.c:491
MW_name
char MW_name[LIGOMETA_SOURCE_MAX]
Definition: inspinj.c:569
inclStd
REAL4 inclStd
Definition: inspinj.c:514
maxZ
REAL4 maxZ
Definition: inspinj.c:495
ifos
char * ifos
Definition: inspinj.c:499
spinInjections
int spinInjections
Definition: inspinj.c:523
norm
REAL8 norm
Definition: inspinj.c:572
randParams
static RandomParams * randParams
Definition: inspinj.c:553
maxabsKappa1
REAL4 maxabsKappa1
Definition: inspinj.c:536
maxMass1
REAL4 maxMass1
Definition: inspinj.c:502
srcCompleteDist
REAL8 srcCompleteDist
Definition: inspinj.c:593
fudge
REAL8 fudge
Definition: inspinj.c:566
maxKappa1
REAL4 maxKappa1
Definition: inspinj.c:534
probability_redshift
REAL8 probability_redshift(REAL8 rshift)
Definition: inspinj.c:605
meanSpin1
REAL4 meanSpin1
Definition: inspinj.c:527
maxMtotal
REAL4 maxMtotal
Definition: inspinj.c:506
draw_log10_snr
static REAL8 draw_log10_snr(REAL8 snrmin, REAL8 snrmax)
Definition: inspinj.c:4666
redshift_mass
REAL8 redshift_mass(REAL8 mass, REAL8 z)
Definition: inspinj.c:653
old_source_data
struct @3 * old_source_data
pntMass2
INT4 pntMass2
Definition: inspinj.c:540
numExtTriggers
INT4 numExtTriggers
Definition: inspinj.c:554
PROGRAM_NAME
#define PROGRAM_NAME
Definition: inspinj.c:396
ninjaMass
INT4 ninjaMass
Definition: inspinj.c:486
massFileName
char * massFileName
Definition: inspinj.c:477
lDistr
SkyLocationDistribution lDistr
Definition: inspinj.c:469
massStdev2
REAL4 massStdev2
Definition: inspinj.c:510
drawLocationFromExttrig
void drawLocationFromExttrig(SimInspiralTable *table)
Definition: inspinj.c:1846
adjust_snr_with_psds_real8
void adjust_snr_with_psds_real8(SimInspiralTable *inj, REAL8 target_snr, int num_ifos, const char **ifo_list, REAL8FrequencySeries **psds, REAL8 *start_freqs)
Definition: inspinj.c:903
num_mass
int num_mass
Definition: inspinj.c:574
name
char name[LIGOMETA_SOURCE_MAX]
Definition: inspinj.c:561
read_nr_data
void read_nr_data(char *filename)
Definition: inspinj.c:1225
read_time_data
void read_time_data(char *filename)
Definition: inspinj.c:1168
read_mass_data
void read_mass_data(char *filename)
Definition: inspinj.c:1125
massStdev1
REAL4 massStdev1
Definition: inspinj.c:509
num_nr
int num_nr
Definition: inspinj.c:595
read_IPN_grid_from_file
void read_IPN_grid_from_file(char *fname)
Definition: inspinj.c:1368
localRate
REAL4 localRate
Definition: inspinj.c:496
maxMassRatio
REAL4 maxMassRatio
Definition: inspinj.c:512
Spin2Std
REAL4 Spin2Std
Definition: inspinj.c:532
minSNR
REAL4 minSNR
Definition: inspinj.c:497
redshift
REAL8 redshift
Definition: inspinj.c:546
drawRedshift
REAL8 drawRedshift(REAL8 zmin, REAL8 zmax, REAL8 pzmax)
Definition: inspinj.c:639
bandPassInj
INT4 bandPassInj
Definition: inspinj.c:543
minabsKappa1
REAL4 minabsKappa1
Definition: inspinj.c:535
snr_in_psd_real8
REAL8 snr_in_psd_real8(const char *ifo, REAL8FrequencySeries *psd, REAL8 start_freq, SimInspiralTable *inj)
Definition: inspinj.c:867
mDistr
MassDistribution mDistr
Definition: inspinj.c:470
haveLoudness
INT4 haveLoudness
Definition: inspinj.c:489
network_snr
REAL8 network_snr(const char *ifos, SimInspiralTable *inj)
Definition: inspinj.c:684
minMass2
REAL4 minMass2
Definition: inspinj.c:503
fracVec
REAL8 * fracVec
Definition: inspinj.c:570
single_IFO_SNR_threshold
REAL8 single_IFO_SNR_threshold
Definition: inspinj.c:548
drawFromSource
void drawFromSource(REAL8 *rightAscension, REAL8 *declination, REAL8 *distance, CHAR name[LIGOMETA_SOURCE_MAX])
Definition: inspinj.c:1779
outCompress
INT4 outCompress
Definition: inspinj.c:485
network_snr_with_psds_real8
REAL8 network_snr_with_psds_real8(int num_ifos, const char **ifo_list, REAL8FrequencySeries **psds, REAL8 *start_freqs, SimInspiralTable *inj)
Definition: inspinj.c:884
minSpin1
REAL4 minSpin1
Definition: inspinj.c:525
drawMassFromSource
void drawMassFromSource(SimInspiralTable *table)
Definition: inspinj.c:1662
minSpin2
REAL4 minSpin2
Definition: inspinj.c:529
psi
REAL4 psi
Definition: inspinj.c:519
iDistr
InclDistribution iDistr
Definition: inspinj.c:471
longitude
REAL4 longitude
Definition: inspinj.c:520
network_snr_real8
REAL8 network_snr_real8(const char *ifos, SimInspiralTable *inj)
Definition: inspinj.c:779
read_source_data
void read_source_data(char *filename)
Definition: inspinj.c:1274
minKappa1
REAL4 minKappa1
Definition: inspinj.c:533
drawMassSpinFromNR
void drawMassSpinFromNR(SimInspiralTable *table)
Definition: inspinj.c:1687
source_data
struct @3 * source_data
n_times
int n_times
Definition: inspinj.c:580
deltaMass2
REAL4 deltaMass2
Definition: inspinj.c:542
draw_uniform_snr
static REAL8 draw_uniform_snr(REAL8 snrmin, REAL8 snrmax)
Definition: inspinj.c:4659
Spin1Std
REAL4 Spin1Std
Definition: inspinj.c:528
alignInj
AlignmentType alignInj
Definition: inspinj.c:545
sourceFileName
char * sourceFileName
Definition: inspinj.c:479
AXIS_MAX
#define AXIS_MAX
Definition: inspinj.c:414
status
static LALStatus status
Definition: inspinj.c:552
simTable
SimInspiralTable * simTable
Definition: inspinj.c:474
maxD
REAL4 maxD
Definition: inspinj.c:493
srcComplete
int srcComplete
Definition: inspinj.c:591
minMassRatio
REAL4 minMassRatio
Definition: inspinj.c:511
pntMass1
INT4 pntMass1
Definition: inspinj.c:539
fixedCoaPhase
REAL4 fixedCoaPhase
Definition: inspinj.c:518
makeCatalog
int makeCatalog
Definition: inspinj.c:592
scale_lalsim_distance
static void scale_lalsim_distance(SimInspiralTable *inj, char **IFOnames, REAL8FrequencySeries **psds, REAL8 *start_freqs, LoudnessDistribution dDistr)
Definition: inspinj.c:4534
vrbflg
int vrbflg
defined in lal/lib/std/LALError.c
meanMass1
REAL4 meanMass1
Definition: inspinj.c:507
luminosity_distance
REAL8 luminosity_distance(REAL8 rshift)
Definition: inspinj.c:616
inj_times
LIGOTimeGPS * inj_times
Definition: inspinj.c:581
adjust_snr
void adjust_snr(SimInspiralTable *inj, REAL8 target_snr, const char *ifos)
Definition: inspinj.c:707
adjust_snr_real8
void adjust_snr_real8(SimInspiralTable *inj, REAL8 target_snr, const char *ifos)
Definition: inspinj.c:804
tDistr
lalinspiral_time_distribution tDistr
Definition: inspinj.c:467
minD
REAL4 minD
Definition: inspinj.c:492
exttrigHead
ExtTriggerTable * exttrigHead
Definition: inspinj.c:555
nrFileName
char * nrFileName
Definition: inspinj.c:478
fixedMass2
REAL4 fixedMass2
Definition: inspinj.c:538
maxSNR
REAL4 maxSNR
Definition: inspinj.c:498
nrSimArray
SimInspiralTable ** nrSimArray
Definition: inspinj.c:597
maxSpin1
REAL4 maxSpin1
Definition: inspinj.c:526
outputFileName
char * outputFileName
Definition: inspinj.c:480
draw_volume_snr
static REAL8 draw_volume_snr(REAL8 snrmin, REAL8 snrmax)
Definition: inspinj.c:4650
ifonames
char ** ifonames
Definition: inspinj.c:549
meanMass2
REAL4 meanMass2
Definition: inspinj.c:508
i
int i
Definition: inspinj.c:596
minZ
REAL4 minZ
Definition: inspinj.c:494
mean_time_step_sfr
REAL8 mean_time_step_sfr(REAL8 zmax, REAL8 rate_local)
Definition: inspinj.c:627
taperInj
LALSimInspiralApplyTaper taperInj
Definition: inspinj.c:544
injtimesFileName
char * injtimesFileName
Definition: inspinj.c:481
print_usage
static void print_usage(char *program)
Definition: inspinj.c:949
next_process_param
ProcessParamsTable * next_process_param(const char *name, const char *type, const char *fmt,...)
Definition: inspinj.c:924
maxSpin2
REAL4 maxSpin2
Definition: inspinj.c:530
spinAligned
int spinAligned
Definition: inspinj.c:524
deltaMass1
REAL4 deltaMass1
Definition: inspinj.c:541
exttrigFileName
char * exttrigFileName
Definition: inspinj.c:482
meanSpin2
REAL4 meanSpin2
Definition: inspinj.c:531
minMtotal
REAL4 minMtotal
Definition: inspinj.c:505
numifos
int numifos
Definition: inspinj.c:550
IPNSkyPositionsFile
char * IPNSkyPositionsFile
Definition: inspinj.c:483
snr_in_ifo_real8
REAL8 snr_in_ifo_real8(const char *ifo, SimInspiralTable *inj)
Definition: inspinj.c:765
coaPhaseFixed
int coaPhaseFixed
Definition: inspinj.c:517
maxMass2
REAL4 maxMass2
Definition: inspinj.c:504
mass1
REAL8 mass1
Definition: inspinj.c:576
dist
REAL8 dist
Definition: inspinj.c:564
max_inc
REAL4 max_inc
Definition: inspinj.c:516
galaxynum
int galaxynum
Definition: inspinj.c:559
temparray
struct @3 * temparray
mass2
REAL8 mass2
Definition: inspinj.c:577
num_source
int num_source
Definition: inspinj.c:557
mass_data
struct @4 * mass_data
ratioVec
REAL8 * ratioVec
Definition: inspinj.c:571
skyPoints
struct @3 * skyPoints
simRingTable
SimRingdownTable * simRingTable
Definition: inspinj.c:475
epsAngle
REAL4 epsAngle
Definition: inspinj.c:522
latitude
REAL4 latitude
Definition: inspinj.c:521
drawMassSpinFromNRNinja2
void drawMassSpinFromNRNinja2(SimInspiralTable *table)
Definition: inspinj.c:1699
numSkyPoints
int numSkyPoints
Definition: inspinj.c:558
fixedMass1
REAL4 fixedMass1
Definition: inspinj.c:537
real8Ninja2
INT4 real8Ninja2
Definition: inspinj.c:487
drawFromIPNsim
void drawFromIPNsim(REAL8 *rightAscension, REAL8 *declination)
Definition: inspinj.c:1819
snr_in_ifo
REAL8 snr_in_ifo(const char *ifo, SimInspiralTable *inj)
Definition: inspinj.c:667
dec
REAL8 dec
Definition: inspinj.c:563
ra
REAL8 ra
Definition: inspinj.c:562
calculate_lalsim_snr
REAL8 calculate_lalsim_snr(SimInspiralTable *inj, char *IFOname, REAL8FrequencySeries *psd, REAL8 start_freq)
Definition: inspiralutils.c:1010
calculate_snr_from_strain_and_psd_real8
REAL8 calculate_snr_from_strain_and_psd_real8(REAL8TimeSeries *strain, REAL8FrequencySeries *psd, REAL8 startFreq, const CHAR ifo[3])
Definition: inspiralutils.c:722
XLALPsdFromFile
int XLALPsdFromFile(REAL8FrequencySeries **psd, const CHAR *filename)
Definition: inspiralutils.c:860
get_FakePsdFromString
void get_FakePsdFromString(REAL8FrequencySeries *PsdFreqSeries, char *FakePsdName, REAL8 StartFreq)
Definition: inspiralutils.c:983
calculate_ligo_snr_from_strain_real8
REAL8 calculate_ligo_snr_from_strain_real8(REAL8TimeSeries *strain, const CHAR ifo[3])
Definition: inspiralutils.c:667
calculate_ligo_snr_from_strain
REAL8 calculate_ligo_snr_from_strain(REAL4TimeVectorSeries *strain, SimInspiralTable *thisInj, const CHAR ifo[3])
Definition: inspiralutils.c:798
AddNumRelStrainModes
void AddNumRelStrainModes(LALStatus *status, REAL4TimeVectorSeries **outStrain, SimInspiralTable *thisinj)
Definition: inspiralutils.c:189
start_freq_from_frame_url
REAL8 start_freq_from_frame_url(CHAR *url)
Definition: inspiralutils.c:620
XLALNRInjectionStrain
REAL8TimeSeries * XLALNRInjectionStrain(const char *ifo, SimInspiralTable *inj)
type
type
head
head
eta
eta
m2
m2
m1
m1
waveform
waveform
fmt
fmt
lalapps_power_pipe.filename
filename
Definition: lalapps_power_pipe.py:266
lalapps_string_contour_plotter.p
p
Definition: lalapps_string_contour_plotter.py:102
lalapps_string_contour_plotter.n
n
Definition: lalapps_string_contour_plotter.py:102
lalapps_string_cs_gamma_largeloops.z
z
Definition: lalapps_string_cs_gamma_largeloops.py:99
c
c
seglen
seglen
psd
psd
int
int
line
string line
strain
strain
rand_seed
rand_seed
log
log
N
int N
atol
float atol
xmlfp
LIGOLwXMLStream * xmlfp
Definition: spininj.c:210
fname
CHAR fname[256]
Definition: spininj.c:211
FakeGalaxy::ra
REAL8 ra
Definition: inspinj.c:585
FakeGalaxy::name
char name[LIGOMETA_SOURCE_MAX]
Definition: inspinj.c:584
FakeGalaxy::dec
REAL8 dec
Definition: inspinj.c:586
FakeGalaxy::next
struct FakeGalaxy * next
Definition: inspinj.c:590
FakeGalaxy::lum
REAL8 lum
Definition: inspinj.c:587
FakeGalaxy::fudge
REAL8 fudge
Definition: inspinj.c:589
FakeGalaxy::dist
REAL8 dist
Definition: inspinj.c:588
LALVCSInfo::vcsDate
const char *const vcsDate
LALVCSInfo::vcsStatus
const char *const vcsStatus
LALVCSInfo::vcsId
const char *const vcsId
LALoption::flag
int * flag
LIGOTimeGPS::gpsNanoSeconds
INT4 gpsNanoSeconds
ProcessParamsTable::type
CHAR type[LIGOMETA_TYPE_MAX]
ProcessParamsTable::param
CHAR param[LIGOMETA_PARAM_MAX]
ProcessParamsTable::value
CHAR value[LIGOMETA_VALUE_MAX]
ProcessParamsTable::next
struct tagProcessParamsTable * next
ProcessParamsTable::program
CHAR program[LIGOMETA_PROGRAM_MAX]
ProcessTable::start_time
LIGOTimeGPS start_time
ProcessTable::end_time
LIGOTimeGPS end_time
ProcessTable::comment
CHAR comment[LIGOMETA_COMMENT_MAX]
REAL4TimeVectorSeries::data
REAL4VectorSequence * data
REAL8Vector::data
REAL8 * data
SimInspiralTable::geocent_end_time
LIGOTimeGPS geocent_end_time
SimInspiralTable::source
CHAR source[LIGOMETA_SOURCE_MAX]
SimInspiralTable::next
struct tagSimInspiralTable * next
SimInspiralTable::taper
CHAR taper[LIGOMETA_INSPIRALTAPER_MAX]
SimInspiralTable::waveform
CHAR waveform[LIGOMETA_WAVEFORM_MAX]
SimInspiralTable::numrel_data
CHAR numrel_data[LIGOMETA_STRING_MAX]
SimRingdownTable::next
struct tagSimRingdownTable * next
fp
FILE * fp
alpha
REAL4 alpha
Definition: tmpltbank.c:432
srate
double srate