added javer's modifications

src/juobs_obsmod.jl

+@doc raw"""
+meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bool=false )     
+
+meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false)
+
+Computes effective mass for a given correlator corr at a given plateau plat.
+Correlator can be passed as an Corr struct or Vector{uwreal}.
+    
+The flags pl and data allow to show the plots and return data as an extra result.
+
+```@example
+data = read_mesons(path, "G5", "G5")
+corr_pp = corr_obs.(data)
+m = meff(corr_pp[1], [50, 60], pl=false)
+```
+"""
+function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bool=false, mu::Union{Vector{Float64}, Nothing}=nothing)     
+    dim = length(corr)                                                                                    
+    aux = 0.5 .* log.((corr[2:dim-2] ./ corr[3:dim-1]).^2)     
+    mass = plat_av(aux, plat)                                                                            
+    uwerr(mass)                       
+    if pl == true
+        x = 1:length(aux)
+        y = value.(aux)
+        dy = err.(aux)
+        v = value(mass)
+        e = err(mass)
+        
+        figure()
+        fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75)
+        errorbar(x, y, dy, fmt="x", color="black")
+        ylabel(L"$m_\mathrm{eff}$")
+        xlabel(L"$x_0$")
+    
+        if !isnothing(mu)
+            title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2]))
+        end
+        display(gcf())
+    end                                                  
+    if data == false                   
+        return mass                                        
+    else               
+        return (mass, aux)     
+    end                    
+end   
+meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false) = 
+    meff(corr.obs, plat, pl=pl, data=data, mu=corr.mu)
+
+## Decay constants                                                                                       
+@doc raw"""
+dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false)meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false)
+
+dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false)
+
+Computes decay constant using the PCVC relation for twisted mass fermions. The decay constant is computed in the plateau plat.
+Correlator can be passed as an Corr struct or Vector{uwreal}. If it is passed as a uwreal vector, vector of twisted masses mu and source position y0
+must be specified.
+
+The flags pl and data allow to show the plots and return data as an extra result.
+
+```@example
+data = read_mesons(path, "G5", "G5")
+corr_pp = corr_obs.(data)
+m = meff(corr_pp[1], [50, 60], pl=false)
+f = dec_const_pcvc(corr_pp[1], [50, 60], m, pl=false)
+```
+"""
+function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false)
+    """
+    compute the decay constant when the source is far from the boundaries
+    """
+    corr_pp = corr[2:end-1]
+    dim = length(corr_pp)
+                                                                 
+    aux = exp.((collect(1:dim) .- y0 ) .* [m for k in 1:dim])
+    R = (aux .* corr_pp).^0.5
+    R_av = plat_av(R, plat)
+    f = sqrt(2) * (mu[1] + mu[2]) *R_av / m^1.5
+    uwerr(f)
+    if pl == true
+        uwerr(R_av)
+        v = value(R_av)
+        e = err(R_av)                                
+        uwerr.(R)
+        figure()
+        fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75)                                            
+        errorbar(1:length(R), value.(R), err.(R), fmt="x", color="black")
+        ylabel(L"$R$")
+        xlabel(L"$x_0$")
+        title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2]))
+        display(gcf())          
+    end 
+    if data == false
+        return f
+    else
+        return (f, R)                   
+    end 
+end 
+dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false) =
+    dec_const_pcvc(corr.obs, plat, m, corr.mu, corr.y0, pl=pl, data=data)
+
+#t0
+@doc raw"""
+
+comp_t0(y::YData, plat::Vector{Int64}; pl::Bool=false, L::Int64=1)
+
+comp_t0(Y::Vector{YData}, plat::Vector{Int64}; pl::Bool=false, L::Int64=1)
+
+Computes t0 using the energy density of the action Ysl(Yang-Mills action).
+t0 is computed in the plateau plat.
+A (linear) interpolation in t is performed to find t0.
+
+The flag pl allows to show the plot.
+
+```@example
+#Single replica
+Y = read_ms(path)
+rw = read_ms(path_rw)
+
+t0 = comp_t0(Y, [38, 58], L=32)
+t0_r = comp_t0(Y, [38, 58], L=32, rw=rw)
+
+#Two replicas
+Y1 = read_ms(path1)
+Y2 = read_ms(path2)
+rw1 = read_ms(path_rw1)
+rw2 = read_ms(path_rw2)
+
+t0 = comp_t0([Y1, Y2], [38, 58], L=32, pl=true)
+t0_r = comp_t0(Y, [38, 58], L=32, rw=[rw1, rw2], pl=true)
+
+```
+"""
+function comp_t0(Y::YData, plat::Vector{Int64}; pl::Bool=false, L::Int64=1, rw::Union{Matrix{Float64}, Nothing}=nothing)
+    Ysl = Y.Ysl
+    t = Y.t
+    id = Y.id
+
+    Ysl = isnothing(rw) ? Ysl : apply_rw(Ysl, rw)
+    
+    xmax = size(mean(Ysl, dims=3))[1]
+    nt0 = t0_guess(t, Ysl, plat, L)
+    
+    Y = Matrix{uwreal}(undef, xmax, 3)
+    for i = 1:xmax
+        k = 1
+        for j = nt0-1:nt0+1
+            Y[i, k] = uwreal(Ysl[i, j, :], id)
+            k = k + 1
+        end
+    end
+    x = t[nt0-1:nt0+1]
+    t2E = [plat_av(Y[:, j], plat) for j=1:3] .* x.^2 / L^3
+    
+    par, chi2exp = lin_fit(x, t2E)
+
+    t0 = x_lin_fit(par, 0.3)
+    if pl
+        uwerr(t0)
+        uwerr.(t2E)
+        v = value.(t2E)
+        e = err.(t2E)
+
+        figure()
+        errorbar(x, v, e, fmt="x")
+        
+        errorbar(value(t0), 0.3, xerr=err(t0), fmt="x")
+        ylabel(L"$t^2E$")
+        xlabel(L"$t/a^2$")
+        display(gcf())        
+        end
+    return t0
+end
+
+function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; pl::Bool=false, L::Int64=1, rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing)
+    nr = length(Y)
+    Ysl = getfield.(Y, :Ysl)
+    t = getfield.(Y, :t)
+    id = getfield.(Y, :id)
+    vtr = getfield.(Y, :vtr)
+    replica = length.(vtr)
+    
+
+    if !all(id .== id[1])
+        println("IDs are not equal")
+        return nothing
+    end
+
+    Ysl = isnothing(rw) ? Ysl : apply_rw.(Ysl, rw)
+
+    xmax = size(mean(Ysl[1], dims=3))[1]
+
+    tmp = Ysl[1]
+    [tmp = cat(tmp, Ysl[k], dims=3) for k = 2:nr]
+    nt0 = t0_guess(t[1], tmp, plat, L)
+
+    Y = Matrix{uwreal}(undef, xmax, 3)
+    for i = 1:xmax
+        k = 1
+        for j = nt0-1:nt0+1
+            Y[i, k] = uwreal(tmp[i, j, :], id[1], replica)
+            k = k + 1
+        end
+    end
+    x = t[1][nt0-1:nt0+1]
+    t2E = [plat_av(Y[:, j], plat) for j=1:3] .* x.^2 / L^3
+    println(t2E)
+    par, chi2exp = lin_fit(x, t2E)
+
+    t0 = x_lin_fit(par, 0.3)
+    if pl
+        uwerr(t0)
+        uwerr.(t2E)
+        v = value.(t2E)
+        e = err.(t2E)
+
+        figure()
+        errorbar(x, v, e, fmt="x")
+        
+        errorbar(value(t0), 0.3, xerr=err(t0), fmt="x")
+        ylabel(L"$t^2E$")
+        xlabel(L"$t/a^2$")
+        display(gcf())        
+    end
+    return t0
+end
+function t0_guess(t::Vector{Float64}, Ysl::Array{Float64, 3}, plat::Vector{Int64}, L::Int64)
+    t2E_ax = t.^2 .* mean(mean(Ysl[plat[1]:plat[2], :, :], dims=1), dims=3)[1, :, 1] / L^3
+    t0_aux = minimum(abs.(t2E_ax .- 0.3))
+    nt0 = findfirst(x-> abs(x - 0.3) == t0_aux, t2E_ax) #index closest value to t0
+    return nt0
+end
\ No newline at end of file

