knitr::opts_chunk$set(echo = TRUE) library(hadron)
First we need to create some artificial data in order to get a feeling whether this function actually works. We will take two states and add some noise.
time_extent <- 96 e0_val <- 0.12 e1_val <- 0.0 rel_amplitude <- 0.4 boot_R <- 99 n_meas <- 100
Next we generate bootstrap samples for the energies and amplitudes that enter here.
a0_val <- 1.0 a1_val <- rel_amplitude a0_err <- 0.001 a1_err <- 0.001 a0_orig <- rnorm(n_meas, a0_val, a0_err) a1_orig <- rnorm(n_meas, a1_val, a1_err) e0_err <- 0.001 e1_err <- 0.001 e0_orig <- rnorm(n_meas, e0_val, e0_err) e1_orig <- rnorm(n_meas, e1_val, e1_err)
The signal is just the sum of two cosh-terms.
t <- 0:(time_extent - 1) signal_val <- a0_val * (exp(-e0_val * t) + exp(-e0_val * (time_extent - t))) + a1_val * (exp(-e1_val * t) + exp(-e1_val * (time_extent - t))) signal_orig <- do.call(cbind, lapply(t, function (t) a0_orig * (exp(-e0_orig * t) + exp(-e0_orig * (time_extent - t))) + a1_orig * (exp(-e1_orig * t) + exp(-e1_orig * (time_extent - t))))) signal_orig <- signal_orig + matrix(rnorm(length(signal_orig), c(signal_orig), 0.1), ncol = ncol(signal_orig)) signal_err <- apply(signal_orig, 2, sd)
We take a look at our artificial data.
plotwitherror(t, signal_val, signal_err, log = 'y', main = 'Artificial pure signal')
All we have are the fake measurements, so we construct a cf
object from them and bootstrap that.
corr <- cf_orig(cf_meta(Time = time_extent), cf = signal_orig) corr_boot <- bootstrap.cf(symmetrise.cf(corr), boot.R = boot_R)
Now it looks like a correlator with very tiny errors.
plot(corr_boot, log = 'y', main = 'Symmetrized correlator with noise')
In the effective mass we can see that we are hopelessly far away from the actual mass that I want to get out.
effmass_solve <- bootstrap.effectivemass(corr_boot, type = 'solve') plot(effmass_solve, ylim = range(effmass_solve$t0, e0_val, na.rm = TRUE)) abline(h = e0_val)
Shifting the correlator gives us a sensible effective mass because it removes the constant term.
corr_shifted <- takeTimeDiff.cf(corr_boot) effmass_shifted <- bootstrap.effectivemass(corr_shifted, type = 'shifted') plot(effmass_shifted) abline(h = e0_val)
We first try the “single” model, that will of course fail miserably.
fit_single <- new_matrixfit( corr_boot, t1 = 1, t2 = 47, useCov = TRUE, model = 'single', fit.method = 'lm') plot(fit_single, log = 'y', main = 'Fit with “single” model') residual_plot(fit_single)
The “shifted” model will of course work.
fit_shifted <- new_matrixfit( corr_shifted, t1 = 1, t2 = 47, model = 'shifted', useCov = TRUE, fit.method = 'lm') plot(fit_shifted, log = 'y', main = 'Fit with “shifted” model') residual_plot(fit_shifted)
And with the “n_particle” model we can also get the result. It is not entirely correct because I must not use the original data for $E_1$ but rather the bootstrapped one. But I do not want to bootstrap that vector since I lack a function to do it nicely for me. So I just generate it uncorrelated.
fit_n_particle <- new_matrixfit( corr_boot, t1 = 1, t2 = 47, model = 'n_particles', fit.method = 'lm', useCov = TRUE, higher_states = list(val = c(e1_val), ampl = c(1), boot = matrix(rnorm(boot_R, e1_val, e1_err / sqrt(n_meas)), ncol = 1))) plot(fit_n_particle, log = 'y', main = 'Fit with “n_particle” model') residual_plot(fit_n_particle)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.