In etmc/hadron: Analysis Framework for Monte Carlo Simulation Data in Physics

knitr::opts_chunk$set(warning = FALSE, message = FALSE)

This is a companion document for the manuscript arxiv.org/abs/2411.14981. It explains the PGEVM and Lanczos implementations available in the hadron R package. The hadron package version required is version 3.3.1 or newer. It is available on github for download.

You first need to install the hadron package, and load it then as follows

library(hadron)

The code snippets given below can be copied one-to-one to the R command line and executed. Alternatively, you may install the rmarkdown package and render the corresponding Rmd file, which was used to generate this PDF file as follows

render("pgevm.Rmd", output_format="pdf_document")

hadron ships with a pion correlation function, which was also used as the pion example in the aforementioned manuscript.

The Standard Effective Mass analysis

After loading hadron as a library, one first needs to load the pion correlator data

data(pscor.sample)

and store it in an object of class cf (the hadron correlation function class) for further treatment:

pioncf <- cf_orig(cf=t(array(pscor.sample[,2], dim=c(48, 316))))
pioncf <- cf_meta(pioncf, nrObs=1, Time=48, symmetrised=FALSE)

Next, we (block) bootstrap and double bootstrap this correlation function. Note that depending on the choice of boot.R and dbboot.R the execution of the following lines might take long.

boot.R <- 29
dbboot.R <- 19
pioncf.boot <- bootstrap.cf(pioncf, boot.R=boot.R, boot.l=2)
pioncf.boot <- double_bootstrap.cf(pioncf.boot, dbboot.R=dbboot.R)

hadron provides functionality to compute the standard effective mass for a $\cosh$-like correlation function (specified by the type argument to the function bootstrap.effectivemass)

pion.efm <- bootstrap.effectivemass(pioncf.boot, type="solve")
pion.efm <- fit.effectivemass(pion.efm, t1=8, t2=17, useCov=TRUE)

Note that you may get documentation on the different functions e.g. by

?bootstrap.effectivemass

The such computed effective masses and the corresponding fit can be visualised as follows

plot(pion.efm, ylab="Meff", xlab="t/a", xlim=c(0,24), ylim=c(0.1,0.2))

A call to

summary(pion.efm)

would also summarize the fit results.

PGEVM

Based on the (double) bootstrapped correlation function data, the PGEVM is applied as follows:

pion.pgevm <- bootstrap.pgevm(pioncf.boot, N=pioncf.boot$Time-1)
pion.pgevm.efm <- pgevm2effectivemass(pion.pgevm, errortype="dbboot",
                                      bias_correction=TRUE)

Note that you might get one or more warning from the above commands, which are due to the fact that for some combination of data the Hankel matrix cannot be inverted or the eigenvalue compuation fails. But this is taken into account in the analysis by missing values.

The second line translates the PGEVM result into an effective mass object. Thus, plotting it works in an identical manner as for the true effective mass

plot(pion.pgevm.efm, ylab="Meff", xlab="n", xlim=c(0,22), ylim=c(0.1,0.2), col="blue")

Also otherwise, pion.pgevm.efm can be treated like a class of type effectivemass. For instance, one can also perform a plateau fit:

pion.pgevm.efm.fit <- fit.effectivemass(pion.pgevm.efm, t1=10, t2=21, useCov=TRUE)

where now t1 and t2 are interpreted as $n_1$ and $n_2$. The corresponding plot looks as follows

plot(pion.pgevm.efm.fit, ylab="Meff", xlab="n", xlim=c(0,22), ylim=c(0.1,0.2), col="blue")

The fit result reads $E =r pion.pgevm.efm.fit$effmassfit$t0[1] \pm r pion.pgevm.efm.fit$effmassfit$se[1]$.

Lanczos

Lanczos works similarly. It is applied to the double bootstrapped correlation data with the function bootstrap.lanczos as follows

pion.lanczos <- bootstrap.lanczos(pioncf.boot, N=pioncf.boot$Time-1,
                                  bias_correction=TRUE, errortype="dbboot", pivot=TRUE)

and can can be plotted similarly

plot(pion.lanczos, ylab="Meff", xlab="n", xlim=c(0,22), ylim=c(0.1,0.2), col="red")

PGEVM for periodic correlators

pion.pgevm2 <- bootstrap.pgevm(pioncf.boot, N=pioncf.boot$Time-1, ndep.Delta=TRUE)
pion.pgevm.efm2 <- pgevm2effectivemass(pion.pgevm2, errortype="dbboot",
                                       bias_correction=TRUE)
pion.pgevm.efm3 <- pgevm2effectivemass(pion.pgevm2, errortype="dbboot",
                                       bias_correction=TRUE, average.negE=TRUE)

plot(pion.pgevm.efm, ylab="Meff", xlab="n", xlim=c(0,24), ylim=c(0.1,0.2), col="blue", pch=21)
plot(pion.pgevm.efm2, rep=TRUE, col="red", pch=22, xshift=0.2)
plot(pion.pgevm.efm3, rep=TRUE, col="darkorange", pch=23, xshift=0.4)