src/juobs_readermod.jl

+function read_global_header(path::String)
+    data = open(path, "r")
+    g_header = zeros(Int32, 4)
+    read!(data, g_header)
+    close(data)
+    a = GHeader(g_header)
+    return a
+end
+
+function read_CHeader(path::String)
+    gh = read_global_header(path)
+    
+    data = open(path, "r")
+    seek(data, gh.hsize)
+    
+    aux_f = zeros(Float64, 6)
+    aux_i = zeros(Int32, 4)
+    theta = zeros(Float64, 6)
+    a = Vector{CHeader}(undef, gh.ncorr)
+    for k = 1:gh.ncorr
+        read!(data, aux_f)
+        read!(data, theta)
+
+        qs1 = read(data, Int32)
+        if qs1 != 0
+            qn1 = read(data, Int32)
+            qeps1 = read(data, Float64)
+            q1 = Sm(qs1, qn1, qeps1, 1)
+        else
+            q1 = Sm(qs1, 1)
+        end
+
+        qs2 = read(data, Int32)
+        if qs2 != 0
+            qn2 = read(data, Int32)
+            qeps2 = read(data, Float64)
+            q2 = Sm(qs2, qn2, qeps2, 1)
+        else
+            q2 = Sm(qs2, 1)
+        end
+
+
+        gs1 = read(data, Int32)
+        if gs1 != 0 && gs1 != 3 && gs1 != 4
+            gn1 = read(data, Int32)
+            geps1 = read(data, Float64)
+            g1 = Sm(gs1, gn1, geps1, 2)
+        elseif gs1 == 3 || gs1 == 4
+            g1 = Sm(gs1, q1.niter, q1.eps, 2)
+        else
+            g1 = Sm(gs1, 2)
+        end
+        
+        gs2 = read(data, Int32)
+        if gs2 != 0 && gs2 != 3 && gs2 != 4
+            gn2 = read(data, Int32)
+            geps2 = read(data, Float64)
+            g2 = Sm(gs2, gn2, geps2, 2)
+        elseif gs1 == 3 || gs1 == 4
+            g2 = Sm(gs2, q2.niter, q2.eps, 2)
+        else
+            g2 = Sm(gs2, 2)
+        end
+        
+        read!(data, aux_i)
+        a[k] = CHeader(aux_f, aux_i, theta, [q1, q2, g1, g2])
+    end
+    close(data)
+    return a
+end
+
+@doc raw"""
+read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing)
+
+This faction read a mesons dat file at a given path and returns a vector of CData structures for different masses and Dirac structures.
+Dirac structures g1 and/or g2 can be passed as string arguments in order to filter correaltors.
+ADerrors id can be specified as argument. If is not specified, the id is fixed according to the ensemble name (example: "H400"-> id = 400)
+
+Examples:
+```@example
+read_mesons(path)
+read_mesons(path, "G5")
+read_mesons(path, nothing, "G5")
+read_mesons(path, "G5", "G5")
+read_mesons(path, "G5", "G5", id=1)
+```
+"""
+function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing)
+    
+    isnothing(g1) ? t1=nothing : t1 = findfirst(x-> x==g1, gamma_name) - 1
+    isnothing(g2) ? t2=nothing : t2 = findfirst(x-> x==g2, gamma_name) - 1
+    if isnothing(id)
+        bname = basename(path)
+        m = findfirst(r"[A-Z][0-9]{3}r[0-9]{3}", bname)
+        id = parse(Int64, bname[m[2:4]])
+    end
+
+    data = open(path, "r")
+    g_header = read_global_header(path)
+    c_header = read_CHeader(path)
+
+    ncorr = g_header.ncorr
+    tvals = g_header.tvals
+    nnoise = g_header.nnoise
+
+    fsize = filesize(path)
+
+    datsize = 4 + sum(c.dsize for c in c_header)*tvals*nnoise #datasize / ncnfg
+    ncfg =  Int32((fsize-g_header.hsize-sum(c.hsize for c in c_header)) / datsize)
+
+    corr_match = findall(x-> (x.type1==t1 || isnothing(t1)) && (x.type2==t2 || isnothing(t2)), c_header)
+
+        
+    seek(data, g_header.hsize + sum(c.hsize for c in c_header))
+
+    res = Array{CData}(undef, length(corr_match))
+
+    data_re = Array{Float64}(undef, length(corr_match), ncfg, tvals)
+    data_im = zeros(length(corr_match), ncfg, tvals)
+    vcfg = Array{Int32}(undef, ncfg)
+
+    for icfg = 1:ncfg
+        vcfg[icfg] = read(data, Int32)
+        c=1
+        for k = 1:ncorr
+            if k in corr_match 
+                if c_header[k].is_real == 1
+                    tmp = Array{Float64}(undef, tvals*nnoise)
+                    read!(data, tmp)
+                    tmp2 = reshape(tmp, (nnoise, tvals))
+                    tmp2 = mean(tmp2, dims=1)
+                    data_re[c, icfg, :] = tmp2[1, :]
+                elseif c_header[k].is_real == 0
+                    tmp = Array{Float64}(undef, 2*tvals*nnoise)
+                    read!(data, tmp)
+                    tmp2 = reshape(tmp, (2, nnoise, tvals))
+                    tmp2 = mean(tmp2, dims=2)
+                    data_re[c, icfg, :] = tmp2[1, 1, :]
+                    data_im[c, icfg, :] = tmp2[2, 1, :]
+
+                end
+                c += 1
+            else
+                seek(data, position(data)  + c_header[k].dsize*tvals*nnoise)
+            end
+        
+
+        end
+    end
+
+    for c in 1:length(corr_match)
+        res[c] = CData(c_header[corr_match[c]], vcfg, data_re[c, :, :], data_im[c, :, :], id)
+    end
+    close(data)
+
+    return res
+end
+
+@doc raw"""
+read_ms1(path::String; v::String="1.2")
+
+Reads openQCD ms1 dat files at a given path. This method returns a matrix W[irw, icfg] that contains the reweighting factors, where
+irw is the rwf index and icfg the configuration number.
+The function is compatible with the output files of openQCD v=1.2, 1.4 and 1.6. Version can be specified as argument.
+
+Examples:
+```@example    
+read_ms1(path)
+read_ms1(path, v="1.4")
+read_ms1(path, v="1.6")
+```
+"""
+function read_ms1(path::String; v::String="1.2")
+    data = open(path, "r")
+    nrw = read(data, Int32)
+
+    if v == "1.4" || v =="1.6"
+        nfct = Array{Int32}(undef, nrw)
+        read!(data, nfct)
+        nfct_inheader = 1
+    elseif v=="1.2"
+        nfct = ones(Int32, nrw)
+        nfct_inheader = 0
+    else
+        println("Error: Version not supported")
+        return nothing  
+    end
+    nsrc = Array{Int32}(undef, nrw)
+    read!(data, nsrc)
+    glob_head_size = 4 + 4*nrw*(nfct_inheader + 1)
+    datsize = 4 + 2*8*sum(nsrc .* nfct)
+
+    fsize = filesize(path)
+
+    ncnfg = Int32((fsize - glob_head_size)/datsize)
+    r_data = Array{Array{Float64}}(undef, nrw)
+    
+    [r_data[k] = zeros(Float64, nfct[k], nsrc[k], ncnfg) for k = 1:nrw]
+    vcfg = Array{Int32}(undef, ncnfg)
+    for icfg = 1:ncnfg
+        vcfg[icfg] = read(data, Int32)
+        for irw = 1:nrw
+            for ifct = 1:nfct[irw]
+                tmp = zeros(Float64, nsrc[irw])
+                seek(data, position(data) + 8 * nsrc[irw])
+                read!(data, tmp)
+                r_data[irw][ifct, :, icfg] = tmp[:]
+            end
+        end
+    end
+    
+    W = zeros(Float64, nrw, ncnfg)
+    [W[k, :] = prod(mean(exp.(.-r_data[k]), dims=2), dims=1) for k = 1:nrw]
+    close(data)
+    return W
+end
+
+@doc raw"""
+read_ms(path::String)
+
+Reads openQCD ms dat files at a given path. This method return: 
+    t(t): flow time values
+    Wsl(x0, t, icfg): the time-slice sums of the densities of the Wilson plaquette action
+    Ysl(x0, t, icfg): the time-slice sums of the densities of the Yang-Mills action 
+    Qsl(x0, t, icfg): the time-slice sums of the densities of the topological charge
+
+Examples:
+```@example    
+t, W, Y, Q = read_ms(path)
+```
+"""
+function read_ms(path::String; id::Union{Int64, Nothing}=nothing)  
+    if isnothing(id)
+        bname = basename(path)
+        m = findfirst(r"[A-Z][0-9]{3}r[0-9]{3}", bname)
+        id = parse(Int64, bname[m[2:4]])
+    end  
+    data = open(path, "r")
+
+    dn = read(data, Int32)
+    nn = read(data, Int32)
+    tvals = read(data, Int32)
+    eps = read(data, Float64)
+
+    fsize = filesize(path)
+    datsize=4 + 3*8*(nn + 1) * tvals # measurement size of each cnfg
+
+    ntr = Int32((fsize - 3*4 - 8) / datsize)
+
+    vntr = Vector{Int32}(undef, ntr)
+    
+    # x0, t, cfg 
+    Wsl = Array{Float64}(undef, tvals, nn + 1, ntr) 
+    Ysl = Array{Float64}(undef, tvals, nn + 1, ntr)
+    Qsl = Array{Float64}(undef, tvals, nn + 1, ntr)
+
+    for itr = 1:ntr
+        vntr[itr] = read(data, Int32)
+        for iobs = 1:3
+            for inn = 0:nn
+                tmp = Vector{Float64}(undef, tvals)
+                read!(data, tmp)
+                if iobs == 1
+                    Wsl[:, inn + 1, itr] = tmp
+                elseif iobs == 2
+                    Ysl[:, inn + 1, itr] = tmp
+                elseif iobs == 3
+                    Qsl[:, inn + 1, itr] = tmp
+                end
+
+            end
+        end
+    end
+    close(data)
+    t = Float64.(0:nn) .* dn .* eps
+    return YData(vntr, t, Ysl, id)
+end

