Merge branch 'javier'

29ae8c09 · Javier · 08e2b69f · 3467bd56 · 29ae8c09 · 29ae8c09
Commit 29ae8c09 authored Mar 23, 2021 by Javier
10 changed files
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,5 @@
 analysis/plat\.txt
 *.pdf
+analysis/*.txt
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -2,17 +2,23 @@
 [[ADerrors]]
 deps = ["BDIO", "FFTW", "FastGaussQuadrature", "ForwardDiff", "LaTeXStrings", "LinearAlgebra", "Nettle", "PGFPlotsX", "Plots", "Printf", "QuadGK", "Roots", "Statistics", "Test", "UnicodePlots"]
-git-tree-sha1 = "f1c9954c17c0e87f185d532163e9e630efd12fe3"
+git-tree-sha1 = "2d3b701de1b206f497c632413b5a783fbd65d714"
-repo-rev = "9c41ec2e161e4826bfe0752be8bf2a0bc025759c"
+repo-rev = "master"
 repo-url = "https://gitlab.ift.uam-csic.es/alberto/aderrors.jl"
 uuid = "5e92007d-7bf1-471c-8ceb-4591b8b567a9"
 version = "0.1.0"
 [[AbstractFFTs]]
 deps = ["LinearAlgebra"]
-git-tree-sha1 = "8ed9de2f1b1a9b1dee48582ad477c6e67b83eb2c"
+git-tree-sha1 = "485ee0867925449198280d4af84bdb46a2a404d0"
 uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
-version = "1.0.0"
+version = "1.0.1"
+[[Adapt]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "ffcfa2d345aaee0ef3d8346a073d5dd03c983ebe"
+uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+version = "3.2.0"
 [[ArgCheck]]
 git-tree-sha1 = "dedbbb2ddb876f899585c4ec4433265e3017215a"
@@ -197,9 +203,9 @@ version = "4.3.1+4"
 [[FFTW]]
 deps = ["AbstractFFTs", "FFTW_jll", "IntelOpenMP_jll", "Libdl", "LinearAlgebra", "MKL_jll", "Reexport"]
-git-tree-sha1 = "c31e446bf3b12aad2ec8fc500fe19528c148d811"
+git-tree-sha1 = "1b48dbde42f307e48685fa9213d8b9f8c0d87594"
 uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
-version = "1.3.1"
+version = "1.3.2"
 [[FFTW_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -275,9 +281,9 @@ version = "6.1.2+6"
 [[GR]]
 deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "LinearAlgebra", "Pkg", "Printf", "Random", "Serialization", "Sockets", "Test", "UUIDs"]
-git-tree-sha1 = "b90b826782cb3ac5b7a7f41b3fd0113180257ed4"
+git-tree-sha1 = "12d971c928b7ecf19b748a2c7df6a365690dbf2c"
 uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
-version = "0.53.0"
+version = "0.55.0"
 [[GR_jll]]
 deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt_jll", "Zlib_jll", "libpng_jll"]
@@ -287,9 +293,9 @@ version = "0.53.0+0"
 [[GeometryBasics]]
 deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
-git-tree-sha1 = "a28d728c2d825285fe27f38ca322399d35d1a5b9"
+git-tree-sha1 = "4d4f72691933d5b6ee1ff20e27a102c3ae99d123"
 uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
-version = "0.3.6"
+version = "0.3.9"
 [[Gettext_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
@@ -309,10 +315,10 @@ uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe"
 version = "1.0.0"
 [[HTTP]]
-deps = ["Base64", "Dates", "IniFile", "MbedTLS", "Sockets"]
+deps = ["Base64", "Dates", "IniFile", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
-git-tree-sha1 = "c7ec02c4c6a039a98a15f955462cd7aea5df4508"
+git-tree-sha1 = "c9f380c76d8aaa1fa7ea9cf97bddbc0d5b15adc2"
 uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
-version = "0.8.19"
+version = "0.9.5"
 [[IOCapture]]
 deps = ["Logging"]
@@ -386,6 +392,18 @@ git-tree-sha1 = "3a0084cec7bf157edcb45a67fac0647f88fe5eaf"
 uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 version = "0.14.7"
+[[LazyArtifacts]]
+deps = ["Pkg"]
+git-tree-sha1 = "4bb5499a1fc437342ea9ab7e319ede5a457c0968"
+uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
+version = "1.3.0"
+[[LeastSquaresOptim]]
+deps = ["FiniteDiff", "ForwardDiff", "LinearAlgebra", "Optim", "Printf", "SparseArrays", "Statistics", "SuiteSparse"]
+git-tree-sha1 = "b5a1931bf37616820c4ef71629b6308d67f0b393"
+uuid = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891"
+version = "0.8.1"
 [[LibGit2]]
 deps = ["Printf"]
 uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
@@ -467,10 +485,10 @@ uuid = "2fda8390-95c7-5789-9bda-21331edee243"
 version = "0.12.0"
 [[MKL_jll]]
-deps = ["IntelOpenMP_jll", "Libdl", "Pkg"]
+deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
-git-tree-sha1 = "eb540ede3aabb8284cb482aa41d00d6ca850b1f8"
+git-tree-sha1 = "c253236b0ed414624b083e6b72bfe891fbd2c7af"
 uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
-version = "2020.2.254+0"
+version = "2021.1.1+1"
 [[MacroTools]]
 deps = ["Markdown", "Random"]
@@ -531,6 +549,11 @@ git-tree-sha1 = "b3b169cec36b850cfacbefc0b0a5f68a6089eb70"
 uuid = "4c82536e-c426-54e4-b420-14f461c4ed8b"
 version = "3.4.1+2"
+[[NetworkOptions]]
+git-tree-sha1 = "ed3157f48a05543cce9b241e1f2815f7e843d96e"
+uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+version = "1.2.0"
 [[Ogg_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "a42c0f138b9ebe8b58eba2271c5053773bde52d0"
@@ -614,9 +637,9 @@ uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 [[PlotThemes]]
 deps = ["PlotUtils", "Requires", "Statistics"]
-git-tree-sha1 = "c6f5ea535551b3b16835134697f0c65d06c94b91"
+git-tree-sha1 = "a3a964ce9dc7898193536002a6dd892b1b5a6f1d"
 uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
-version = "2.0.0"
+version = "2.0.1"
 [[PlotUtils]]
 deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "Statistics"]
@@ -626,9 +649,9 @@ version = "1.0.10"
 [[Plots]]
 deps = ["Base64", "Contour", "Dates", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "JSON", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs"]
-git-tree-sha1 = "4797acb266b8d9ff316f4581924e71c6709f152d"
+git-tree-sha1 = "142dd04f5060c04de91cc10ca76dffb291a02426"
 uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-version = "1.10.1"
+version = "1.10.6"
 [[PositiveFactorizations]]
 deps = ["LinearAlgebra"]
@@ -679,9 +702,9 @@ version = "1.1.1"
 [[RecipesPipeline]]
 deps = ["Dates", "NaNMath", "PlotUtils", "RecipesBase"]
-git-tree-sha1 = "9ea2f5bf1b26918b16e9f885bb8e05206bfc2144"
+git-tree-sha1 = "c4d54a78e287de7ec73bbc928ce5eb3c60f80b24"
 uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c"
-version = "0.2.1"
+version = "0.3.1"
 [[Reexport]]
 deps = ["Pkg"]
@@ -709,9 +732,9 @@ version = "0.2.2+1"
 [[Roots]]
 deps = ["Printf"]
-git-tree-sha1 = "8f743e4f4368d1d753f3806bf635899dad6b4847"
+git-tree-sha1 = "369e25546984dff5df351bc056fccc30de615080"
 uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
-version = "1.0.7"
+version = "1.0.8"
 [[SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
@@ -777,10 +800,10 @@ uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 version = "0.9.6"
 [[StructArrays]]
-deps = ["DataAPI", "Tables"]
+deps = ["Adapt", "DataAPI", "Tables"]
-git-tree-sha1 = "ad1f5fd155426dcc879ec6ede9f74eb3a2d582df"
+git-tree-sha1 = "26ea43b4be7e919a2390c3c0f824e7eb4fc19a0a"
 uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
-version = "0.4.2"
+version = "0.5.0"
 [[SuiteSparse]]
 deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
@@ -794,14 +817,19 @@ version = "1.0.0"
 [[Tables]]
 deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
-git-tree-sha1 = "240d19b8762006ff04b967bdd833269ad642d550"
+git-tree-sha1 = "a716dde43d57fa537a19058d044b495301ba6565"
 uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
-version = "1.2.2"
+version = "1.3.2"
 [[Test]]
 deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
 uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+[[URIs]]
+git-tree-sha1 = "7855809b88d7b16e9b029afd17880930626f54a2"
+uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
+version = "1.2.0"
 [[UUIDs]]
 deps = ["Random", "SHA"]
 uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
@@ -983,9 +1011,9 @@ version = "1.2.11+18"
 [[Zstd_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "6f1abcb0c44f184690912aa4b0ba861dd64f11b9"
+git-tree-sha1 = "2c1332c54931e83f8f94d310fa447fd743e8d600"
 uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
-version = "1.4.5+2"
+version = "1.4.8+0"
 [[libass_jll]]
 deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]

--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ version = "0.1.0"
 ADerrors = "5e92007d-7bf1-471c-8ceb-4591b8b567a9"
 BDIO = "375f315e-f2c4-11e9-2ef9-134f02f79e27"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
+LeastSquaresOptim = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"

--- a/analysis/analysis_pcac.jl
+++ b/analysis/analysis_pcac.jl
+using juobs, ADerrors, DelimitedFiles, PyPlot, LaTeXStrings
+const path = "/home/javier/Lattice/charm/production_2"
+const path_plat = "/home/javier/Lattice/juobs/analysis/plat.txt"
+const ensembles = ["H400", "N202", "N200", "N203", "N300", "J303"]
+const deg = [true, true, false, false, true, false]
+const L = [32, 48, 48, 48, 48, 64]
+const beta = [3.46 , 3.55, 3.55, 3.55, 3.70, 3.70]
+const R = [["H400r001", "H400r002"], "N202r001", ["N200r000", "N200r001"], ["N203r000", "N203r001"], "N300r002", "J303r003"]
+include("/home/javier/Lattice/juobs/analysis/functions.jl")
+include("/home/javier/Lattice/juobs/constants/juobs_const.jl")
+m_ll = Vector{uwreal}(undef, length(ensembles))
+m_ss = Vector{uwreal}(undef, length(ensembles))
+m_hh = Vector{Vector{uwreal}}(undef, length(ensembles))
+mu_pp = Vector{Vector{Vector{Float64}}}(undef, length(ensembles))
+for iens = 1:length(ensembles)
+    pp = read_dat(R[iens], "G5", "G5")
+    a0p = read_dat(R[iens], "G5", "G0G5")
+    rew = read_rew(R[iens])
+    pp_obs = corr_obs.(pp, L=L[iens], rw=rew)
+    a0p_obs = corr_obs.(a0p, L=L[iens], rw=rew)
+    mu_pp[iens] = getfield.(pp_obs, :mu)
+    m = comp_pcac(a0p_obs, pp_obs, deg[iens], ensembles[iens])
+    m_ll[iens] = get_ll(mu_pp[iens], m, deg[iens])
+    m_ss[iens] = deg[iens] ? m_ll[iens] : get_ss(mu_pp[iens], m, deg[iens])
+    m_hh[iens] = get_hh(mu_pp[iens], m, deg[iens])
+end
+mm = get_mu.(mu_pp, deg)
+mul_pp = getindex.(mm, 1)
+mus_pp = getindex.(mm, 2)
+muh_pp = getindex.(mm, 3)
+cot_ll = Vector{uwreal}(undef, length(ensembles))
+cot_ss = Vector{uwreal}(undef, length(ensembles))
+cot_hh = Vector{Vector{uwreal}}(undef, length(ensembles))
+for iens = 1:length(ensembles)
+    cot_ll[iens] = za(beta[iens]) * m_ll[iens] / mul_pp[iens]
+    cot_ss[iens] = za(beta[iens]) * m_ss[iens] / mus_pp[iens]
+    cot_hh[iens] = fill(za(beta[iens]), 3) .* m_hh[iens] ./ muh_pp[iens]
+end
\ No newline at end of file
--- a/analysis/functions.jl
+++ b/analysis/functions.jl
@@ -8,14 +8,19 @@ function read_dat(rep::Vector{String}, g1::String="G5", g2::String="G5")
    p = joinpath.(path, rep, "info")
    f = [filter(x-> contains(x, ".dat"), readdir(p[k]))[1] for k = 1:length(p)] 
    f = joinpath.(p, f)
-    aux = read_mesons.(f, g1, g2)
+    return read_mesons(f, g1, g2)
-    res = []
+end
-    for j = 1:length(aux[1])
+function read_rew(rep::String)
-        aux2 = [aux[i][j] for i = 1:length(aux)]
+    p = joinpath(path, rep, "rew")
-        push!(res, aux2)
+    f = filter(x-> contains(x, ".dat"), readdir(p))
+    f = joinpath(p, f[1])
+    try
+        return read_ms1(f)
+    catch
+        return read_ms1(f, v="1.4")
    end
-    return res
 end
+read_rew(rep::Vector{String}) = read_rew.(rep)
 function get_mu(mu_list::Vector{Vector{Float64}}, deg::Bool)
    mu_sorted = unique(sort(minimum.(mu_list)))
@@ -48,6 +53,16 @@ function get_ls(mu_list, obs::Vector{uwreal}, deg::Bool)
    end
 end
+function get_ss(mu_list, obs::Vector{uwreal}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mus in mu && mu[1] == mu[2] #s-s
+            return obs[k]
+        end
+    end
+end
 function get_lh(mu_list, obs::Vector{uwreal}, deg::Bool)
    mul, mus, muh = get_mu(mu_list, deg)
    obs_lh = Vector{uwreal}(undef, 0)
@@ -72,6 +87,17 @@ function get_sh(mu_list, obs::Vector{uwreal}, deg::Bool)
    return obs_sh
 end
+function get_hh(mu_list, obs::Vector{uwreal}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    obs_hh = Vector{uwreal}(undef, 0)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mu[1] == mu[2] && !(mul in mu) && !(mus in mu)
+            push!(obs_hh, obs[k])
+        end
+    end
+    return obs_hh
+end
 function select_plateau(ens::String, mu_list, deg::Bool)
    mul, mus, muh = get_mu(mu_list, deg)
    f = readdlm(path_plat)
@@ -116,6 +142,11 @@ function comp_meff(pp_obs::Vector{juobs.Corr}, deg::Bool, ens::String; pl::Bool=
    plat = select_plateau(ens, mu_list, deg)
    return meff.(pp_obs, plat, pl=pl)
 end
+function comp_pcac(a0p_obs::Vector{juobs.Corr}, pp_obs::Vector{juobs.Corr}, deg::Bool, ens::String; pl::Bool=false)
+    mu_list = getfield.(pp_obs, :mu)
+    plat = select_plateau(ens, mu_list, deg)
+    return mpcac.(a0p_obs, pp_obs, plat, pl=pl)
+end
 function comp_f(pp_obs::Vector{juobs.Corr}, m::Vector{uwreal}, deg::Bool, ens::String; pl::Bool=false)
    mu_list = getfield.(pp_obs, :mu)
    plat = select_plateau(ens, mu_list, deg)

--- a/src/juobs.jl
+++ b/src/juobs.jl
 module juobs
-using ADerrors, PyPlot, LaTeXStrings, LinearAlgebra, LsqFit, Optim
+using ADerrors, PyPlot, LaTeXStrings, LinearAlgebra, LsqFit, LeastSquaresOptim
 import Statistics: mean
 include("juobs_types.jl")
@@ -11,6 +11,6 @@ include("juobs_obs.jl")
 export read_mesons, read_ms1, read_ms, read_md, truncate_data!
 export get_matrix, energies, uwdot, uweigvals, uweigvecs, uweigen, invert, getall_eigvals, getall_eigvecs
 export corr_obs, md_sea, md_val, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
-export meff, dec_const_pcvc, comp_t0
+export meff, mpcac, dec_const, dec_const_pcvc, comp_t0
 end # module
--- a/src/juobs_obs.jl
+++ b/src/juobs_obs.jl
--- a/src/juobs_reader.jl
+++ b/src/juobs_reader.jl
@@ -81,13 +81,13 @@ function read_CHeader(path::String; legacy::Bool=false)
 end
 @doc raw"""
-    read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing, legacy::Bool=false)
+    read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{String, Nothing}=nothing, legacy::Bool=false)
-    read_mesons(path::Vector{String}, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing, legacy::Bool=false)
+    read_mesons(path::Vector{String}, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{String, Nothing}=nothing, legacy::Bool=false)
 This function 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)
+ADerrors id can be specified as argument. If is not specified, the `id` is fixed according to the ensemble name (example: "H400"-> id = "H400")
 *For the old version (without smearing, distance preconditioning and theta) set legacy=true.
@@ -97,17 +97,18 @@ read_mesons(path)
 read_mesons(path, "G5")
 read_mesons(path, nothing, "G5")
 read_mesons(path, "G5", "G5")
-read_mesons(path, "G5", "G5", id=1)
+read_mesons(path, "G5", "G5", id="H100")
 read_mesons(path, "G5_d2", "G5_d2", legacy=true)
 ```
 """
-function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing, legacy::Bool=false)    
+function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{String, Nothing}=nothing, legacy::Bool=false)    
    t1 = isnothing(g1) ? nothing : findfirst(x-> x==g1, gamma_name) - 1
    t2 = isnothing(g2) ? nothing : 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]])
+        id = bname[m[1:4]]
+        #id = parse(Int64, bname[m[2:4]])
    end
    data = open(path, "r")
@@ -120,8 +121,8 @@ function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union
    fsize = filesize(path)
-    datsize = 4 + sum(c.dsize for c in c_header)*tvals*nnoise #datasize / ncnfg
+    datsize = 4 + sum(getfield.(c_header, :dsize)) * tvals * nnoise #data_size / ncnfg
-    ncfg =  Int32((fsize-g_header.hsize-sum(c.hsize for c in c_header)) / datsize)
+    ncfg = div(fsize - g_header.hsize - sum(getfield.(c_header, :hsize)), datsize) #(total size - header_size) / data_size
    corr_match = findall(x-> (x.type1==t1 || isnothing(t1)) && (x.type2==t2 || isnothing(t2)), c_header)
@@ -171,7 +172,7 @@ function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union
    return res
 end
-function read_mesons(path::Vector{String}, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing, legacy::Bool=false)
+function read_mesons(path::Vector{String}, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{String, Nothing}=nothing, legacy::Bool=false)
    res = read_mesons.(path, g1, g2, id=id, legacy=legacy)
    nrep = length(res)
    ncorr = length(res[1])
@@ -274,28 +275,32 @@ function read_md(path::String)
 end
 @doc raw"""
-    read_ms(path::String; id::Union{Int64, Nothing}=nothing, dtr::Int64=1)  
+    read_ms(path::String; id::Union{String, Nothing}=nothing, dtr::Int64=1, obs::String="Y")  
 Reads openQCD ms dat files at a given path. This method return YData: 
 - `t(t)`: flow time values
- `Ysl(icfg, x0, t)`: the time-slice sums of the densities of the Yang-Mills action 
+- `obs(icfg, x0, t)`: the time-slice sums of the densities of the observable (Wsl, Ysl or Qsl)
 - `vtr`: vector that contains trajectory number
 - `id`: ensmble id
+`dtr` = `dtr_cnfg` / `dtr_ms`, where `dtr_cnfg` is the number of trajectories computed before saving the configuration. `dtr_ms`
+is the same but applied to the ms.dat file.
 Examples:
 ```@example    
 Y = read_ms(path)
 ```
 """
-function read_ms(path::String; id::Union{Int64, Nothing}=nothing, dtr::Int64=1)  
+function read_ms(path::String; id::Union{String, Nothing}=nothing, dtr::Int64=1 , obs::String="Y")  
    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]])
+        id = bname[m[1:4]]
+        #id = parse(Int64, bname[m[2:4]])
    end  
    data = open(path, "r")
@@ -346,7 +351,17 @@ function read_ms(path::String; id::Union{Int64, Nothing}=nothing, dtr::Int64=1)
    end
    close(data)
    t = Float64.(0:nn) .* dn .* eps
-    return YData(vntr, t, Ysl, id)
+    if obs == "W"
+        return YData(vntr, t, Wsl, id)
+    elseif obs == "Y"
+        return YData(vntr, t, Ysl, id)
+    elseif obs == "Q"
+        return YData(vntr, t, Qsl, id)
+    else
+        println("obs = ", obs," is not valid")
+        return nothing
+    end
 end
 @doc raw"""
@@ -377,7 +392,7 @@ truncate_data!(Y, [nc1, nc2])
 """
 function truncate_data!(data::YData, nc::Int64)
    data.vtr = data.vtr[1:nc]
-    data.Ysl = data.Ysl[1:nc, :, :]
+    data.obs = data.obs[1:nc, :, :]
    return nothing
 end
 function truncate_data!(data::Vector{YData}, nc::Vector{Int64})

--- a/src/juobs_tools.jl
+++ b/src/juobs_tools.jl
@@ -29,21 +29,6 @@ function apply_rw(data::Vector{<:Array{Float64}}, W::Vector{Matrix{Float64}})
    return data_r
 end
-function dobsdp(a::uwreal, p::uwreal) # Compute da / dp
-    if count(p.prop .== true) != 1
-        error("I do not know how to compute this")
-    end
-    for i = 1:min(length(a.der), length(p.der))
-        if (p.prop[i] && a.prop[i])
-            return a.der[i] / p.der[i]
-        end 
-    end 
-    return 0.0
-end
 function check_corr_der(obs::Corr, derm::Vector{Corr})   
    g1 = Vector{String}(undef, 0)
    g2 = Vector{String}(undef, 0)
@@ -141,7 +126,7 @@ function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array
    return Corr(obs, cdata)
 end    
-#TODO: VECTORIZE
+#TODO: VECTORIZE, uwreal?
 @doc raw"""
    md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADerrors.wsg)
@@ -195,8 +180,8 @@ function md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADer
    if !all(id .== id[1])
        error("ids do not match")
    end
-    id = id[1]
+    id = ws.id2str[id[1]]
    ivrep = getfield.(ws.fluc[p], :ivrep)
    ivrep1 = fill(ivrep[1], length(ivrep))
    if !all(ivrep .== ivrep1)
@@ -208,6 +193,7 @@ function md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADer
        error("Nr obs != Nr md")
    end
+    #md_aux as a Matrix + Automatic truncation
    md_aux = md[1][:, 1:ivrep[1]]
    for k = 2:length(md)
        md_aux = cat(md_aux, md[k][:, 1:ivrep[k]], dims=2)
@@ -215,27 +201,21 @@ function md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADer
    fluc_obs = getfield.(ws.fluc[p], :delta)
    fluc_md = md_aux .- mean(md_aux, dims=2)
+    uwerr(a)
+    fluc_obs = mchist(a, id)
    nrw = size(fluc_md, 1)
-    d = a.der[p]
-    nobs = sum(a.prop)
    if nrw == 1
-        der1 = uwreal(0)
+        der1 = uwreal(-fluc_md[1, :] .* fluc_obs, id, ivrep)
-        for k = 1:nobs
-            der1 = der1 - d[k] * uwreal(fluc_md[1, :] .* fluc_obs[k], id, ivrep)
-        end
        return (der1, der1)
    elseif nrw == 2
-        der1 = uwreal(0)
+        der1 = uwreal(-fluc_md[1, :] .* fluc_obs, id, ivrep)
-        der2 = uwreal(0)
+        der2 = uwreal(-fluc_md[2, :] .* fluc_obs, id, ivrep)
-        for k = 1:nobs
-            der1 = der1 - d[k] * uwreal(fluc_md[1, :] .* fluc_obs[k], id, ivrep)
-            der2 = der2 - d[k] * uwreal(fluc_md[2, :] .* fluc_obs[k], id, ivrep)
-        end
        return (der1, der2)
    else
        return nothing
    end
 end
 @doc raw"""
@@ -286,7 +266,7 @@ function md_val(a::uwreal, obs::Corr, derm::Vector{Corr})
    corr = getfield(obs, :obs)
-    der = [dobsdp(a, corr[k]) for k = 1:length(corr)]
+    der = [derivative(a, corr[k]) for k = 1:length(corr)]
    derm1, derm2 = derm
    return (sum(der .* derm1.obs), sum(der .* derm2.obs))
 end
@@ -442,14 +422,18 @@ function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}
        C = isnothing(wpm) ? [ADerrors.cov(aux[k]) for k = 1:Ndata] : [ADerrors.cov(aux[k], wpm) for k = 1:Ndata]
        chisq_full_cov(p, d) = get_chi2_cov(model, d, C, p, Nalpha)
        min_fun_cov(t) = chisq_full_cov(t, dat)
-        sol = optimize(min_fun_cov, vcat(fit.param, dat[1:Nalpha*Ndata]), method=LBFGS())
+        sol = optimize(min_fun_cov, vcat(fit.param, dat[1:Nalpha*Ndata]), LevenbergMarquardt())
-        (upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq_full_cov, Optim.minimizer(sol), data) : fit_error(chisq_full_cov, Optim.minimizer(sol), data, wpm)
+        (upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq_full_cov, sol.minimizer, data) : fit_error(chisq_full_cov, sol.minimizer, data, wpm)
+        println("Chisq / chiexp: ", min_fun_cov(sol.minimizer), " / ", chi2_exp, " (dof: ", length(ydata) - param,")")
    else
        chisq_full(p, d) = get_chi2(model, d, ddat, p, Nalpha)
        min_fun(t) = chisq_full(t, dat)
-        sol = optimize(min_fun, vcat(fit.param, dat[1:Nalpha*Ndata]), method=LBFGS())
+        sol = optimize(min_fun, vcat(fit.param, dat[1:Nalpha*Ndata]), LevenbergMarquardt())
-        (upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq_full, Optim.minimizer(sol), data) : fit_error(chisq_full, Optim.minimizer(sol), data, wpm)
+        (upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq_full, sol.minimizer, data) : fit_error(chisq_full, sol.minimizer, data, wpm)
+        println("Chisq / chiexp: ", min_fun(sol.minimizer), " / ", chi2_exp, " (dof: ", length(ydata) - param,")")
    end
    #### chisq_full, min_fun out of conditional ->
@@ -461,7 +445,6 @@ function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}
        print("\n Fit parameter: ", i, ": ")
        details(upar[i])
    end
-    println("Chisq / chiexp: ", sol.minimum, " / ", chi2_exp, " (dof: ", length(ydata) - param,")")
    return upar
 end

--- a/src/juobs_types.jl
+++ b/src/juobs_types.jl
@@ -133,7 +133,7 @@ mutable struct CData
    vcfg::Array{Int32}
    re_data::Array{Float64}
    im_data::Array{Float64}
-    id::Int64
+    id::String
    CData(a, b, c, d, e) = new(a, b, c, d, e)
 end
@@ -172,8 +172,8 @@ end
 mutable struct YData
    vtr::Vector{Int32}
    t::Vector{Float64}
-    Ysl::Array{Float64, 3}
+    obs::Array{Float64, 3}
-    id::Int64
+    id::String
    YData(a, b, c, d) = new(a, b, c, d)
 end
@@ -214,4 +214,14 @@ end
 function Base.show(io::IO, a::CData)
    print(io, a.header)
+end
+function Base.show(io::IO, a::Corr)
+    fnames = fieldnames(Corr)
+    for k in fnames
+        f = getfield(a, k)
+        if k != :obs
+            print(io, "$k = $f\t")
+        end
+    end
 end
\ No newline at end of file