00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include <stdio.h>
00032 #include <stdlib.h>
00033 #include <config.h>
00034
00035 #include <math.h>
00036 #include <ctype.h>
00037 #include <assert.h>
00038 #include <string.h>
00039 #include <getopt.h>
00040 #include <time.h>
00041 #include <lalapps.h>
00042 #include <processtable.h>
00043 #include <lal/LALStdio.h>
00044 #include <lal/LALStdlib.h>
00045 #include <lal/LALConstants.h>
00046 #include <lal/LIGOMetadataTables.h>
00047 #include <lal/LIGOMetadataUtils.h>
00048 #include <lal/LIGOLwXML.h>
00049 #include <lal/Units.h>
00050 #include <lal/Date.h>
00051 #include <lal/Inject.h>
00052 #include <lal/InspiralInjectionParams.h>
00053 #include <lal/GenerateInspiral.h>
00054 #include <lal/GenerateInspRing.h>
00055 #include <lal/FindChirp.h>
00056 #include <lal/GenerateRing.h>
00057 #include <lal/Ring.h>
00058 #include <lal/LALNoiseModels.h>
00059 #include <lal/RealFFT.h>
00060 #include <lal/FrequencySeries.h>
00061 #include <lal/TimeSeries.h>
00062 #include <lal/TimeFreqFFT.h>
00063 #include <lal/VectorOps.h>
00064
00065 RCSID( "$Id: blindinj.c,v 1.17 2008/07/03 23:50:25 dfazi Exp $" );
00066 #define CVS_ID_STRING "$Id: blindinj.c,v 1.17 2008/07/03 23:50:25 dfazi Exp $"
00067 #define CVS_NAME_STRING "$Name: $"
00068 #define CVS_REVISION "$Revision: 1.17 $"
00069 #define CVS_SOURCE "$Source: /usr/local/cvs/lscsoft/lalapps/src/inspiral/blindinj.c,v $"
00070 #define CVS_DATE "$Date: 2008/07/03 23:50:25 $"
00071 #define PROGRAM_NAME "blindinj"
00072
00073 #define USAGE \
00074 "lalapps_blindinj [options]\n"\
00075 "\nDefaults are shown in brackets\n\n" \
00076 " --help display this message\n"\
00077 " --version print version information and exit\n"\
00078 " --verbose be verbose\n"\
00079 " --gps-start-time TIME start time of injection\n"\
00080 " --injection-type TYPE type of injection, must be one of \n"\
00081 " (strain, etmx, etmy)\n"\
00082 " --seed SEED seed random number generator with SEED (1)\n"\
00083 " --debug-level LEVEL set the LAL debug level to LEVEL\n"\
00084 "\n"
00085
00086
00087
00088 typedef enum
00089 {
00090 noResponse,
00091 unityResponse,
00092 LIGOdesign,
00093 actuationX,
00094 actuationY
00095 }
00096 ResponseFunction;
00097
00098 typedef struct actuationparameters {
00099 REAL4 ETMXcal;
00100 REAL4 ETMYcal;
00101 REAL4 pendFX;
00102 REAL4 pendFY;
00103 REAL4 pendQX;
00104 REAL4 pendQY;
00105 REAL4 length;
00106 } ActuationParameters;
00107
00108 static ProcessParamsTable *next_process_param(
00109 const char *name,
00110 const char *type,
00111 const char *fmt, ... );
00112
00113 static void destroyCoherentGW( CoherentGW *waveform );
00114
00115 static REAL4TimeSeries *injectWaveform(
00116 LALStatus *status,
00117 SimInspiralTable *inspInj,
00118 SimRingdownTable *ringdownevents,
00119 ResponseFunction responseType,
00120 InterferometerNumber ifoNumber,
00121 LIGOTimeGPS epoch);
00122
00123 extern int vrbflg;
00124
00125
00126 ActuationParameters actuationParams[LAL_NUM_IFO];
00127 REAL8 dynRange = 1.0/3.0e-23;
00128 INT4 duration = 100.0;
00129 INT4 sampleRate = 16384;
00130
00131
00132 REAL4 fLower = 30;
00133 REAL4 mergerLength = 10;
00134 REAL8 longestSignal = 95.0;
00135 REAL8 timeWindow = 4;
00136
00137 UINT4 numInjections = 0;
00138 INT4 randSeed = 1;
00139 REAL4 minNSMass = 1.0;
00140 REAL4 maxNSMass = 2.0;
00141 REAL4 minBHMass = 2.0;
00142 REAL4 maxBHMass = 30.0;
00143 REAL4 minTotalMass = 0.0;
00144 REAL4 maxTotalMass = 35.0;
00145
00146 REAL4 minNSSpin = 0.0;
00147 REAL4 maxNSSpin = 0.2;
00148 REAL4 minBHSpin = 0.0;
00149 REAL4 maxBHSpin = 1.0;
00150
00151 REAL4 BNSfrac = 0.35;
00152 REAL4 BBHfrac = 0.35;
00153
00154
00155 REAL4 bnsSnrMean = 9.0;
00156 REAL4 bnsSnrStd = 0.5;
00157 REAL4 snrMean = 12.0;
00158 REAL4 snrStd = 1.0;
00159
00160 REAL4Vector *normalDev;
00161
00162
00163
00164 ProcessParamsTable *next_process_param(
00165 const char *name,
00166 const char *type,
00167 const char *fmt, ... )
00168 {
00169 ProcessParamsTable *pp;
00170 va_list ap;
00171 pp = calloc( 1, sizeof( *pp ) );
00172
00173 if ( ! pp )
00174 {
00175 perror( "next_process_param" );
00176 exit( 1 );
00177 }
00178 strncpy( pp->program, PROGRAM_NAME, LIGOMETA_PROGRAM_MAX );
00179 if ( ! strcmp( name, "userTag" ) || ! strcmp( name, "user-tag" ) )
00180 LALSnprintf( pp->param, LIGOMETA_PARAM_MAX, "-userTag" );
00181 else
00182 LALSnprintf( pp->param, LIGOMETA_PARAM_MAX, "--%s", name );
00183 strncpy( pp->type, type, LIGOMETA_TYPE_MAX );
00184 va_start( ap, fmt );
00185 LALVsnprintf( pp->value, LIGOMETA_VALUE_MAX, fmt, ap );
00186 va_end( ap );
00187
00188 return pp;
00189 }
00190
00191 static void destroyCoherentGW( CoherentGW *waveform )
00192 {
00193 if ( waveform->h )
00194 {
00195 XLALDestroyREAL4VectorSequence( waveform->h->data );
00196 LALFree( waveform->a );
00197 }
00198 if ( waveform->a )
00199 {
00200 XLALDestroyREAL4VectorSequence( waveform->a->data );
00201 LALFree( waveform->a );
00202 }
00203 if ( waveform->phi )
00204 {
00205 XLALDestroyREAL8Vector( waveform->phi->data );
00206 LALFree( waveform->phi );
00207 }
00208 if ( waveform->f )
00209 {
00210 XLALDestroyREAL4Vector( waveform->f->data );
00211 LALFree( waveform->f );
00212 }
00213 if ( waveform->shift )
00214 {
00215 XLALDestroyREAL4Vector( waveform->shift->data );
00216 LALFree( waveform->shift );
00217 }
00218
00219 return;
00220 }
00221
00222 static REAL4TimeSeries *injectWaveform(
00223 LALStatus *status,
00224 SimInspiralTable *inspInj,
00225 SimRingdownTable *ringdownevents,
00226 ResponseFunction responseType,
00227 InterferometerNumber ifoNumber,
00228 LIGOTimeGPS epoch)
00229 {
00230 REAL4TimeSeries *chan;
00231 COMPLEX8FrequencySeries *resp;
00232 COMPLEX8Vector *unity;
00233 CHAR name[LALNameLength];
00234 CHAR ifo[LIGOMETA_IFO_MAX];
00235 PPNParamStruc ppnParams;
00236 SimRingdownTable *thisRingdownEvent = NULL;
00237
00238 INT8 waveformStartTime;
00239 DetectorResponse detector;
00240 CoherentGW waveform, *wfm;
00241 ActuationParameters actData = actuationParams[ifoNumber];
00242 UINT4 i,k;
00243 int injectSignalType = imr_inject;
00244 const LALUnit strainPerCount = {0,{0,0,0,0,0,1,-1},{0,0,0,0,0,0,0}};
00245 FILE *fp = NULL;
00246 char fileName[FILENAME_MAX];
00247
00248
00249 XLALReturnIFO( ifo, ifoNumber );
00250 LALSnprintf( name, LALNameLength * sizeof(CHAR), "%s:INJECT", ifo );
00251 chan = XLALCreateREAL4TimeSeries( name, &epoch, 0, 1./sampleRate,
00252 &lalADCCountUnit, sampleRate * duration );
00253 if ( ! chan )
00254 {
00255 exit( 1 );
00256 }
00257
00258 memset( chan->data->data, 0, chan->data->length * sizeof(REAL4) );
00259
00260 thisRingdownEvent = ringdownevents;
00261
00262
00263
00264
00265
00266
00267
00268
00269 memset( &ppnParams, 0, sizeof(PPNParamStruc) );
00270 ppnParams.deltaT = chan->deltaT;
00271 ppnParams.lengthIn = 0;
00272 ppnParams.ppn = NULL;
00273
00274 memset( &waveform, 0, sizeof(CoherentGW) );
00275
00276 LAL_CALL( LALGenerateInspiral(status, &waveform, inspInj, &ppnParams),
00277 status);
00278
00279
00280 wfm = XLALGenerateInspRing( &waveform, inspInj, thisRingdownEvent,
00281 injectSignalType );
00282
00283 if ( !wfm )
00284 {
00285 fprintf( stderr, "Failed to generate the waveform \n" );
00286 if (xlalErrno == XLAL_EFAILED)
00287 {
00288 fprintf( stderr, "Too much merger\n");
00289 XLALDestroyREAL4TimeSeries( chan );
00290 xlalErrno = XLAL_SUCCESS;
00291 return ( NULL );
00292 }
00293 else exit ( 1 );
00294 }
00295
00296 waveform = *wfm;
00297
00298
00299 if ( ifoNumber == LAL_IFO_H1 && vrbflg )
00300 {
00301 fprintf( stdout,
00302 "Writing out A+, Ax, f, phi, shift for waveform "
00303 "to file named INSPIRAL_WAVEFORM.dat\n");
00304 LALSnprintf( fileName, FILENAME_MAX, "INSPIRAL_WAVEFORM.dat");
00305 fp = fopen(fileName, "w");
00306 for ( i = 0; i < waveform.phi->data->length; i++ )
00307 {
00308 if ( waveform.shift ) fprintf( fp, "%e\t %e\t %f\t %f\t %f\n",
00309 waveform.a->data->data[2*i],
00310 waveform.a->data->data[2*i+1], waveform.f->data->data[i],
00311 waveform.phi->data->data[i] , waveform.shift->data->data[i] );
00312
00313 else fprintf( fp, "%e\t %e\t %f\t %f\n",
00314 waveform.a->data->data[2*i],
00315 waveform.a->data->data[2*i+1], waveform.f->data->data[i],
00316 waveform.phi->data->data[i] );
00317 }
00318 fclose( fp );
00319 }
00320
00321
00322
00323
00324
00325
00326
00327 resp = XLALCreateCOMPLEX8FrequencySeries( chan->name, &(chan->epoch), 0,
00328 1.0 / ( duration ), &strainPerCount, ( sampleRate * duration / 2 + 1 ) );
00329
00330 switch ( responseType )
00331 {
00332 case noResponse:
00333 fprintf( stderr, "Must specify the response function\n" );
00334 exit( 1 );
00335 break;
00336
00337 case unityResponse:
00338
00339 for ( k = 0; k < resp->data->length; ++k )
00340 {
00341 resp->data->data[k].re = (REAL4) (1.0 /dynRange);
00342 resp->data->data[k].im = 0.0;
00343 }
00344 break;
00345
00346 case LIGOdesign:
00347
00348 for ( k = 0; k < resp->data->length; ++k )
00349 {
00350 REAL8 sim_psd_freq = (REAL8) k * resp->deltaF;
00351 REAL8 sim_psd_value;
00352 LALLIGOIPsd( NULL, &sim_psd_value, sim_psd_freq );
00353 resp->data->data[k].re = (REAL4) pow( sim_psd_value, 0.5 ) / dynRange;
00354 resp->data->data[k].im = 0.0;
00355 }
00356 break;
00357
00358 case actuationX:
00359
00360 resp = generateActuation( resp, actData.ETMXcal / actData.length,
00361 actData.pendFX, actData.pendQX);
00362 break;
00363
00364 case actuationY:
00365
00366 resp = generateActuation( resp, actData.ETMYcal / actData.length,
00367 actData.pendFY, actData.pendQY);
00368 break;
00369 }
00370
00371
00372
00373
00374
00375
00376
00377
00378 memset( &detector, 0, sizeof( DetectorResponse ) );
00379 detector.site = (LALDetector *) LALMalloc( sizeof(LALDetector) );
00380 XLALReturnDetector( detector.site, ifoNumber );
00381 detector.transfer = XLALCreateCOMPLEX8FrequencySeries( chan->name,
00382 &(chan->epoch), 0, 1.0 / ( duration ), &strainPerCount,
00383 ( sampleRate * duration / 2 + 1 ) );
00384
00385 XLALUnitInvert( &(detector.transfer->sampleUnits), &(resp->sampleUnits) );
00386
00387
00388 unity = XLALCreateCOMPLEX8Vector( resp->data->length );
00389 for ( k = 0; k < unity->length; ++k )
00390 {
00391 unity->data[k].re = 1.0;
00392 unity->data[k].im = 0.0;
00393 }
00394
00395 XLALCCVectorDivide( detector.transfer->data, unity, resp->data );
00396 XLALDestroyCOMPLEX8Vector( unity );
00397
00398
00399 waveformStartTime = XLALGPStoINT8( &(inspInj->geocent_end_time));
00400 waveformStartTime -= (INT8) ( 1000000000.0 * ppnParams.tc );
00401
00402 XLALINT8toGPS( &(waveform.a->epoch), waveformStartTime );
00403 memcpy(&(waveform.f->epoch), &(waveform.a->epoch),
00404 sizeof(LIGOTimeGPS) );
00405 memcpy(&(waveform.phi->epoch), &(waveform.a->epoch),
00406 sizeof(LIGOTimeGPS) );
00407
00408
00409 LAL_CALL( LALSimulateCoherentGW( status, chan, &waveform, &detector ),
00410 status);
00411
00412 destroyCoherentGW( &waveform );
00413 XLALDestroyCOMPLEX8FrequencySeries( detector.transfer );
00414 if ( detector.site ) LALFree( detector.site );
00415
00416 XLALDestroyCOMPLEX8FrequencySeries( resp );
00417 return( chan );
00418 }
00419
00420
00421
00422
00423
00424
00425
00426 int main( int argc, char *argv[] )
00427 {
00428 LALStatus status = blank_status;
00429 LALLeapSecAccuracy accuracy = LALLEAPSEC_LOOSE;
00430 LIGOTimeGPS gpsStartTime = {0, 0};
00431 LIGOTimeGPS earliestEndTime = {0, 0};
00432 ResponseFunction injectionResponse = noResponse;
00433
00434
00435 RandomParams *randParams = NULL;
00436 REAL4 massPar = 0;
00437 MassDistribution mDist = logComponentMass;
00438 InterferometerNumber ifoNumber = LAL_UNKNOWN_IFO;
00439 REAL4 desiredSnr = 0;
00440 REAL4 snrsqAt1Mpc = 0;
00441 REAL4TimeSeries *chan = NULL;
00442 REAL4FFTPlan *pfwd;
00443 COMPLEX8FrequencySeries *fftData;
00444
00445
00446 CHAR waveform[LIGOMETA_WAVEFORM_MAX];
00447
00448
00449 CHAR fname[FILENAME_MAX];
00450 MetadataTable proctable;
00451 MetadataTable procparams;
00452 MetadataTable inspInjections;
00453 MetadataTable ringInjections;
00454 ProcessParamsTable *this_proc_param;
00455 SimInspiralTable *inj = NULL;
00456 SimRingdownTable *ringList = NULL;
00457 LIGOLwXMLStream xmlfp;
00458 FILE *fp = NULL;
00459
00460
00461 struct option long_options[] =
00462 {
00463 {"help", no_argument, 0, 'h'},
00464 {"verbose", no_argument, &vrbflg, 1 },
00465 {"version", no_argument, 0, 'V'},
00466 {"gps-start-time", required_argument, 0, 'a'},
00467 {"injection-type", required_argument, 0, 't'},
00468 {"seed", required_argument, 0, 's'},
00469 {"debug-level", required_argument, 0, 'z'},
00470 {0, 0, 0, 0}
00471 };
00472 int c;
00473
00474
00475
00476
00477 actuationParams[LAL_IFO_H1].ETMXcal = -0.795e-9;
00478 actuationParams[LAL_IFO_H1].pendFX = 0.767;
00479 actuationParams[LAL_IFO_H1].pendQX = 10.0;
00480 actuationParams[LAL_IFO_H1].ETMYcal = -0.827e-9;
00481 actuationParams[LAL_IFO_H1].pendFY = 0.761;
00482 actuationParams[LAL_IFO_H1].pendQY = 10.0;
00483 actuationParams[LAL_IFO_H1].length = 4000.0;
00484
00485
00486
00487 actuationParams[LAL_IFO_H2].ETMXcal = -0.876e-9;
00488 actuationParams[LAL_IFO_H2].pendFX = 0.749;
00489 actuationParams[LAL_IFO_H2].pendQX = 10.0;
00490 actuationParams[LAL_IFO_H2].ETMYcal = -0.912e-9;
00491 actuationParams[LAL_IFO_H2].pendFY = 0.764;
00492 actuationParams[LAL_IFO_H2].pendQY = 10.0;
00493 actuationParams[LAL_IFO_H2].length = 2000.0;
00494
00495
00496
00497 actuationParams[LAL_IFO_L1].ETMXcal = -0.447e-9;
00498 actuationParams[LAL_IFO_L1].pendFX = 0.766;
00499 actuationParams[LAL_IFO_L1].pendQX = 100.0;
00500 actuationParams[LAL_IFO_L1].ETMYcal = -0.438e-9;
00501 actuationParams[LAL_IFO_L1].pendFY = 0.756;
00502 actuationParams[LAL_IFO_L1].pendQY = 100.0;
00503 actuationParams[LAL_IFO_L1].length = 4000.0;
00504
00505
00506 lal_errhandler = LAL_ERR_EXIT;
00507 set_debug_level( "33" );
00508
00509
00510
00511 proctable.processTable = (ProcessTable *)
00512 calloc( 1, sizeof(ProcessTable) );
00513 LAL_CALL( LALGPSTimeNow ( &status, &(proctable.processTable->start_time),
00514 &accuracy ), &status );
00515 LAL_CALL( populate_process_table( &status, proctable.processTable,
00516 PROGRAM_NAME, CVS_REVISION, CVS_SOURCE, CVS_DATE ), &status );
00517 LALSnprintf( proctable.processTable->comment, LIGOMETA_COMMENT_MAX, " " );
00518 this_proc_param = procparams.processParamsTable = (ProcessParamsTable *)
00519 calloc( 1, sizeof(ProcessParamsTable) );
00520
00521
00522 memset( waveform, 0, LIGOMETA_WAVEFORM_MAX * sizeof(CHAR) );
00523 LALSnprintf( waveform, LIGOMETA_WAVEFORM_MAX * sizeof(CHAR),
00524 "SpinTaylorthreePN");
00525
00526
00527
00528
00529
00530
00531
00532
00533 while ( 1 )
00534 {
00535
00536 int option_index = 0;
00537 long int gpsinput;
00538
00539 c = getopt_long_only( argc, argv, "a:hs:t:V", long_options, &option_index );
00540
00541
00542 if ( c == - 1 )
00543 {
00544 break;
00545 }
00546
00547 switch ( c )
00548 {
00549 case 0:
00550
00551 if ( long_options[option_index].flag != 0 )
00552 {
00553 break;
00554 }
00555 else
00556 {
00557 fprintf( stderr, "error parsing option %s with argument %s\n",
00558 long_options[option_index].name, optarg );
00559 exit( 1 );
00560 }
00561 break;
00562
00563 case 'a':
00564 gpsinput = atol( optarg );
00565 if ( gpsinput < 441417609 )
00566 {
00567 fprintf( stderr, "invalid argument to --%s:\n"
00568 "GPS start time is prior to "
00569 "Jan 01, 1994 00:00:00 UTC:\n"
00570 "(%ld specified)\n",
00571 long_options[option_index].name, gpsinput );
00572 exit( 1 );
00573 }
00574 if ( gpsinput > 999999999 )
00575 {
00576 fprintf( stderr, "invalid argument to --%s:\n"
00577 "GPS start time is after "
00578 "Sep 14, 2011 01:46:26 UTC:\n"
00579 "(%ld specified)\n",
00580 long_options[option_index].name, gpsinput );
00581 exit( 1 );
00582 }
00583 gpsStartTime.gpsSeconds = gpsinput;
00584
00585 this_proc_param = this_proc_param->next =
00586 next_process_param( long_options[option_index].name, "int",
00587 "%ld", gpsinput );
00588 break;
00589
00590 case 's':
00591 randSeed = atoi( optarg );
00592 this_proc_param = this_proc_param->next =
00593 next_process_param( long_options[option_index].name, "int",
00594 "%d", randSeed );
00595 break;
00596
00597 case 't':
00598
00599 if ( ! strcmp( "strain", optarg ) )
00600 {
00601 injectionResponse = unityResponse;
00602 }
00603 else if ( ! strcmp( "etmx", optarg ) )
00604 {
00605 injectionResponse = actuationX;
00606 }
00607 else if ( ! strcmp( "etmy", optarg ) )
00608 {
00609 injectionResponse = actuationY;
00610 }
00611 else
00612 {
00613 fprintf( stderr, "invalid argument to --%s:\n"
00614 "unknown injection type specified: "
00615 "%s (must be strain, etmx or etmy)\n",
00616 long_options[option_index].name, optarg );
00617 exit( 1 );
00618 }
00619 next_process_param( long_options[option_index].name, "string", "%s",
00620 optarg );
00621 break;
00622
00623 case 'V':
00624
00625 fprintf( stdout, "blind hardware injection generation routine\n"
00626 "Stephen Fairhurst\n"
00627 "CVS Version: " CVS_ID_STRING "\n"
00628 "CVS Tag: " CVS_NAME_STRING "\n" );
00629 exit( 0 );
00630 break;
00631
00632 case 'z':
00633 set_debug_level( optarg );
00634 next_process_param( long_options[option_index].name, "int", "%d",
00635 optarg );
00636 break;
00637
00638 case 'h':
00639 case '?':
00640 fprintf( stderr, USAGE );
00641 exit( 0 );
00642 break;
00643
00644 default:
00645 fprintf( stderr, "unknown error while parsing options\n" );
00646 fprintf( stderr, USAGE );
00647 exit( 1 );
00648 }
00649 }
00650
00651 if ( optind < argc )
00652 {
00653 fprintf( stderr, "extraneous command line arguments:\n" );
00654 while ( optind < argc )
00655 {
00656 fprintf ( stderr, "%s\n", argv[optind++] );
00657 }
00658 exit( 1 );
00659 }
00660
00661
00662 if ( injectionResponse == noResponse )
00663 {
00664 fprintf( stderr, "Must specify the --injection-type\n" );
00665 exit( 1 );
00666 }
00667
00668 if ( gpsStartTime.gpsSeconds <= 0 )
00669 {
00670 fprintf( stderr, "Must specify the --gps-start-time\n" );
00671 exit( 1 );
00672 }
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682 randParams = XLALCreateRandomParams( randSeed );
00683
00684
00685 memset( &inspInjections, 0, sizeof(MetadataTable) );
00686 memset( &ringInjections, 0, sizeof(MetadataTable) );
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696 inspInjections.simInspiralTable = inj = (SimInspiralTable *)
00697 LALCalloc( 1, sizeof(SimInspiralTable) );
00698
00699 ringInjections.simRingdownTable = ringList = (SimRingdownTable *)
00700 LALCalloc( 1, sizeof(SimRingdownTable) );
00701
00702
00703
00704 earliestEndTime = gpsStartTime;
00705 earliestEndTime = *XLALGPSAdd( &earliestEndTime, longestSignal );
00706 inj = XLALRandomInspiralTime( inj, randParams, earliestEndTime,
00707 timeWindow );
00708
00709
00710 inj->distance = 1.0;
00711
00712
00713 massPar = XLALUniformDeviate( randParams );
00714
00715 if ( vrbflg) fprintf( stdout,
00716 "Random variable to determine inj type = %f\n",
00717 massPar);
00718
00719 while ( numInjections == 0 )
00720 {
00721
00722 normalDev = XLALCreateREAL4Vector( 1 );
00723 XLALNormalDeviates(normalDev, randParams);
00724
00725 if ( massPar < BNSfrac )
00726 {
00727 if ( vrbflg ) fprintf( stdout, "Generating a BNS injection\n" );
00728 inj = XLALRandomInspiralMasses( inj, randParams, mDist,
00729 minNSMass, maxNSMass, minNSMass, maxNSMass, minTotalMass, maxTotalMass );
00730 inj = XLALRandomInspiralSpins( inj, randParams, minNSSpin,
00731 maxNSSpin, minNSSpin, maxNSSpin, -1.0, 1.0, 0.0, 0.1 );
00732 desiredSnr = bnsSnrMean + bnsSnrStd * normalDev->data[0];
00733 }
00734 else if ( massPar < (BNSfrac + BBHfrac) )
00735 {
00736 if ( vrbflg ) fprintf( stdout, "Generating a BBH injection\n" );
00737 inj = XLALRandomInspiralMasses( inj, randParams, mDist,
00738 minBHMass, maxBHMass, minBHMass, maxBHMass, minTotalMass, maxTotalMass );
00739 inj = XLALRandomInspiralSpins( inj, randParams, minBHSpin,
00740 maxBHSpin, minBHSpin, maxBHSpin , -1.0, 1.0, 0.0, 0.1);
00741 desiredSnr = snrMean + snrStd * normalDev->data[0];
00742 }
00743 else
00744 {
00745 if ( vrbflg ) fprintf( stdout, "Generating an NS - BH injection\n" );
00746 inj = XLALRandomInspiralMasses( inj, randParams, mDist,
00747 minNSMass, maxNSMass, minBHMass, maxBHMass, minTotalMass, maxTotalMass );
00748 inj = XLALRandomInspiralSpins( inj, randParams, minNSSpin,
00749 maxNSSpin, minBHSpin, maxBHSpin , -1.0, 1.0, 0.0, 0.1);
00750 desiredSnr = snrMean + snrStd * normalDev->data[0];
00751 }
00752 XLALDestroyVector( normalDev );
00753
00754
00755
00756 InclDistribution iDist = uniformInclDist ;
00757 inj = XLALRandomInspiralSkyLocation( inj, randParams );
00758 inj = XLALRandomInspiralOrientation( inj, randParams, iDist, 0);
00759
00760
00761 LALSnprintf( inj->source, LIGOMETA_SOURCE_MAX * sizeof(CHAR), "???" );
00762 memcpy( inj->waveform, waveform, LIGOMETA_WAVEFORM_MAX * sizeof(CHAR));
00763
00764
00765 inj = XLALPopulateSimInspiralSiteInfo( inj );
00766
00767
00768 inj->f_lower = fLower;
00769
00770
00771
00772
00773
00774
00775
00776
00777 if ( vrbflg ) fprintf( stdout, "Injection details\n"
00778 "mass 1 = %.2f;\t spin 1 = (%.2f, %.2f, %.2f)\n"
00779 "mass 2 = %.2f;\t spin 2 = (%.2f, %.2f, %.2f)\n"
00780 "Hanford effective distance = %.2f Mpc\n"
00781 "Livingston effective distance = %.2f Mpc\n",
00782 inj->mass1, inj->spin1x, inj->spin1y, inj->spin1z,
00783 inj->mass2, inj->spin2x, inj->spin2y, inj->spin2z,
00784 inj->eff_dist_h,
00785 inj->eff_dist_l );
00786
00787 for ( ifoNumber = 0; ifoNumber < LAL_NUM_IFO; ifoNumber++ )
00788 {
00789 if ( ifoNumber == LAL_IFO_H1 || ifoNumber == LAL_IFO_L1 )
00790 {
00791 ResponseFunction responseType = unityResponse;
</