diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 72897ad8ca2b633aace5fc1510d3635161d6caf5..1acb1355e46274468a9022f6b891f5847f0a8b0d 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -42,6 +42,7 @@ python-3:
- coverage --version
- coverage erase
- coverage run --debug=trace --include=/opt/conda/lib/python3.6/site-packages/tupak* -a test/conversion_tests.py
+ - coverage run --include=/opt/conda/lib/python3.6/site-packages/tupak* -a test/calibration_tests.py
- coverage run --include=/opt/conda/lib/python3.6/site-packages/tupak* -a test/detector_tests.py
- coverage run --include=/opt/conda/lib/python3.6/site-packages/tupak* -a test/utils_tests.py
- coverage run --include=/opt/conda/lib/python3.6/site-packages/tupak* -a test/prior_tests.py
diff --git a/examples/injection_examples/calibration_example.py b/examples/injection_examples/calibration_example.py
new file mode 100644
index 0000000000000000000000000000000000000000..d2cf310554365a4e66e7d70bae285837d3531508
--- /dev/null
+++ b/examples/injection_examples/calibration_example.py
@@ -0,0 +1,86 @@
+#!/bin/python
+"""
+Tutorial to demonstrate running parameter estimation with calibration
+uncertainties included.
+"""
+from __future__ import division, print_function
+
+import numpy as np
+import tupak
+
+# Set the duration and sampling frequency of the data segment
+# that we're going to create and inject the signal into.
+
+duration = 4.
+sampling_frequency = 2048.
+
+# Specify the output directory and the name of the simulation.
+outdir = 'outdir'
+label = 'calibration'
+tupak.core.utils.setup_logger(outdir=outdir, label=label)
+
+# Set up a random seed for result reproducibility. This is optional!
+np.random.seed(88170235)
+
+# We are going to inject a binary black hole waveform. We first establish a
+# dictionary of parameters that includes all of the different waveform
+# parameters, including masses of the two black holes (mass_1, mass_2),
+# spins of both black holes (a, tilt, phi), etc.
+injection_parameters = dict(
+ mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
+ phi_12=1.7, phi_jl=0.3, luminosity_distance=2000., iota=0.4, psi=2.659,
+ phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
+
+# Fixed arguments passed into the source model
+waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
+ reference_frequency=50.)
+
+# Create the waveform_generator using a LAL BinaryBlackHole source function
+waveform_generator = tupak.gw.WaveformGenerator(
+ duration=duration, sampling_frequency=sampling_frequency,
+ frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole,
+ parameters=injection_parameters,waveform_arguments=waveform_arguments)
+
+# Set up interferometers. In this case we'll use three interferometers
+# (LIGO-Hanford (H1), LIGO-Livingston (L1), and Virgo (V1)).
+# These default to their design sensitivity
+ifos = tupak.gw.detector.InterferometerList(['H1', 'L1', 'V1'])
+for ifo in ifos:
+ injection_parameters.update({
+ 'recalib_{}_amplitude_{}'.format(ifo.name, ii): 0.1 for ii in range(5)})
+ injection_parameters.update({
+ 'recalib_{}_phase_{}'.format(ifo.name, ii): 0.01 for ii in range(5)})
+ ifo.calibration_model = tupak.gw.calibration.CubicSpline(
+ prefix='recalib_{}_'.format(ifo.name),
+ minimum_frequency=ifo.minimum_frequency,
+ maximum_frequency=ifo.maximum_frequency, n_points=5)
+ifos.set_strain_data_from_power_spectral_densities(
+ sampling_frequency=sampling_frequency, duration=duration)
+ifos.inject_signal(parameters=injection_parameters,
+ waveform_generator=waveform_generator)
+
+# Set up prior, which is a dictionary
+# Here we fix the injected cbc parameters and most of the calibration parameters
+# to the injected values.
+# We allow a subset of the calibration parameters to vary.
+priors = injection_parameters.copy()
+for key in injection_parameters:
+ if 'recalib' in key:
+ priors[key] = injection_parameters[key]
+for name in ['recalib_H1_amplitude_0', 'recalib_H1_amplitude_1']:
+ priors[name] = tupak.prior.Gaussian(
+ mu=0, sigma=0.2, name=name, latex_label='H1 $A_{}$'.format(name[-1]))
+
+# Initialise the likelihood by passing in the interferometer data (IFOs) and
+# the waveform generator
+likelihood = tupak.gw.GravitationalWaveTransient(
+ interferometers=ifos, waveform_generator=waveform_generator)
+
+# Run sampler. In this case we're going to use the `dynesty` sampler
+result = tupak.run_sampler(
+ likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000,
+ injection_parameters=injection_parameters, outdir=outdir, label=label)
+
+# make some plots of the outputs
+result.plot_corner()
+
diff --git a/examples/injection_examples/marginalized_likelihood.py b/examples/injection_examples/marginalized_likelihood.py
index 89fb01ea7a383019b70bc4ce5962c14d7b76f3ba..08b9fdf2bea14d2cb5102c1158ebc0d54f00cc14 100644
--- a/examples/injection_examples/marginalized_likelihood.py
+++ b/examples/injection_examples/marginalized_likelihood.py
@@ -36,7 +36,7 @@ IFOs = [tupak.gw.detector.get_interferometer_with_fake_noise_and_injection(
# Set up prior
priors = tupak.gw.prior.BBHPriorSet()
# These parameters will not be sampled
-for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'phase', 'iota', 'ra', 'dec', 'geocent_time']:
+for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'iota', 'ra', 'dec', 'geocent_time']:
priors[key] = injection_parameters[key]
# Initialise GravitationalWaveTransient
diff --git a/examples/open_data_examples/GW150914.py b/examples/open_data_examples/GW150914.py
index 06fa0b810728a246e2751fbd9351d0cadd4f67fc..d619dd7e1928a7dadb95a6cdf2ba5eca93aaa16b 100644
--- a/examples/open_data_examples/GW150914.py
+++ b/examples/open_data_examples/GW150914.py
@@ -36,13 +36,15 @@ prior = tupak.gw.prior.BBHPriorSet(filename='GW150914.prior')
# creates the frequency-domain strain. In this instance, we are using the
# `lal_binary_black_hole model` source model. We also pass other parameters:
# the waveform approximant and reference frequency.
-waveform_generator = tupak.gw.WaveformGenerator(frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole,
- waveform_arguments={'waveform_approximant': 'IMRPhenomPv2',
- 'reference_frequency': 50})
+waveform_generator = tupak.WaveformGenerator(
+ frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole,
+ waveform_arguments={'waveform_approximant': 'IMRPhenomPv2',
+ 'reference_frequency': 50})
# In this step, we define the likelihood. Here we use the standard likelihood
# function, passing it the data and the waveform generator.
-likelihood = tupak.gw.likelihood.GravitationalWaveTransient(interferometers, waveform_generator)
+likelihood = tupak.gw.likelihood.GravitationalWaveTransient(
+ interferometers, waveform_generator)
# Finally, we run the sampler. This function takes the likelihood and prior
# along with some options for how to do the sampling and how to save the data
diff --git a/test/calibration_tests.py b/test/calibration_tests.py
new file mode 100644
index 0000000000000000000000000000000000000000..b481e2936ed88ec6d4cf68136c495f0ef4d25539
--- /dev/null
+++ b/test/calibration_tests.py
@@ -0,0 +1,53 @@
+from tupak.gw import calibration
+import unittest
+import numpy as np
+
+
+class TestBaseClass(unittest.TestCase):
+
+ def setUp(self):
+ self.model = calibration.Recalibrate()
+
+ def tearDown(self):
+ del self.model
+
+ def test_calibration_factor(self):
+ frequency_array = np.linspace(20, 1024, 1000)
+ cal_factor = self.model.get_calibration_factor(frequency_array)
+ assert np.alltrue(cal_factor.real == np.ones_like(frequency_array))
+
+
+class TestCubicSpline(unittest.TestCase):
+
+ def setUp(self):
+ self.model = calibration.CubicSpline(
+ prefix='recalib_', minimum_frequency=20, maximum_frequency=1024,
+ n_points=5)
+ self.parameters = {'recalib_{}_{}'.format(param, ii): 0.0
+ for ii in range(5)
+ for param in ['amplitude', 'phase']}
+
+ def tearDown(self):
+ del self.model
+ del self.parameters
+
+ def test_calibration_factor(self):
+ frequency_array = np.linspace(20, 1024, 1000)
+ cal_factor = self.model.get_calibration_factor(frequency_array,
+ **self.parameters)
+ assert np.alltrue(cal_factor.real == np.ones_like(frequency_array))
+
+
+class TestCubicSplineRequiresFourNodes(unittest.TestCase):
+
+ def test_cannot_instantiate_with_too_few_nodes(self):
+ for ii in range(6):
+ if ii < 4:
+ with self.assertRaises(ValueError):
+ calibration.CubicSpline('test', 1, 10, ii)
+ else:
+ calibration.CubicSpline('test', 1, 10, ii)
+
+
+if __name__ == '__main__':
+ unittest.main()
diff --git a/tupak/core/prior.py b/tupak/core/prior.py
index a589c438f7e7c3a172e10607037f0c84874c410e..8ea049c6befa343e9c0845d124f12f699e9c2f1f 100644
--- a/tupak/core/prior.py
+++ b/tupak/core/prior.py
@@ -26,7 +26,11 @@ class PriorSet(dict):
dict.__init__(self)
if type(dictionary) is dict:
self.update(dictionary)
- elif filename:
+ elif type(dictionary) is str:
+ logger.debug('Argument "dictionary" is a string.'
+ + ' Assuming it is intended as a file name.')
+ self.read_in_file(dictionary)
+ elif type(filename) is str:
self.read_in_file(filename)
def write_to_file(self, outdir, label):
@@ -41,7 +45,7 @@ class PriorSet(dict):
"""
utils.check_directory_exists_and_if_not_mkdir(outdir)
- prior_file = os.path.join(outdir, "{}_prior.txt".format(label))
+ prior_file = os.path.join(outdir, "{}.prior".format(label))
logger.debug("Writing priors to {}".format(prior_file))
with open(prior_file, "w") as outfile:
for key in self.keys():
diff --git a/tupak/gw/__init__.py b/tupak/gw/__init__.py
index 9969e9d5f35434cc0247c1470c00bde6ee7f45ed..47063bf1e2790e5ef41ad096f9723905994e9cd9 100644
--- a/tupak/gw/__init__.py
+++ b/tupak/gw/__init__.py
@@ -6,6 +6,8 @@ import tupak.gw.prior
import tupak.gw.source
import tupak.gw.utils
import tupak.gw.waveform_generator
+from . import calibration
+
from tupak.gw.waveform_generator import WaveformGenerator
from tupak.gw.likelihood import GravitationalWaveTransient
diff --git a/tupak/gw/calibration.py b/tupak/gw/calibration.py
new file mode 100644
index 0000000000000000000000000000000000000000..e54757a86a9c60bfa029271e176b75a23c714112
--- /dev/null
+++ b/tupak/gw/calibration.py
@@ -0,0 +1,110 @@
+""" Functions for adding calibration factors to waveform templates.
+"""
+
+import numpy as np
+from scipy.interpolate import UnivariateSpline
+
+
+class Recalibrate(object):
+
+ name = 'none'
+
+ def __init__(self, prefix='recalib_'):
+ """
+ Base calibration object. This applies no transformation
+
+ Parameters
+ ----------
+ prefix: str
+ Prefix on parameters relating to the calibration.
+ """
+ self.params = dict()
+ self.prefix = prefix
+
+ def get_calibration_factor(self, frequency_array, **params):
+ """Apply calibration model
+
+ This method should be overwritten by subclasses
+
+ Parameters
+ ----------
+ frequency_array: array-like
+ The frequency values to calculate the calibration factor for.
+ params : dict
+ Dictionary of sampling parameters which includes
+ calibration parameters.
+
+ Returns
+ -------
+ calibration_factor : array-like
+ The factor to multiply the strain by.
+ """
+ self.set_calibration_parameters(**params)
+ return np.ones_like(frequency_array)
+
+ def set_calibration_parameters(self, **params):
+ self.params.update({key[len(self.prefix):]: params[key] for key in params
+ if self.prefix in key})
+
+
+class CubicSpline(Recalibrate):
+
+ name = 'cubic_spline'
+
+ def __init__(self, prefix, minimum_frequency, maximum_frequency, n_points):
+ """
+ Cubic spline recalibration
+
+ see https://dcc.ligo.org/DocDB/0116/T1400682/001/calnote.pdf
+
+ This assumes the spline points follow
+ np.logspace(np.log(minimum_frequency), np.log(maximum_frequency), n_points)
+
+ Parameters
+ ----------
+ prefix: str
+ Prefix on parameters relating to the calibration.
+ minimum_frequency: float
+ minimum frequency of spline points
+ maximum_frequency: float
+ maximum frequency of spline points
+ n_points: int
+ number of spline points
+ """
+ Recalibrate.__init__(self, prefix=prefix)
+ if n_points < 4:
+ raise ValueError('Cubic spline calibration requires at least 4 spline nodes.')
+ self.n_points = n_points
+ self.spline_points = np.logspace(np.log10(minimum_frequency), np.log10(maximum_frequency), n_points)
+
+ def get_calibration_factor(self, frequency_array, **params):
+ """Apply calibration model
+
+ Parameters
+ ----------
+ frequency_array: array-like
+ The frequency values to calculate the calibration factor for.
+ prefix: str
+ Prefix for calibration parameter names
+ params : dict
+ Dictionary of sampling parameters which includes
+ calibration parameters.
+
+ Returns
+ -------
+ calibration_factor : array-like
+ The factor to multiply the strain by.
+ """
+ self.set_calibration_parameters(**params)
+ amplitude_parameters = [self.params['amplitude_{}'.format(ii)] for ii in range(self.n_points)]
+ amplitude_spline = UnivariateSpline(self.spline_points, amplitude_parameters)
+ delta_amplitude = amplitude_spline(frequency_array)
+
+ phase_parameters = [self.params['phase_{}'.format(ii)] for ii in range(self.n_points)]
+ phase_spline = UnivariateSpline(self.spline_points, phase_parameters)
+ delta_phase = phase_spline(frequency_array)
+
+ calibration_factor = (1 + delta_amplitude) * (2 + 1j * delta_phase) / (2 - 1j * delta_phase)
+
+ return calibration_factor
+
diff --git a/tupak/gw/detector.py b/tupak/gw/detector.py
index 842fc37333a9b1a1f538f23f8cd90e1de0736b84..f20c188ab0e70cd5580d7a76e12cc72ff91ab200 100644
--- a/tupak/gw/detector.py
+++ b/tupak/gw/detector.py
@@ -4,12 +4,13 @@ import os
import matplotlib.pyplot as plt
import numpy as np
-import scipy
+from scipy.signal.windows import tukey
from scipy.interpolate import interp1d
import tupak.gw.utils
from tupak.core import utils
from tupak.core.utils import logger
+from .calibration import Recalibrate
try:
import gwpy
@@ -85,7 +86,7 @@ class InterferometerList(list):
`waveform_generator.frequency_domain_strain()`. If
`waveform_generator` is also given, the injection_polarizations will
be calculated directly and this argument can be ignored.
- waveform_generator: tupak.gw.waveform_generator
+ waveform_generator: tupak.gw.waveform_generator.WaveformGenerator
A WaveformGenerator instance using the source model to inject. If
`injection_polarizations` is given, this will be ignored.
@@ -353,7 +354,7 @@ class InterferometerStrainData(object):
self.roll_off = roll_off
elif alpha is not None:
self.roll_off = alpha * self.duration / 2
- window = scipy.signal.windows.tukey(len(self._time_domain_strain), alpha=self.alpha)
+ window = tukey(len(self._time_domain_strain), alpha=self.alpha)
self.window_factor = np.mean(window**2)
return window
@@ -741,9 +742,9 @@ class InterferometerStrainData(object):
class Interferometer(object):
"""Class for the Interferometer """
- def __init__(self, name, power_spectral_density, minimum_frequency,
- maximum_frequency, length, latitude, longitude, elevation,
- xarm_azimuth, yarm_azimuth, xarm_tilt=0., yarm_tilt=0.):
+ def __init__(self, name, power_spectral_density, minimum_frequency, maximum_frequency,
+ length, latitude, longitude, elevation, xarm_azimuth, yarm_azimuth,
+ xarm_tilt=0., yarm_tilt=0., calibration_model=Recalibrate()):
"""
Instantiate an Interferometer object.
@@ -774,6 +775,9 @@ class Interferometer(object):
ellipsoid earth model in LIGO-T980044-08.
yarm_tilt: float, optional
Tilt of the y arm in radians above the horizontal.
+ calibration_model: Recalibration
+ Calibration model, this applies the calibration correction to the
+ template, the default model applies no correction.
"""
self.__x_updated = False
self.__y_updated = False
@@ -791,6 +795,7 @@ class Interferometer(object):
self.yarm_tilt = yarm_tilt
self.power_spectral_density = power_spectral_density
self.time_marginalization = False
+ self.calibration_model = calibration_model
self._strain_data = InterferometerStrainData(
minimum_frequency=minimum_frequency,
maximum_frequency=maximum_frequency)
@@ -1187,6 +1192,9 @@ class Interferometer(object):
signal_ifo = signal_ifo * np.exp(
-1j * 2 * np.pi * dt * self.frequency_array)
+ signal_ifo *= self.calibration_model.get_calibration_factor(
+ self.frequency_array, prefix='recalib_{}_'.format(self.name), **parameters)
+
return signal_ifo
def inject_signal(self, parameters=None, injection_polarizations=None,
diff --git a/tupak/gw/prior.py b/tupak/gw/prior.py
index 3346a155762bd8a8d7edb8a786c9ae00917b3f68..3df269375b05c38ee03197c98a39792182556437 100644
--- a/tupak/gw/prior.py
+++ b/tupak/gw/prior.py
@@ -1,7 +1,8 @@
import os
-from tupak.core.prior import PriorSet, FromFile, Prior
-
-from tupak.core.utils import logger
+from ..core.prior import PriorSet, FromFile, Prior, DeltaFunction, Gaussian
+from ..core.utils import logger
+import numpy as np
+from scipy.interpolate import UnivariateSpline
class UniformComovingVolume(FromFile):
@@ -114,3 +115,168 @@ Prior._default_latex_labels = {
'psi': '$\psi$',
'phase': '$\phi$',
'geocent_time': '$t_c$'}
+
+
+class CalibrationPriorSet(PriorSet):
+
+ def __init__(self, dictionary=None, filename=None):
+ """ Initialises a Prior set for Binary Black holes
+
+ Parameters
+ ----------
+ dictionary: dict, optional
+ See superclass
+ filename: str, optional
+ See superclass
+ """
+ if dictionary is None and filename is not None:
+ filename = os.path.join(os.path.dirname(__file__),
+ 'prior_files', filename)
+ PriorSet.__init__(self, dictionary=dictionary, filename=filename)
+ self.source = None
+
+ def write_to_file(self, outdir, label):
+ """
+ Write the prior to file. This includes information about the source if
+ possible.
+
+ Parameters
+ ----------
+ outdir: str
+ Output directory.
+ label: str
+ Label for prior.
+ """
+ PriorSet.write_to_file(self, outdir=outdir, label=label)
+ if self.source is not None:
+ prior_file = os.path.join(outdir, "{}.prior".format(label))
+ with open(prior_file, "a") as outfile:
+ outfile.write("# prior source file is {}".format(self.source))
+
+ @staticmethod
+ def from_envelope_file(envelope_file, minimum_frequency,
+ maximum_frequency, n_nodes, label):
+ """
+ Load in the calibration envelope.
+
+ This is a text file with columns:
+ frequency median-amplitude median-phase -1-sigma-amplitude
+ -1-sigma-phase +1-sigma-amplitude +1-sigma-phase
+
+ Parameters
+ ----------
+ envelope_file: str
+ Name of file to read in.
+ minimum_frequency: float
+ Minimum frequency for the spline.
+ maximum_frequency: float
+ Minimum frequency for the spline.
+ n_nodes: int
+ Number of nodes for the spline.
+ label: str
+ Label for the names of the parameters, e.g., recalib_H1_
+
+ Returns
+ -------
+ prior: PriorSet
+ Priors for the relevant parameters.
+ This includes the frequencies of the nodes which are _not_ sampled.
+ """
+ calibration_data = np.genfromtxt(envelope_file).T
+ frequency_array = calibration_data[0]
+ amplitude_median = 1 - calibration_data[1]
+ phase_median = calibration_data[2]
+ amplitude_sigma = (calibration_data[4] - calibration_data[2]) / 2
+ phase_sigma = (calibration_data[5] - calibration_data[3]) / 2
+
+ nodes = np.logspace(np.log10(minimum_frequency),
+ np.log10(maximum_frequency), n_nodes)
+
+ amplitude_mean_nodes =\
+ UnivariateSpline(frequency_array, amplitude_median)(nodes)
+ amplitude_sigma_nodes =\
+ UnivariateSpline(frequency_array, amplitude_sigma)(nodes)
+ phase_mean_nodes =\
+ UnivariateSpline(frequency_array, phase_median)(nodes)
+ phase_sigma_nodes =\
+ UnivariateSpline(frequency_array, phase_sigma)(nodes)
+
+ prior = CalibrationPriorSet()
+ for ii in range(n_nodes):
+ name = "recalib_{}_amplitude_{}".format(label, ii)
+ latex_label = "$A^{}_{}$".format(label, ii)
+ prior[name] = Gaussian(mu=amplitude_mean_nodes[ii],
+ sigma=amplitude_sigma_nodes[ii],
+ name=name, latex_label=latex_label)
+ for ii in range(n_nodes):
+ name = "recalib_{}_phase_{}".format(label, ii)
+ latex_label = "$\\phi^{}_{}$".format(label, ii)
+ prior[name] = Gaussian(mu=phase_mean_nodes[ii],
+ sigma=phase_sigma_nodes[ii],
+ name=name, latex_label=latex_label)
+ for ii in range(n_nodes):
+ name = "recalib_{}_frequency_{}".format(label, ii)
+ latex_label = "$f^{}_{}$".format(label, ii)
+ prior[name] = DeltaFunction(peak=nodes[ii], name=name,
+ latex_label=latex_label)
+ prior.source = os.path.abspath(envelope_file)
+ return prior
+
+ @staticmethod
+ def constant_uncertainty_spline(
+ amplitude_sigma, phase_sigma, minimum_frequency, maximum_frequency,
+ n_nodes, label):
+ """
+ Make prior assuming constant in frequency calibration uncertainty.
+
+ This assumes Gaussian fluctuations about 0.
+
+ Parameters
+ ----------
+ amplitude_sigma: float
+ Uncertainty in the amplitude.
+ phase_sigma: float
+ Uncertainty in the phase.
+ minimum_frequency: float
+ Minimum frequency for the spline.
+ maximum_frequency: float
+ Minimum frequency for the spline.
+ n_nodes: int
+ Number of nodes for the spline.
+ label: str
+ Label for the names of the parameters, e.g., recalib_H1_
+
+ Returns
+ -------
+ prior: PriorSet
+ Priors for the relevant parameters.
+ This includes the frequencies of the nodes which are _not_ sampled.
+ """
+ nodes = np.logspace(np.log10(minimum_frequency),
+ np.log10(maximum_frequency), n_nodes)
+
+ amplitude_mean_nodes = [0] * n_nodes
+ amplitude_sigma_nodes = [amplitude_sigma] * n_nodes
+ phase_mean_nodes = [0] * n_nodes
+ phase_sigma_nodes = [phase_sigma] * n_nodes
+
+ prior = CalibrationPriorSet()
+ for ii in range(n_nodes):
+ name = "recalib_{}_amplitude_{}".format(label, ii)
+ latex_label = "$A^{}_{}$".format(label, ii)
+ prior[name] = Gaussian(mu=amplitude_mean_nodes[ii],
+ sigma=amplitude_sigma_nodes[ii],
+ name=name, latex_label=latex_label)
+ for ii in range(n_nodes):
+ name = "recalib_{}_phase_{}".format(label, ii)
+ latex_label = "$\\phi^{}_{}$".format(label, ii)
+ prior[name] = Gaussian(mu=phase_mean_nodes[ii],
+ sigma=phase_sigma_nodes[ii],
+ name=name, latex_label=latex_label)
+ for ii in range(n_nodes):
+ name = "recalib_{}_frequency_{}".format(label, ii)
+ latex_label = "$f^{}_{}$".format(label, ii)
+ prior[name] = DeltaFunction(peak=nodes[ii], name=name,
+ latex_label=latex_label)
+
+ return prior
diff --git a/tupak/gw/prior_files/GW150914.prior b/tupak/gw/prior_files/GW150914.prior
index 6a6c9390bafd93cf8e31e980e6766c10ad98fde3..276a6b9f710b3835206441279b2b3ffbe7644b32 100644
--- a/tupak/gw/prior_files/GW150914.prior
+++ b/tupak/gw/prior_files/GW150914.prior
@@ -14,3 +14,35 @@ iota = Sine(name='iota')
psi = Uniform(name='psi', minimum=0, maximum=2 * np.pi)
phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi)
geocent_time = Uniform(1126259462.322, 1126259462.522, name='geocent_time')
+# These are the calibration parameters as described in
+# https://journals.aps.org/prx/abstract/10.1103/PhysRevX.6.041015
+# recalib_H1_frequency_0 = 20
+# recalib_H1_frequency_1 = 54
+# recalib_H1_frequency_2 = 143
+# recalib_H1_frequency_3 = 383
+# recalib_H1_frequency_4 = 1024
+# recalib_H1_amplitude_0 = Gaussian(mu=0, sigma=0.048, name='recalib_H1_amplitude_0), '$\\delta A_{H0}$'
+# recalib_H1_amplitude_1 = Gaussian(mu=0, sigma=0.048, name='recalib_H1_amplitude_1), '$\\delta A_{H1}$'
+# recalib_H1_amplitude_2 = Gaussian(mu=0, sigma=0.048, name='recalib_H1_amplitude_2), '$\\delta A_{H2}$'
+# recalib_H1_amplitude_3 = Gaussian(mu=0, sigma=0.048, name='recalib_H1_amplitude_3), '$\\delta A_{H3}$'
+# recalib_H1_amplitude_4 = Gaussian(mu=0, sigma=0.048, name='recalib_H1_amplitude_4), '$\\delta A_{H4}$'
+# recalib_H1_phase_0 = Gaussian(mu=0, sigma=0.056, name='recalib_H1_phase_0', '$\\delta \\phi_{H0}$')
+# recalib_H1_phase_1 = Gaussian(mu=0, sigma=0.056, name='recalib_H1_phase_1', '$\\delta \\phi_{H1}$')
+# recalib_H1_phase_2 = Gaussian(mu=0, sigma=0.056, name='recalib_H1_phase_2', '$\\delta \\phi_{H2}$')
+# recalib_H1_phase_3 = Gaussian(mu=0, sigma=0.056, name='recalib_H1_phase_3', '$\\delta \\phi_{H3}$')
+# recalib_H1_phase_4 = Gaussian(mu=0, sigma=0.056, name='recalib_H1_phase_4', '$\\delta \\phi_{H4}$')
+# recalib_L1_frequency_0 = 20
+# recalib_L1_frequency_1 = 54
+# recalib_L1_frequency_2 = 143
+# recalib_L1_frequency_3 = 383
+# recalib_L1_frequency_4 = 1024
+# recalib_L1_amplitude_0 = Gaussian(mu=0, sigma=0.082, name='recalib_L1_amplitude_0), '$\\delta A_{L0}$'
+# recalib_L1_amplitude_1 = Gaussian(mu=0, sigma=0.082, name='recalib_L1_amplitude_1), '$\\delta A_{L1}$'
+# recalib_L1_amplitude_2 = Gaussian(mu=0, sigma=0.082, name='recalib_L1_amplitude_2), '$\\delta A_{L2}$'
+# recalib_L1_amplitude_3 = Gaussian(mu=0, sigma=0.082, name='recalib_L1_amplitude_3), '$\\delta A_{L3}$'
+# recalib_L1_amplitude_4 = Gaussian(mu=0, sigma=0.082, name='recalib_L1_amplitude_4), '$\\delta A_{L4}$'
+# recalib_L1_phase_0 = Gaussian(mu=0, sigma=0.073, name='recalib_L1_phase_0', '$\\delta \\phi_{L0}$')
+# recalib_L1_phase_1 = Gaussian(mu=0, sigma=0.073, name='recalib_L1_phase_1', '$\\delta \\phi_{L1}$')
+# recalib_L1_phase_2 = Gaussian(mu=0, sigma=0.073, name='recalib_L1_phase_2', '$\\delta \\phi_{L2}$')
+# recalib_L1_phase_3 = Gaussian(mu=0, sigma=0.073, name='recalib_L1_phase_3', '$\\delta \\phi_{L3}$')
+# recalib_L1_phase_4 = Gaussian(mu=0, sigma=0.073, name='recalib_L1_phase_4', '$\\delta \\phi_{L4}$')
diff --git a/tupak/gw/utils.py b/tupak/gw/utils.py
index 741053af92a70644730260b410587778aef4b70b..c146a204634423cf3daec234b38135c4be5fec6e 100644
--- a/tupak/gw/utils.py
+++ b/tupak/gw/utils.py
@@ -2,9 +2,9 @@ from __future__ import division
import os
import numpy as np
-from scipy import signal
-from tupak.core.utils import gps_time_to_gmst, ra_dec_to_theta_phi, speed_of_light, nfft, logger
+from ..core.utils import (gps_time_to_gmst, ra_dec_to_theta_phi, speed_of_light,
+ nfft, logger)
try:
from gwpy.signal import filter_design