merged with IFT pc version

constants/juobs_const.jl

@@ -112,7 +112,6 @@ The default value of acceptance has no real justification at the moment
 If verbose = true, it prints a warning if the g_0^2 is > 1.8  
 
 """
-
 function ca(beta::Float64; verbose::Bool = true, acceptance::Float64 = 10^-2)
   g02 = 6/beta; 
   if (1.8-g02)>-acceptance # for g_0^2 about 1.8 and below it work. I decided on a 10^-2 acceptance range

src/juobs_obs.jl

@@ -80,15 +80,18 @@ function meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false,
     end
 end
 @doc raw"""
-    mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+    mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,savepl::Bool=false, label::Union{String,Nothing}=nothing)
 
-    mpcac(a0p::Corr, pp::Corr, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+    mpcac(a0p::Corr, pp::Corr, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,savepl::Bool=false, label::Union{String,Nothing}=nothing)
+
+    mpcac(a0p::Vector{Corr},pp::Vector{Corr},plat::Vector{Vector{int}}; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, savepl::Bool=false, label::Union{String,Nothing}=nothing)
 
 Computes the bare PCAC mass for a given correlator `a0p` and `pp` 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. The `ca` variable allows to compute `mpcac` using
-the improved axial current.
+the improved axial current. The flag `savepl` allow to save automatically the plot, the variable `label` if given is used as lalel for the plot 
+and in the filename.
 
 ```@example
 data_pp = read_mesons(path, "G5", "G5")
@@ -102,7 +105,7 @@ p1 = -13.9847
 g2 = 1.73410
 ca = -0.006033 * g2 *( 1 + exp(p0 + p1/g2))
 
-m12 = mpcac(corr_a0p, corr_pp, [50, 60], pl=false, ca=ca)
+m12 = mpcac(corr_a0p, corr_pp, [50, 60], pl=false, ca=ca, savepl=true, label="with_ca_improvement")
 ```
 """
 function mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false, 
@@ -192,6 +195,17 @@ function mpcac(a0p::Corr, pp::Corr, plat::Vector{Int64}; ca::Float64=0.0, pl::Bo
     end
 end
 
+function mpcac(a0p::Vector{Corr},pp::Vector{Corr},plat::Vector{Vector{Int64}}; ca::Float64=0.0, pl::Bool=true, data::Bool=false,
+    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
+    savepl::Bool=false, label::Union{String,Nothing}=nothing)
+    if !data 
+        return mpcac.(a0p,pp,plat,ca=ca,pl=pl,data=data,wpm=wpm,savepl=savepl,label=label)
+    else
+        aux = mpcac.(a0p,pp,plat,ca=ca,pl=pl,data=data,wpm=wpm,savepl=savepl,label=label)
+        return (getfield.(aux,1),getfield.(aux,2))
+    end
+end
+#= 
 @doc raw"""
     mpcac_ren(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, za_zp::uwreal; ca::Float64=0.0,ba_bp_z::uwreal=0.0,pl::Bool=true, data::Bool=false,
             kappa::Union{Vector{Float64}, Nothing}=nothing, mu::Union{Vector{Float64}, Nothing}=nothing, 

src/juobs_tools.jl

@@ -1104,12 +1104,12 @@ function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}
     
     ## CODE HERE THE UPDATED VERSION OF CORRELATED FITS
     if isnothing(inv_cov)
-        @info("Uncorrelated fit")
+        # @info("Uncorrelated fit")
         chisq = gen_chisq(model, xdata, yer)
         fit = curve_fit(model, xdata, yval, 1.0 ./ yer.^2, fill(0.5, param))
         (upar, chi_exp) = isnothing(wpm) ? fit_error(chisq, coef(fit), ydata) : fit_error(chisq, coef(fit), ydata, wpm)
     else
-        @info("Correlated fit")
+        # @info("Correlated fit")
         #covar_inv = inv(get_covariance(ydata)) 
         #covar_inv = (covar_inv + covar_inv')/2
         #covar_inv = inv(Symmetric(get_covariance(ydata)))
@@ -1126,8 +1126,8 @@ function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}
     # end
     # println("\n")
 
-    println("chi2 from fit residual = ", chi2_fit_res)
-    println("chi2 from chi2 function = ", chisq(coef(fit), ydata))
+    # println("chi2 from fit residual = ", chi2_fit_res)
+    # println("chi2 from chi2 function = ", chisq(coef(fit), ydata))
 
 
     #Info
@@ -1142,9 +1142,9 @@ function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}
     chis2_corrected = (length(yval) - param) * chi2_fit_res / chi_exp
 
     #println("Chisq / chiexp: ", chisq(coef(fit), ydata), " / ", chi_exp, " (dof: ", length(yval) - param,")")
-    println("Chisq / chiexp: ", chi2_fit_res, " / ", chi_exp, " (dof: ", length(yval) - param,")")
-    println("Chisq corrected: ", chis2_corrected)
-    println("Return: params, chi2, chi2exp")
+    # println("Chisq / chiexp: ", chi2_fit_res, " / ", chi_exp, " (dof: ", length(yval) - param,")")
+    # println("Chisq corrected: ", chis2_corrected)
+    # println("Return: params, chi2, chi2exp")
     return upar, chi2_fit_res,  chi_exp
 end
 
@@ -1604,6 +1604,8 @@ function pvalue(chisq::Function,
         x = x / eig[1]
         #dQ = juobs.mean((x .> 0) .* exp.(-x * 0.5) * 0.5 ./ sqrt.(abs.(x)))
         #dQ = err(cse)/value(cse) * dQ
+    else
+        error("negative degree of freedom!")
     end
     return Q
 end
@@ -1668,6 +1670,8 @@ function pvalue(chisq::Function,
         x = x / eig[1]
         #dQ = juobs.mean((x .> 0) .* exp.(-x * 0.5) * 0.5 ./ sqrt.(abs.(x)))
         #dQ = err(cse)/value(cse) * dQ
+    else
+        error("negative degree of freedom!")
     end
     return Q
 end