read_md added

read_rw included. read_rw reads ms1.dat and pbp.dat

src/juobs.jl

@@ -7,7 +7,7 @@ include("juobs_reader.jl")
 include("juobs_tools.jl")
 include("juobs_obs.jl")
 
-export read_mesons, read_ms1, read_ms
+export read_mesons, read_ms1, read_ms, read_md
 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

src/juobs_reader.jl

@@ -156,21 +156,7 @@ function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union
     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")
+function read_rw(path::String; v::String="1.2")
     data = open(path, "r")
     nrw = read(data, Int32)
 
@@ -208,12 +194,51 @@ function read_ms1(path::String; v::String="1.2")
             end
         end
     end
+    close(data)
+    return r_data
+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")
+    r_data = read_rw(path, v=v)
+    nrw = length(r_data)
+    ncnfg = size(r_data[1])[3]
     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_md(path::String)
+
+Reads openQCD  pbp.dat files at a given path. This method returns a matrix md[irw, icfg] that contains the derivatives dS/dm, where
+md[irw=1] = dS/dm_l and md[irw=2] = dS/dm_s
+
+Examples:
+```@example    
+md = read_md(path)
+```
+"""
+function read_md(path::String)
+    r_data = read_rw(path, v="1.4")
+    nrw = length(r_data)
+    ncnfg = size(r_data[1])[3]
+    md = zeros(Float64, nrw, ncnfg)
+    [md[k, :] = prod(mean(r_data[k], dims=2), dims=1) for k = 1:nrw]
+    return md
+end
 
 @doc raw"""
 read_ms(path::String)