168 KiB
import numpy as np
import matplotlib.pyplot as plt
import pyerrors as pe
plt.style.use('./base_style.mplstyle')
plt.rc('text', usetex=True)
We can load data from a preprocessed file which contains a list of pyerror Obs:
correlator_data = pe.input.json.load_json("./data/correlator_test")
With this list a Corr object can be initialised
my_correlator = pe.Corr(correlator_data)
my_correlator.print([8, 14])
The show method can display the correlator
my_correlator.show()
Manipulating correlators¶
Corr objects can be shifted
shifted_correlator = my_correlator.roll(20)
shifted_correlator.tag = r'Correlator shifted by $x_0/a=20$'
or symmetrised
symmetrised_correlator = my_correlator.symmetric()
symmetrised_correlator.tag = 'Symmetrised correlator'
We can compare different Corr objects by passing comp to the show method
shifted_correlator.show(comp=symmetrised_correlator, logscale=True)
Effective mass¶
The effective mass of the correlator can be obtained by calling the m_eff method
m_eff = symmetrised_correlator.m_eff()
m_eff.tag = 'Effective mass'
We can also use the priodicity of the lattice in order to obtain the cosh effective mass
periodic_m_eff = symmetrised_correlator.m_eff('periodic')
periodic_m_eff.tag = 'Cosh effective mass'
We can compare the two and see how the standard effective mass deviates form the plateau at the center of the lattice
periodic_m_eff.show([4,47], comp=m_eff, ylabel=r'$am_\mathrm{eff}$')
Arithmetic operations and mathematical functions are also overloaded for the Corr class. We can compute the difference between the two variants of the effective mass as follows.
difference_m_eff = np.abs(periodic_m_eff - m_eff)
difference_m_eff.show([0, 47], logscale=True)
Derivatives¶
We can obtain derivatives of correlators in the following way
first_derivative = symmetrised_correlator.deriv()
first_derivative.tag = 'First derivative'
second_derivative = symmetrised_correlator.second_deriv()
second_derivative.tag = 'Second derivative'
symmetrised_correlator.show([5, 20], comp=[first_derivative, second_derivative], y_range=[-500, 1300])
Missing Values¶
Apart from the build-in functions, there is another reason, why one should use a Corr instead of a list of Obs. Missing values are handled for you. We will create a second correlator with missing values.
new_content=[(my_correlator.content[i] if i not in [6,8,9,12,14,15,20] else None ) for i in range(my_correlator.T) ] # We reuse the old example and replace a few values with None
correlator_incomplete=pe.Corr(new_content)
correlator_incomplete
We see that this is still a valid correlator. It is just missing some values. When we perform operations, which generate new correlators, the missing values are handled automatically.
Here is an example:
some_new_corr=np.sin(my_correlator+2*correlator_incomplete)
some_new_corr
Some functions might also return correlators with missing values. We already looked at the forward derivative. The forward derivative is not defined for the last value.
The important thing is that, whatever you do, correlators keep their length T. So there will never be confusion about how you count timeslices.
assert first_derivative.T==my_correlator.T==len(first_derivative.content)==len(my_correlator.content)
assert first_derivative.content[-1]==None
You can also take a plateau or perform a fit, even though some values might be missing.
There is a range of addtional methods of the Corr class which can be found in the documentation.
Matrix Correlators¶
A correlator can not only contain a list of Obs, but also a list of matrices of obs. This is useful, if there are multiple sources and sinks used. In our example, the sources have a different Gaussian smearing applied.
We will load such a correlator.
matrix_V1V1= pe.input.json.load_json("./data/matrix_correlator_V1V1")
print(matrix_V1V1.content[0])
We printed out the content at timeslice 0. As we can see, it is a matrix of Obs.
Let us try to get the effective mass.
try:
matrix_V1V1.m_eff() #This does not work!
except:
print("Something is wrong")
Many methods we could use for regular correlators do not work with matrix-correlators. In order to get the effective mass, we need to convert to a regular correlator first.
One way to do it, is to pick a smearing out of the matrix:
single_smearing=matrix_V1V1.smearing(0,0)
single_smearing
.smearing(i,j) picks the element [i,j] from every matrix and returns a correlator containing one Obs per timeslice. But there is a more usefull way to retrieve a single value per timeslice. We might want a linear combination of different sources and sinks. We can formalize this as
$$C_{\textrm{projected}}(t)=v_1^T \underline{C}(t) v_2$$If we choose the vectors to be $v_1=v_2=(0,1,0,0)$, we should get the same correlator as in the cell above.
Thinking about it this way is usefull in the Context of the generalized eigenvalue problem (GEVP), used to find the source-sink combination, which best describes a certain energy eigenstate. A good introduction is found in https://arxiv.org/abs/0902.1265.
vec=matrix_V1V1.GEVP(t0=3,ts=6,state=0)
assert len(vec)==matrix_V1V1.N
vec
As we see, the eigenvector is of length matrix_V1V1.N and contains regular floats. We can use it to project the correlator.
matrix_V1V1.projected(vec).m_eff().show(comp=single_smearing.m_eff())
There is a lot going on in this line of code. We start with our matrix correlator and we project it, using the vector we got from the GEVP routine.
This gives us a new correlator with one Obs per timeslice. We then calculate its effective mass and plot it. We tell the .plot() function to show another correlator as a comparison.
We can see, that the projected correlator (blue) converges to a mass plateau much quicker than the single smearing.
Example Analysis¶
We can use what we learned so far to perform an actually usefull analysis. The correlator matrix_V1V1 we looked at corresponds to vector-charmonium.
We might be interested in the mass of the $J/\Psi$ state.
# We do not just have V1V1, but also the two other spacial directions. We can average over them for better statistics.
matrix_V2V2= pe.input.json.load_json("./data/matrix_correlator_V2V2")
matrix_V3V3= pe.input.json.load_json("./data/matrix_correlator_V3V3")
matrix_VnVn=(matrix_V1V1+matrix_V2V2+matrix_V3V3)/3.
#We then solve the GEVP to get eigenvectors corresponding to the ground state.
vec_ground=matrix_VnVn.GEVP(t0=3,ts=6,state=0)
#Now we project the matrix-correlators to get new correlators belonging to the ground state.
corr_ground=matrix_VnVn.projected(vec_ground)
# We get the effective mass using the periodic cosh method.
m_eff_Jpsi=corr_ground.m_eff(variant="cosh")
m_eff_Jpsi.show([5,25])
#From the plot we can pick a plateau range and get a single value for the mass.
m_Jpsi=m_eff_Jpsi.plateau([8,18])
# Since the lattice spacing is known, we can multiply with hbar*c/a to see it in physical units
m_Jpsi=m_Jpsi*197/0.0653
#As a last step we call the gamma method to get the error
m_Jpsi.gamma_method()
print("--- The mass was calculated to be" , m_Jpsi, "MeV ---")