diff --git a/README.md b/README.md
index 3e5b51a6841eb8e248159f9d6c4ce88c52a560d8..7f6b1d1d78c7cddc37883cb918aa6a21401b5f80 100644
--- a/README.md
+++ b/README.md
@@ -10,7 +10,8 @@ All examples can be run from their source scripts in [examples](examples), or
for each example there is descriptive documentation:
* [Making fake data with and without glitches](docs/make_fake_data.md)
-* [Fully coherent searches MCMC](docs/fully_coherent_search.md)
+* [Fully coherent MCMC search](docs/fully_coherent_search.md)
+* [Fully coherent MCMC search on data containing glitching signals](docs/fully_coherent_search_on_glitching_data.md)
## Installation
diff --git a/docs/fully_coherent_search_on_glitching_data.md b/docs/fully_coherent_search_on_glitching_data.md
new file mode 100644
index 0000000000000000000000000000000000000000..8e5ca97f42c9437a5198fa18969455307137c09b
--- /dev/null
+++ b/docs/fully_coherent_search_on_glitching_data.md
@@ -0,0 +1,79 @@
+# Fully coherent search on glitching data using MCMC
+
+This example applies the basic [fully coherent
+search](fully_coherent_search.md), to the glitching signal data set created in
+[make fake data](make_fake_data.md]). The aim here is to illustrate the effect
+such a signal can have on a fully-coherent search. The complete script for this
+example canbe found
+[here](../example/fully_cohrent_search_on_glitching_data.py).
+
+
+We use the same prior as in the basic fully-coherent search, except that we
+will modify the prior on `F0` to a flat uniform prior. The reason for this is
+to highlight the multimodal nature of the posterior which results from the
+glitch (a normal prior centered on one of the modes will bias one mode over
+the other). So our initial set up is
+
+```
+from pyfstat import MCMCSearch
+
+F0 = 30.0
+F1 = -1e-10
+F2 = 0
+Alpha = 5e-3
+Delta = 6e-2
+tref = 362750407.0
+
+tstart = 1000000000
+duration = 100*86400
+tend = tstart = duration
+
+theta_prior = {'F0': {'type': 'unif', 'lower': F0-5e-5,
+ 'upper': F0+5e-5},
+ 'F1': {'type': 'norm', 'loc': F1, 'scale': abs(1e-6*F1)},
+ 'F2': F2,
+ 'Alpha': Alpha,
+ 'Delta': Delta
+ }
+```
+
+Next, we will use 10 temperatures and a larger number of walkers - these have
+been tuned to illustrate the bimodal nature of the posterior
+
+```
+ntemps = 10
+betas = np.logspace(0, -30, ntemps)
+nwalkers = 500
+nsteps = [100, 100, 100]
+
+mcmc = MCMCSearch('fully_coherent_on_glitching_data', 'data',
+ sftlabel='glitch', sftdir='data',
+ theta_prior=theta_prior, tref=tref, tstart=tstart, tend=tend,
+ nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps, betas=betas,
+ scatter_val=1e-6)
+
+mcmc.run()
+mcmc.plot_corner(add_prior=True)
+```
+
+Running this takes slightly longer than the basic example (a few minutes) and
+produces a multimodal posterior:
+
+
+Clearly one central peak pertains to the original frequency `F0=30`. At
+`30+0.4e-5`, we find the second largest peak - this is the mode corresponding
+to the second half of the data. In reality, we would expect both peaks to be
+of equal size since the glitch occurs exactly half way through (see [how the
+data was made](make_fake_data.md)). We will confirm this later on by performing
+a grid search.
+
+Finally, the maximum twoF value found is
+
+```
+>>> mcmc.print_summary()
+Max twoF: 411.595
+```
+That is, compared to the basic search (on a smooth signal) which had a twoF of
+`1756.44177246` (in agreement with the predicted twoF), we have lost a large
+fraction of the SNR due to the glitch.
+
diff --git a/docs/img/fully_coherent_on_glitching_data_corner.png b/docs/img/fully_coherent_on_glitching_data_corner.png
new file mode 100644
index 0000000000000000000000000000000000000000..bfb47e5b294063edc8b3576ca7115495a053831d
Binary files /dev/null and b/docs/img/fully_coherent_on_glitching_data_corner.png differ
diff --git a/docs/make_fake_data.md b/docs/make_fake_data.md
index df730ee2f02910ed84c2c4d3e68ebeba6e819265..0869431aeffa2182457384c6bee3e751edc7f75c 100644
--- a/docs/make_fake_data.md
+++ b/docs/make_fake_data.md
@@ -74,6 +74,7 @@ data.run_makefakedata()
In fact, the previous two commands are wrapped together by a single call to
`data.make_data()` which we will use from now on.
+
## Glitching signal
We now want to generate a set of data which contains a *glitching signal*. We
@@ -83,8 +84,8 @@ another `Writer` instance called `glitch_data`, and then run `make_data()`
```
dtglitch = duration/2.0
-delta_F0 = 1e-6 * F0
-delta_F1 = 1e-5 * F1
+delta_F0 = 0.4e-5
+delta_F1 = 0
glitch_data = Writer(
label='glitch', outdir='data', tref=tref, tstart=tstart, F0=F0, F1=F1,
@@ -117,9 +118,9 @@ Delta = 5.999999999999999778e-02
h0 = 9.999999999999999604e-24
cosi = 0.000000000000000000e+00
psi = 0.000000000000000000e+00
-phi0 = -1.222261350197196007e+05
-Freq = 3.000003064156959098e+01
-f1dot = -1.000009999999999993e-10
+phi0 = -1.612440256772935390e+04
+Freq = 3.000000400000000056e+01
+f1dot = -1.000000000000000036e-10
f2dot = 0.000000000000000000e+00
refTime = 362750407.000000
transientWindowType=rect
@@ -130,4 +131,14 @@ transientTauDays=50.000
The glitch config file uses transient windows to create two non-overlapping,
but continuous signals.
+## Expected twoF
+
+Finally, the `Writer` class also provides a wrapper of `lalapps_PredictFstat`.
+So calling
+
+```
+>>> print data.predict_fstat()
+1721.1
+```
+Notice that the predicted value will be the same for both sets of data.
diff --git a/examples/fully_coherent_search_on_glitching_data.py b/examples/fully_coherent_search_on_glitching_data.py
new file mode 100644
index 0000000000000000000000000000000000000000..2ea029b0af205bd5ff8d7a5f0c35ade1eee9687a
--- /dev/null
+++ b/examples/fully_coherent_search_on_glitching_data.py
@@ -0,0 +1,35 @@
+from pyfstat import MCMCSearch
+import numpy as np
+
+F0 = 30.0
+F1 = -1e-10
+F2 = 0
+Alpha = 5e-3
+Delta = 6e-2
+tref = 362750407.0
+
+tstart = 1000000000
+duration = 100*86400
+tend = tstart = duration
+
+theta_prior = {'F0': {'type': 'unif', 'lower': F0-5e-5,
+ 'upper': F0+5e-5},
+ 'F1': {'type': 'norm', 'loc': F1, 'scale': abs(1e-6*F1)},
+ 'F2': F2,
+ 'Alpha': Alpha,
+ 'Delta': Delta
+ }
+
+ntemps = 10
+betas = np.logspace(0, -30, ntemps)
+nwalkers = 500
+nsteps = [100, 100, 100]
+
+mcmc = MCMCSearch('fully_coherent_on_glitching_data', 'data',
+ sftlabel='glitch', sftdir='data',
+ theta_prior=theta_prior, tref=tref, tstart=tstart, tend=tend,
+ nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps, betas=betas,
+ scatter_val=1e-6)
+mcmc.run()
+mcmc.plot_corner(add_prior=True)
+mcmc.print_summary()
diff --git a/examples/make_fake_data.py b/examples/make_fake_data.py
index 6f74c3e9a68cb5f6bf1e8956a85638e4e656d56e..b6315e1140817fc92e581fcbef39e63572dd2e51 100644
--- a/examples/make_fake_data.py
+++ b/examples/make_fake_data.py
@@ -27,8 +27,8 @@ data.make_data()
# Next, taking the same signal parameters, we include a glitch half way through
dtglitch = duration/2.0
-delta_F0 = 1e-6 * F0
-delta_F1 = 1e-5 * F1
+delta_F0 = 0.4e-5
+delta_F1 = 0
glitch_data = Writer(
label='glitch', outdir='data', tref=tref, tstart=tstart, F0=F0, F1=F1,
@@ -36,3 +36,7 @@ glitch_data = Writer(
dtglitch=dtglitch, delta_F0=delta_F0, delta_F1=delta_F1)
glitch_data.make_data()
+
+# The predicted twoF, given by lalapps_predictFstat can be accsessed by
+
+print data.predict_fstat()