src/juobs_typesmod.jl

+ ##############UTILITY NAMES##############
+  #=  STARTING AT 0
+    noise_name=['Z2', 'GAUSS', 'U1', 'Z4']
+    gamma_name = ['G0', 'G1', 'G2', 'G3',
+'invalid', 'G5', '1', 'G0G1', 'G0G2', 'G0G3',
+'G0G5', 'G1G2', 'G1G3', 'G1G5', 'G2G3', 'G2G5', 'G3G5']
+
+    # quark and gauge Smearing types
+    qs_name=['Local', 'Wuppertal', '3D Gradient Flow', 'Gradient Flow']
+    gs_name=['Local', 'APE', '3D Wilson Flow', 'Quark 3D Gradient Flow', 'Quark Gradient Flow']
+=#
+const noise_name=["Z2", "GAUSS", "U1", "Z4"]
+const gamma_name = ["G0", "G1", "G2", "G3",
+"invalid", "G5", "1", "G0G1", "G0G2", "G0G3",
+"G0G5", "G1G2", "G1G3", "G1G5", "G2G3", "G2G5", "G3G5"]
+const qs_name=["Local", "Wuppertal", "3D Gradient Flow", "Gradient Flow"]
+const gs_name=["Local", "APE", "3D Wilson Flow", "Quark 3D Gradient Flow", "Quark Gradient Flow"]
+
+mutable struct GHeader
+    ncorr::Int32
+    nnoise::Int32
+    tvals::Int32
+    noisetype::Int32
+
+    hsize::Int32 #headersize
+    GHeader(a, b, c, d) = new(a, b, c, d, 4*4)
+    function GHeader(x::Vector{Int32})
+        a = new(x[1], x[2], x[3], x[4], 4*4)
+        return a
+    end
+end
+
+mutable struct Sm
+    type::Int32
+    niter::Int32
+    eps::Float64
+    qg::Int32 #1->fermionic 2->gluonic
+    Sm(t,q) = new(t, 0, 0.0, q)
+    Sm(t, n, e, q) = new(t, n, e, q)
+end
+
+mutable struct CHeader
+    k1::Float64
+    k2::Float64
+    mu1::Float64
+    mu2::Float64
+    dp1::Float64
+    dp2::Float64
+    theta1::Vector{Float64}
+    theta2::Vector{Float64}
+    q1::Sm
+    q2::Sm
+    g1::Sm
+    g2::Sm
+    type1::Int32
+    type2::Int32
+    x0::Int32
+    is_real::Int32
+
+    hsize::Int32 #headersize
+    dsize::Int32 #datasize / (nnoise * T * ncfg)
+
+    function CHeader(aux_f::Vector{Float64}, aux_i::Vector{Int32}, theta::Vector{Float64}, sm_par::Vector{Sm})
+        a = new()
+        a.k1 = aux_f[1]
+        a.k2 = aux_f[2]
+        a.mu1 = aux_f[3]
+        a.mu2 = aux_f[4]
+        a.dp1 = aux_f[5]
+        a.dp2 = aux_f[6]
+        a.type1 = aux_i[1]
+        a.type2 = aux_i[2]
+        a.x0 = aux_i[3]
+        a.is_real = aux_i[4]
+
+        a.theta1 = theta[1:3]
+        a.theta2 = theta[4:6]
+        a.q1 = sm_par[1]
+        a.q2 = sm_par[2]
+        a.g1 = sm_par[3]
+        a.g2 = sm_par[4]
+        a.hsize = 8*12 + 4*8 #without smearing parameters niter, neps
+        if a.q1.type != 0
+            a.hsize += 8+4
+        end
+        if a.q2.type != 0
+            a.hsize += 8+4
+        end
+        if a.g1.type != 0 && a.g1.type != 3 && a.g1.type != 4
+            a.hsize += 8+4
+        end
+        if a.g2.type != 0 && a.g2.type != 3 && a.g2.type != 4
+            a.hsize += 8+4
+        end
+        a.dsize = 16 - 8* a.is_real 
+        
+        return a
+    end
+
+end
+
+mutable struct CData
+    header::CHeader
+    vcfg::Array{Int32}
+    re_data::Array{Float64}
+    im_data::Array{Float64}
+    id::Int64
+    CData(a, b, c, d, e) = new(a, b, c, d, e)
+
+end
+Base.copy(a::CData) = CData(a.header, a.vcfg, a.re_data, a.im_data, a.id)
+
+mutable struct Corr
+    obs::Vector{uwreal}
+    mu::Vector{Float64}
+    gamma::Vector{String}
+    y0::Int64
+    function Corr(a::Vector{uwreal}, b::CData)
+        h = getfield(b, :header)
+        mu = [h.mu1, h.mu2]
+        gamma = [gamma_name[h.type1+1], gamma_name[h.type2+1]]
+        y0 = Int64(h.x0)
+        return new(a, mu, gamma, y0)
+    end 
+    function Corr(a::Vector{uwreal}, b::Vector{CData})
+        sym = [:mu1, :mu2, :type1, :type2, :x0]
+        h = getfield.(b, :header)
+        for s in sym
+            if !all(getfield.(h, s) .== getfield(h[1], s))
+                println("Corr: Parameter mismatch")
+                return nothing
+            end
+        end
+        mu = [h[1].mu1, h[1].mu2]
+        gamma = [gamma_name[h[1].type1+1], gamma_name[h[1].type2+1]]
+        y0 = Int64(h[1].x0)
+        return new(a, mu, gamma, y0)
+    end 
+end
+
+mutable struct YData
+    vtr::Vector{Int32}
+    t::Vector{Float64}
+    Ysl::Array{Float64, 3}
+    id::Int64
+    YData(a, b, c, d) = new(a, b, c, d)
+end
+
+function Base.show(io::IO, a::GHeader)
+    print(io, "ncorr = ", a.ncorr, "\t")
+    print(io, "nnoise = ", a.nnoise, "\t")
+    print(io, "tvals = ", a.tvals, "\t")
+    print(io, "noisetype = ", noise_name[a.noisetype + 1], "\t")
+    print(io, "\n")
+
+end
+
+function Base.show(io::IO, a::CHeader)
+
+    fnames = fieldnames(CHeader)
+    for k in fnames
+        f = getfield(a, k)
+        if k != :type1 && k!= :type2
+            print(io, "$k = $f\t")
+        else
+            print(io, "$k = ", gamma_name[f + 1], "\t")   
+        end
+    end
+    print(io, "\n")
+end
+
+function Base.show(io::IO, a::Sm)
+    print(io, "\n")
+    if a.qg == 1
+        print(io, "type = ", qs_name[a.type+1], "\t")
+    elseif a.qg == 2
+        print(io, "type = ", gs_name[a.type+1], "\t")
+    end
+    print(io, "niter = ", a.niter, "\t")
+    print(io, "eps = ", a.eps, "\t")
+    print(io, "\n")
+end
+
+function Base.show(io::IO, a::CData)
+    print(io, a.header)
+end
\ No newline at end of file

src/juobsmod.jl

+module juobs
+using ADerrors, PyPlot, Statistics, LaTeXStrings, LinearAlgebra, LsqFit
+
+include("juobs_types.jl")
+include("juobs_linalg.jl")
+include("juobs_reader.jl")
+include("juobs_tools.jl")
+include("juobs_obs.jl")
+
+export read_mesons, read_ms1, read_ms
+export get_matrix, uwgevp_tot, energies, uwdot
+export corr_obs, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
+export meff, dec_const_pcvc, comp_t0
+
+end # module