Commit 082d4a7a authored by Cristina Valeria Torres's avatar Cristina Valeria Torres
Browse files

Merge branch 'master' of git+ssh://cristina.torres@ligo-vcs.phys.uwm.edu/usr/local/git/lalsuite

Original: bde4cd1bd855947735bfa5ac3530e42a401c2e1b
parents 55c89fcc 39beb93d
......@@ -811,22 +811,35 @@ for result_db in result_dbs_cache:
############################################################################
# Compute the upper limits as a function of mass for these types of mass binnings
if 'FULL_DATA' in tag: #for full data, use FAR of loudest event
for bintype in ["MASS1_MASS2","CHIRP_MASS","TOTAL_MASS"]:
output_cache_name = ''.join(sorted(instruments)) + "-SEARCH_VOLUME_BY_"+bintype+"_" + tag + ".cache"
search_volume_node = inspiral.SearchVolumeNode( search_volume_job, result_db.path(), output_tag = tag, output_cache = output_cache_name, veto_segments_name = "VETO_CAT"+str(get_veto_cat_from_tag(tag))+"_CUMULATIVE", bintype = bintype )
if 'FULL_DATA' in tag:
#for full data, use FAR of loudest event
output_cache_name = ''.join(sorted(instruments)) + "-SEARCH_VOLUME_" + tag + ".cache"
search_volume_node = inspiral.SearchVolumeNode( search_volume_job )
search_volume_node.add_database( result_db.path() )
search_volume_node.set_output_tag( tag )
search_volume_node.set_output_cache( output_cache_name )
search_volume_node.set_veto_segments_name( "VETO_CAT"+str(get_veto_cat_from_tag(tag))+"_CUMULATIVE" )
search_volume_node.add_parent( ccfar_node )
search_upper_limit_node = inspiral.SearchUpperLimitNode( search_upper_limit_job, output_cache_name )
search_upper_limit_node = inspiral.SearchUpperLimitNode( search_upper_limit_job )
search_upper_limit_node.set_open_box()
search_upper_limit_node.add_input_cache( output_cache_name )
search_upper_limit_node.add_parent( search_volume_node )
dag.add_node( search_volume_node )
dag.add_node( search_upper_limit_node )
# upper limits for playground (using FAR of expected loudest event)
# make corresponding jobs for playground (using FAR of expected loudest event)
search_volume_node = inspiral.SearchVolumeNode( search_volume_job )
search_volume_node.use_expected_loudest_event()
output_cache_name = output_cache_name.replace("FULL_DATA","PLAYGROUND")
search_volume_node = inspiral.SearchVolumeNode( search_volume_job, result_db.path(), output_tag = tag.replace("FULL_DATA","PLAYGROUND"), output_cache = output_cache_name, use_expected_loudest_event = True, veto_segments_name = "VETO_CAT"+str(get_veto_cat_from_tag(tag))+"_CUMULATIVE", bintype = bintype )
search_volume_node.add_database( result_db.path() )
search_volume_node.set_output_tag( tag.replace("FULL_DATA","PLAYGROUND" ) )
search_volume_node.set_output_cache( output_cache_name )
search_volume_node.set_veto_segments_name( "VETO_CAT"+str(get_veto_cat_from_tag(tag))+"_CUMULATIVE" )
search_volume_node.add_parent( ccfar_node )
search_upper_limit_node = inspiral.SearchUpperLimitNode( search_upper_limit_job, output_cache_name )
search_upper_limit_node = inspiral.SearchUpperLimitNode( search_upper_limit_job )
search_upper_limit_node.add_input_cache( output_cache_name )
search_upper_limit_node.add_parent( search_volume_node )
dag.add_node( search_volume_node )
dag.add_node( search_upper_limit_node )
......
......@@ -2733,7 +2733,7 @@ int main( int argc, char *argv[] )
simTable->bandpass = bandPassInj;
/* populate the sim_ringdown table */
/* populate the sim_ringdown table */
if ( writeSimRing )
{
memcpy( simRingTable->waveform, "Ringdown",
......@@ -2752,14 +2752,17 @@ int main( int argc, char *argv[] )
simRingTable->polarization = simTable->polarization;
simRingTable->phase = 0;
simRingTable->mass = XLALNonSpinBinaryFinalBHMass(simTable->eta, simTable->mass1, simTable->mass2);
simRingTable->spin = XLALNonSpinBinaryFinalBHSpin(simTable->eta);
/* The final spin calc has been generalized so as to allow initially spinning systems*/
/* simRingTable->spin = XLALNonSpinBinaryFinalBHSpin(simTable->eta); */
simRingTable->spin = XLALSpinBinaryFinalBHSpin(simTable->eta, simTable->mass1, simTable->mass2,
simTable->spin1x, simTable->spin2x,simTable->spin1y, simTable->spin2y, simTable->spin1z, simTable->spin2z);
simRingTable->frequency = XLALBlackHoleRingFrequency( simRingTable->mass, simRingTable->spin);
simRingTable->quality = XLALBlackHoleRingQuality(simRingTable->spin);
simRingTable->epsilon = 0.01;
simRingTable->epsilon = 0.01;
simRingTable->amplitude = XLALBlackHoleRingAmplitude( simRingTable->frequency, simRingTable->quality, simRingTable->distance, simRingTable->epsilon );
simRingTable->eff_dist_h = simTable->eff_dist_h;
simRingTable->eff_dist_l = simTable->eff_dist_l;
simRingTable->eff_dist_v = simTable->eff_dist_v;
simRingTable->eff_dist_h = simTable->eff_dist_h;
simRingTable->eff_dist_l = simTable->eff_dist_l;
simRingTable->eff_dist_v = simTable->eff_dist_v;
simRingTable->hrss = XLALBlackHoleRingHRSS( simRingTable->frequency, simRingTable->quality, simRingTable->amplitude, 2., 0. );
// need hplus & hcross in each detector to populate these
simRingTable->hrss_h = 0.; //XLALBlackHoleRingHRSS( simRingTable->frequency, simRingTable->quality, simRingTable->amplitude, 0., 0. );
......
......@@ -3347,30 +3347,26 @@ class SearchVolumeNode(pipeline.SqliteNode):
"""
A search volume node.
"""
def __init__(self, job, database, output_cache = None, output_tag = "SEARCH_VOLUME", bootstrap_iterations=10000, veto_segments_name="vetoes", use_expected_loudest_event = False, bintype = "TOTAL_MASS"):
def __init__(self, job):
"""
@database: the pipedown database containing the injection triggers
@ouptut_cache: name prefix for cache file to be written out by program
@output_tag: a string label for the output files
@use_expected_loudest_event: disables the use of loudest event FAR, use 1./livetime instead
"""
pipeline.SqliteNode.__init__(self, job)
self.add_var_arg(database)
self.add_var_opt("bootstrap-iterations",bootstrap_iterations)
self.add_var_opt("veto-segments-name",veto_segments_name)
if output_cache:
self.add_var_opt("output-cache",output_cache)
if output_tag:
self.add_var_opt("output-name-tag",output_tag)
if use_expected_loudest_event:
self.add_var_arg("--use-expected-loudest-event")
if bintype == "TOTAL_MASS":
self.add_var_arg("--bin-by-total-mass")
if bintype == "CHIRP_MASS":
self.add_var_arg("--bin-by-chirp-mass")
if bintype == "MASS1_MASS2":
self.add_var_arg("--bin-by-m1m2")
def add_database(self, db):
self.add_var_arg(db)
def set_output_cache(self, file):
self.add_var_opt("output-cache", file)
def set_output_tag(self, tag):
self.add_var_opt("output-name-tag",tag)
def set_veto_segments_name(self, name):
self.add_var_opt("veto-segments-name", name)
def use_expected_loudest_event(self):
self.add_var_arg("--use-expected-loudest-event")
class SearchUpperLimitJob(pipeline.SqliteJob):
"""
......@@ -3390,14 +3386,16 @@ class SearchUpperLimitNode(pipeline.SqliteNode):
"""
A search upper limit node.
"""
def __init__(self, job, input_cache):
def __init__(self, job):
"""
@job: a SearchUpperLimitJob
"""
pipeline.SqliteNode.__init__(self, job)
self.add_var_opt("input-cache", input_cache)
self.open_box = False
def add_input_cache(self, input_cache):
self.add_var_arg(input_cache)
def set_open_box(self):
'''
Set the open box flag.
......@@ -3405,4 +3403,3 @@ class SearchUpperLimitNode(pipeline.SqliteNode):
if not self.open_box:
self.open_box = True
self.add_var_arg("--open-box")
......@@ -979,6 +979,11 @@ plot-normalized-energy-range = 0,25.5
; THIS SECTION MUST BE LEFT BLANK!!!!!!!!
[search-volume]
dist-err = 0.197
dist-sys-err = 0.074
mass-bins = 5
dist-err = 0.2
dist-sys-err = 0.1
mass-bins = 11
bin-by-total-mass =
bin-by-component-mass =
bin-by-chirp-mass =
iterations = 10000
bootstrap =
......@@ -930,6 +930,12 @@ plot-normalized-energy-range = 0,25.5
; THIS SECTION MUST BE LEFT BLANK!!!!!!!!
[search-volume]
dist-err = 0.197
dist-sys-err = 0.074
mass-bins = 5
dist-err = 0.2
dist-sys-err = 0.1
total-mass-bins = '2,5,8,11,14,17,20,25'
component-mass1-bins = '1,3,8,13,18,23'
bin-by-chirp-mass =
bin-by-bns-bbh =
iterations = 10000
bootstrap =
......@@ -904,49 +904,14 @@ def write_upperlimit(page, opts,thisSearch='playground',pipedown='pipedown'):
sections.append(section)
htmlAndPlot.append(section)
titles[section] = "Sensitive Distance as a Function of Component Masses."
tags[section] = '*lalapps_cbc_sink_by_mass1_mass2*'+ifartag[0]+'*CAT_4_VETO*'
tags[section] = '*lalapps_cbc_sink_by_component_mass*'+ifartag[0]+'*CAT_4_VETO*'
htmltags[section] = '*'+ifartag[0]+'_CAT_4_VETO_mass1_mass2_range_summary*'
imagetags[section] = ['distance','fractional_error']
captions[section] = """Mean sensitive distance for the CBC search pipeline as a function of the component masses. The sensitive distance is determined by computing the cube root of the mean distance^3 for injections found with a FAR below the FAR of the loudest event. The corresponding fractional error plot gives a measure of the relative uncertainty in the measured sensitive distance arising from uncertainty in calibration and a finite injection population. (CAT2 and CAT3 vetoes applied and CBC hardware injections removed)."""
comments[section] = None
configs[section] = None
images_dirs[section] = PDdir
section= 'upperlimitVmtotal'
sections.append(section)
htmlAndPlot.append(section)
titles[section] = "Upper Limits as a Function of Total Mass."
tags[section] = '*lalapps_cbc_sink_by_total_mass*'+ifartag[0]+'*CAT_4_VETO*'
htmltags[section] = '*'+ifartag[0]+'_CAT_4_VETO_total_mass_range_summary*'
imagetags[section] = ['upper_limit_plot']
captions[section] = """Uncombined upper limits on the inspiral merger rate (in mergers/Mpc^3/yr) as a function of total mass. The plot gives the upper limits with and without including the effect of systematic errors due to calibration uncertainty and finite injection population. (CAT2 and CAT3 vetoes applied and CBC hardware injections removed)."""
comments[section] = None
configs[section] = None
images_dirs[section] = PDdir
section= 'upperlimitVmchirp'
sections.append(section)
htmlAndPlot.append(section)
titles[section] = "Upper Limits as a Function of Chirp Mass."
tags[section] = '*lalapps_cbc_sink_by_chirp_mass*'+ifartag[0]+'*CAT_4_VETO*'
htmltags[section] = '*'+ifartag[0]+'_CAT_4_VETO_chirp_mass_range_summary*'
imagetags[section] = ['upper_limit_plot']
captions[section] = """Uncombined upper limits on the inspiral merger rate (in mergers/Mpc^3/yr) as a function of chirp mass. The plot gives the upper limits with and without including the effect of systematic errors due to calibration uncertainty and finite injection population. (CAT2 and CAT3 vetoes applied and CBC hardware injections removed)."""
imagetags[section] = ['distance']
captions[section] = """Mean sensitive distance for the CBC search pipeline as a function of the component masses. The sensitive distance is determined by computing the cube root of the mean distance^3 for injections found with a FAR below the FAR of the loudest event. (CAT2 and CAT3 vetoes applied and CBC hardware injections removed)."""
comments[section] = None
configs[section] = None
images_dirs[section] = PDdir
section= 'upperlimitVm1m2'
sections.append(section)
htmlAndPlot.append(section)
titles[section] = "Upper Limits as a Function of Component Masses."
tags[section] = '*lalapps_cbc_sink_by_mass1_mass2*'+ifartag[0]+'*CAT_4_VETO*'
htmltags[section] = '*'+ifartag[0]+'_CAT_4_VETO_mass1_mass2_range_summary*'
imagetags[section] = ['upper_limit_plot']
captions[section] = """Uncombined upper limits on the inspiral merger rate (in mergers/Mpc^3/yr) as a function of the component masses. The plot gives the upper limits with and without including the effect of systematic errors due to calibration uncertainty and finite injection population. (CAT2 and CAT3 vetoes applied and CBC hardware injections removed)."""
comments[section] = None
configs[section] = None
images_dirs[section] = PDdir
section= 'combupperlimitVmtotal'
sections.append(section)
......@@ -976,8 +941,8 @@ def write_upperlimit(page, opts,thisSearch='playground',pipedown='pipedown'):
sections.append(section)
htmlAndPlot.append(section)
titles[section] = "Combined Upper Limits as a Function of Component Masses."
tags[section] = '*lalapps_cbc_combined_sink_by_mass1_mass2*'+ifartag[0]+'*CAT_4_VETO*'
htmltags[section] = '*'+ifartag[0]+'_CAT_4_VETO_mass1_mass2_combined_upper_limit*'
tags[section] = '*lalapps_cbc_combined_sink_by_component_mass*'+ifartag[0]+'*CAT_4_VETO*'
htmltags[section] = '*'+ifartag[0]+'_CAT_4_VETO_component_mass_combined_upper_limit*'
imagetags[section] = ['upper_limit_plot']
captions[section] = """Combined upper limits on the inspiral merger rate (in mergers/Mpc^3/yr) as a function of the component masses. The plot gives the upper limits with and without including the effect of systematic errors due to calibration uncertainty and finite injection population. (CAT2 and CAT3 vetoes applied and CBC hardware injections removed)."""
comments[section] = None
......@@ -2784,9 +2749,9 @@ if opts.hardware_injection is True:
if opts.full_data is True: html_sections['full_data'] = "Full Data"
if opts.followup is True: html_sections['followup'] = "Followup"
if opts.upperlimit is True:
html_sections['upperlimit_play'] = "Playground Upper Limit"
html_sections['upperlimit_play'] = "Search Sensitivity"
if opts.upperlimit and opts.full_data:
html_sections['upperlimit_full'] = "Full Data Upper Limit"
html_sections['upperlimit_full'] = "Upper Limits"
if opts.summary_files:
html_sections['summary_files'] = "Summary files"
html_sections['logfile'] = "Log File"
......@@ -2836,7 +2801,7 @@ mainpage.div(id_='header')
mainpage.add('<h1>' + title +' </h1>')
leads = 'Lead: ' + Lead + '&nbsp;&nbsp&nbsp Second: '+ Second
mainpage.add('<h3> ' + leads + ' </h3>')
mainpage.add('<h3>' + Notes + ' </h3')
mainpage.add('<h3>' + Notes + ' </h3>')
mainpage.div.close()
mainpage.div(id_='iframecontent')
mainpage.add('<p id="placeholder">Please select a report section on the left.</p>')
......
......@@ -186,6 +186,10 @@ int main(int argc, char **argv){
inputParams_t inputs;
/* allocate memory for headers structures */
/* head1 = calloc(sizeof(headerRecord1), 1);
head2 = calloc(sizeof(headerRecord2), 1); */
/* get input parameters */
get_input_args(&inputs, argc, argv);
......@@ -286,8 +290,6 @@ inputs.noverlap);
return 1;
}
fprintf(stderr, "sizeof(coeffArray) = %zu\n", sizeof(coeffArray));
if(verbose){
fprintf(stderr, "The array size for ephemeris file %s is %d.\n",
inputs.ephemfile , ARRAY_SIZE);
......@@ -490,25 +492,56 @@ A[1], A[2]);
/* function to determine the record size and initialise the ephemeris */
INT4 fsizer(FILE *fp){
UINT4 mrecl=0;
INT4 kmx=0, khi=0;
INT4 i=0;
UINT4 i=0;
INT4 nd=0;
UINT4 size=0;
mrecl=1500;
/* read in header info */
if( fread(&head1, sizeof(REAL8), mrecl, fp) != mrecl ){
fprintf(stderr, "Error reading in ephemeris header info.\n");
/* read in data structure one part at a time so 32vs64bit alignment
issues don't become a problem */
if( fread(&head1.data.ttl, sizeof(head1.data.ttl), 1, fp) != 1 ){
fprintf(stderr, "Error reading in ephemeris ttl header info.\n");
return 0;
}
if( fread(&head1.data.cnam, sizeof(head1.data.cnam), 1, fp) != 1 ){
fprintf(stderr, "Error reading in ephemeris cnam header info.\n");
return 0;
}
if( fread(&head1.data.ss, sizeof(head1.data.ss), 1, fp) != 1 ){
fprintf(stderr, "Error reading in ephemeris ss header info.\n");
return 0;
}
if( fread(&head1.data.con, sizeof(head1.data.con), 1, fp) != 1 ){
fprintf(stderr, "Error reading in ephemeris con header info.\n");
return 0;
}
if( fread(&head1.data.au, sizeof(head1.data.au), 1, fp) != 1 ){
fprintf(stderr, "Error reading in ephemeris au header info.\n");
return 0;
}
if( fread(&head1.data.emrat, sizeof(head1.data.emrat), 1, fp) != 1 ){
fprintf(stderr, "Error reading in ephemeris emrat header info.\n");
return 0;
}
if( fread(&head1.data.ipt, sizeof(head1.data.ipt), 1, fp) != 1 ){
fprintf(stderr, "Error reading in ephemeris ipt header info.\n");
return 0;
}
if( fread(&head1.data.numde, sizeof(head1.data.numde), 1, fp) != 1 ){
fprintf(stderr, "Error reading in ephemeris numde header info.\n");
return 0;
}
if( fread(&head1.data.libratPtr, sizeof(head1.data.libratPtr), 1, fp) != 1 ){
fprintf(stderr, "Error reading in ephemeris libratPtr header info.\n");
return 0;
}
/* flip bytes of ipt */
endian_swap((CHAR*)&head1.data.ipt, sizeof(INT4), 36);
endian_swap((CHAR*)&head1.data.libratPtr, sizeof(INT4), 3);
/*** Calculate array size in the ephemeris */
for (i=0;i<12;i++){
if(head1.data.ipt[i][0] > kmx){
......@@ -536,28 +569,25 @@ INT4 fsizer(FILE *fp){
}
/*****************************************/
rewind(fp);
/* set file pointer to position of header data 2 */
fseek(fp, size*sizeof(REAL8), SEEK_SET);
/* Initialise to ephemeris to the point at which the coefficient values start
*/
if( fread(&head1, sizeof(REAL8), size, fp) != size ){
fprintf(stderr, "Error initialising ephemeris.\n");
return 0;
}
if( fread(&head2, sizeof(REAL8), size, fp) != size ){
fprintf(stderr, "Error initialising ephemeris.\n");
/* read in values of head2 */
if( fread(&head2.data.cval, sizeof(head2.data.cval), 1, fp) != 1 ){
fprintf(stderr, "Error reading in ephemeris cval header info.\n");
return 0;
}
/* Initialise to ephemeris to the point at which the coefficient values start
*/
fseek(fp, 2*size*sizeof(REAL8), SEEK_SET);
/* flip bytes of values */
endian_swap((CHAR*)&head1.data.ipt, sizeof(INT4), 36);
endian_swap((CHAR*)&head1.data.libratPtr, sizeof(INT4), 3);
endian_swap((CHAR*)&head1.data.au, sizeof(REAL8), 1);
endian_swap((CHAR*)&head1.data.emrat, sizeof(REAL8), 1);
endian_swap((CHAR*)&head1.data.ss, sizeof(REAL8), 3);
endian_swap((CHAR*)&head1.data.con, sizeof(INT4), 1);
endian_swap((CHAR*)&head1.data.numde, sizeof(INT4), 1);
endian_swap((CHAR*)&head2.data.cval, sizeof(REAL8), 400);
if(verbose){
......
......@@ -1503,7 +1503,7 @@ void AddSFTtoList(LALStatus *status,
*sftHead = sftList;
}
else {
(*sftTail)->nextSFT = (struct SFTListElement *)sftList;
(*sftTail)->nextSFT = (SFTListElement *)sftList;
}
*sftTail = sftList;
......@@ -1531,7 +1531,7 @@ void AddPSDtoList(LALStatus *status,
if (!(*psdHead)) {
*psdHead = psdList;
} else {
(*psdTail)->nextPSD = (struct PSDListElement *)psdList;
(*psdTail)->nextPSD = (PSDListElement *)psdList;
}
*psdTail = psdList;
......@@ -1558,7 +1558,7 @@ void AddREAL8toList(LALStatus *status,
if (!(*head)) {
*head = List;
} else {
(*tail)->nextVal = (struct REAL8ListElement *)List;
(*tail)->nextVal = (REAL8ListElement *)List;
}
*tail = List;
......@@ -1587,7 +1587,7 @@ void AddBeamFntoList(LALStatus *status,
if (!(*beamHead)) {
*beamHead = beamList;
} else {
(*beamTail)->nextBeamfn = (struct CrossCorrBeamFnListElement *)beamList;
(*beamTail)->nextBeamfn = (CrossCorrBeamFnListElement *)beamList;
}
*beamTail = beamList;
......
......@@ -776,9 +776,9 @@ int main( int argc, char **argv )
for( ifoNumber = 0; ifoNumber < LAL_NUM_IFO; ifoNumber++)
{
/* get location in three dimensions */
for ( i = 0; i < 3; i++ )
for ( ui = 0; ui < 3; ui++ )
{
detLoc[i] = (double) detectors[ifoNumber]->location[i];
detLoc[ui] = (double) detectors[ifoNumber]->location[ui];
}
/* calculate time offsets */
timeOffsets[ifoNumber] =
......@@ -817,7 +817,7 @@ int main( int argc, char **argv )
pValues, gammaBeta, snrComps,
nullSNR, traceSNR, bankVeto,
autoVeto, chiSquare, PTFM );
verbose( "Generated triggers for segment %d, template %d, sky points %d at %ld \n", j, i, sp, timeval_subtract(&startTime) );
verbose( "Generated triggers for segment %d, template %d, sky point %d at %ld \n", j, i, sp, timeval_subtract(&startTime) );
// Clustering can happen here. The clustering routine needs refining
// coh_PTF_cluster_triggers(params,&eventList,&thisEvent);
......@@ -864,11 +864,14 @@ int main( int argc, char **argv )
}
}
LALFree(chisqOverlaps);
chisqOverlaps = NULL;
}
if (frequencyRangesPlus)
LALFree(frequencyRangesPlus);
frequencyRangesPlus = NULL;
if (frequencyRangesCross)
LALFree(frequencyRangesCross);
frequencyRangesCross = NULL;
}
if ( params->doBankVeto )
{
......@@ -1631,7 +1634,6 @@ void coh_PTF_statistic(
tempqVec = XLALCreateCOMPLEX8VectorSequence ( 1, numPoints );
if (! chisqOverlaps)
{
fprintf( stderr,"Calculating chi2 filters\n" );
chisqOverlaps = LALCalloc(params->numChiSquareBins,sizeof( *chisqOverlaps));
for( j = 0; j < params->numChiSquareBins; j++)
{
......
......@@ -319,9 +319,18 @@ int main( int argc, char **argv )
verbose( "Created segments for null stream at %ld \n",
timeval_subtract(&startTime) );
PTFM[ifoNumber] = XLALCreateREAL8ArrayL( 2, 5, 5 );
PTFN[ifoNumber] = XLALCreateREAL8ArrayL( 2, 5, 5 );
PTFqVec[ifoNumber] = XLALCreateCOMPLEX8VectorSequence ( 5, numPoints );
if (params->approximant == FindChirpPTF)
{
PTFM[ifoNumber] = XLALCreateREAL8ArrayL( 2, 5, 5 );
PTFN[ifoNumber] = XLALCreateREAL8ArrayL( 2, 5, 5 );
PTFqVec[ifoNumber] = XLALCreateCOMPLEX8VectorSequence ( 5, numPoints );
}
else
{
PTFM[ifoNumber] = XLALCreateREAL8ArrayL( 2, 1, 1 );
PTFN[ifoNumber] = XLALCreateREAL8ArrayL( 2, 1, 1 );
PTFqVec[ifoNumber] = XLALCreateCOMPLEX8VectorSequence ( 1, numPoints );
}
}
/*------------------------------------------------------------------------*
......@@ -347,26 +356,56 @@ int main( int argc, char **argv )
fcInitParams = LALCalloc( 1, sizeof( *fcInitParams ) );
fcTmplt = LALCalloc( 1, sizeof( *fcTmplt ) );
fcTmpltParams = LALCalloc( 1, sizeof( *fcTmpltParams ) );
fcTmpltParams->approximant = FindChirpPTF;
fcTmpltParams->approximant = params->approximant;
fcTmpltParams->order = params->order;
/* Note that although non-spinning only uses Q1, the PTF
generator still generates Q1-5, thus size of these vectors */
fcTmplt->PTFQtilde =
XLALCreateCOMPLEX8VectorSequence( 5, numPoints / 2 + 1 );
fcTmpltParams->PTFQ = XLALCreateVectorSequence( 5, numPoints );
fcTmpltParams->PTFphi = XLALCreateVector( numPoints );
fcTmpltParams->PTFomega_2_3 = XLALCreateVector( numPoints );
fcTmpltParams->PTFe1 = XLALCreateVectorSequence( 3, numPoints );
fcTmpltParams->PTFe2 = XLALCreateVectorSequence( 3, numPoints );
if (params->approximant == FindChirpPTF)
{
fcTmplt->PTFQtilde =
XLALCreateCOMPLEX8VectorSequence( 5, numPoints / 2 + 1 );
fcTmpltParams->PTFQ = XLALCreateVectorSequence( 5, numPoints );
fcTmpltParams->PTFphi = XLALCreateVector( numPoints );
fcTmpltParams->PTFomega_2_3 = XLALCreateVector( numPoints );
fcTmpltParams->PTFe1 = XLALCreateVectorSequence( 3, numPoints );
fcTmpltParams->PTFe2 = XLALCreateVectorSequence( 3, numPoints );
}
else
{
fcTmplt->PTFQtilde =
XLALCreateCOMPLEX8VectorSequence( 1, numPoints / 2 + 1 );
fcTmplt->data = XLALCreateCOMPLEX8Vector( numPoints / 2 + 1 );
fcTmpltParams->xfacVec = XLALCreateVector(numPoints / 2 + 1 );
/* Set the values of xfacVec This is k^(-1/3) */
const REAL4 xfacExponent = -1.0/3.0;
REAL4 *xfac = NULL;
xfac = fcTmpltParams->xfacVec->data;
xfac[0] = 0;
for (ui = 1; ui < fcTmpltParams->xfacVec->length; ++ui)
xfac[ui] = pow( (REAL4) ui, xfacExponent );
}
fcTmpltParams->fwdPlan = XLALCreateForwardREAL4FFTPlan( numPoints, 0 );
fcTmpltParams->deltaT = 1.0/params->sampleRate;
fcTmpltParams->fLow = params->lowTemplateFrequency;
for( ifoNumber = 0; ifoNumber < LAL_NUM_IFO; ifoNumber++ )
{
if ( params->haveTrig[ifoNumber] )
{
PTFM[ifoNumber] = XLALCreateREAL8ArrayL( 2, 5, 5 );
PTFN[ifoNumber] = XLALCreateREAL8ArrayL( 2, 5, 5 );
PTFqVec[ifoNumber] = XLALCreateCOMPLEX8VectorSequence ( 5, numPoints );
if ( params->approximant == FindChirpPTF )
{
PTFM[ifoNumber] = XLALCreateREAL8ArrayL( 2, 5, 5 );
PTFN[ifoNumber] = XLALCreateREAL8ArrayL( 2, 5, 5 );
PTFqVec[ifoNumber] = XLALCreateCOMPLEX8VectorSequence ( 5, numPoints );
}
else
{
PTFM[ifoNumber] = XLALCreateREAL8ArrayL( 2, 1, 1 );
PTFN[ifoNumber] = XLALCreateREAL8ArrayL( 2, 1, 1 );
PTFqVec[ifoNumber] = XLALCreateCOMPLEX8VectorSequence ( 1, numPoints );
}
}
}
......@@ -429,9 +468,10 @@ int main( int argc, char **argv )
struct bankTemplateOverlaps *bankNormOverlaps = NULL;
struct bankComplexTemplateOverlaps *bankOverlaps = NULL;
struct bankDataOverlaps *dataOverlaps = NULL;
if ( params->doBankVeto )
{
verbose( "Initializing bank veto filters at %ld \n", timeval_subtract(&startTime));
/* Reads in and initializes the bank veto sub bank */
subBankSize = coh_PTF_read_sub_bank(params,&PTFBankTemplates);
bankNormOverlaps = LALCalloc( subBankSize,sizeof( *bankNormOverlaps));
......@@ -454,7 +494,10 @@ int main( int argc, char **argv )
/* Only store Q1. Structures used in fcTmpltParams will be overwritten */
for ( uj = 0 ; uj < (numPoints/2 +1) ; uj++ )
{
bankFcTmplts[ui].PTFQtilde->data[uj] = fcTmplt->PTFQtilde->data[uj];
if (params->approximant == FindChirpPTF)
bankFcTmplts[ui].PTFQtilde->data[uj] = fcTmplt->PTFQtilde->data[uj];
else
bankFcTmplts[ui].PTFQtilde->data[uj] = fcTmplt->data->data[uj];
}
}
/* Calculate the overlap between templates for bank veto */
......@@ -487,11 +530,13 @@ int main( int argc, char **argv )
*------------------------------------------------------------------------*/
/* Create the structures needed for the auto veto, if necessary */
UINT4 timeStepPoints = 0;
struct bankComplexTemplateOverlaps *autoTempOverlaps = NULL;
if ( params->doAutoVeto )
{
verbose( "Initializing auto veto filters at %ld \n", timeval_subtract(&startTime));
/* Initializations */
autoTempOverlaps = LALCalloc( params->numAutoPoints,
sizeof( *autoTempOverlaps));
......@@ -512,6 +557,7 @@ int main( int argc, char **argv )