Updated dec_const to use plot_data and plot_func, added optional return for...

Updated dec_const to use plot_data and plot_func, added optional return for fig and ax objects, updated documentation

src/juobs_obs.jl

@@ -17,7 +17,7 @@ In case of Vector{Corr} it is better to not broadcast, especially in case of `da
      When `savepl = true`, the plots are saved and closed. When `pl = true` the plot are generated, displayed and the last two returns are always `fig,ax` object to allow more costumability.
      If both are `true`, `pl` is ignored.
     
-     See also: [`plot_data`](@ref), [`plot_funct`](@ref).
+     See also: [`plot_data`](@ref), [`plot_func`](@ref).
   - `data::Bool=false`: if set to `true`, the second return is a vector with the effective mass.
   - `filename::String=""`: if given, overrides the automatically generated filename. `mkpath(dirname(filename))` is called before saving. If not given, a folder "meff/" is made and the correlator `μ` or `κ` are used to name the figure 
 
@@ -115,7 +115,7 @@ function meff(corr::Vector{Corr},plat::Vector{Vector{Int64}}; kwargs...)
   end
 
   # Due to the broadcast, a is a Vector{Tuple{...}} where each element of the vector is the output of the function 
-  # with one set of a0p, pp and plat. The following lines of code reshape a into a Matrix such that each column 
+  # The following lines of code reshape a into a Matrix such that each column 
   # correspond to one Tuple and each row correspond to one component of the Tuple. The it will return each row 
   # separately.
   # 
@@ -152,7 +152,7 @@ The return structure is: mpcac[,data][,figs,axs]
      When `savepl = true`, the plots are saved and closed. When `pl = true` the plot are generated, displayed and the last two returns are always `fig,ax` object to allow more costumability 
      If both are `true`, `pl` is ignored.
     
-     See also: [`plot_data`](@ref), [`plot_funct`](@ref).
+     See also: [`plot_data`](@ref), [`plot_func`](@ref).
   - `data::Bool=false`: if set to `true`, the second return is a vector with the effective PCAC mass.
   - `ca::Float64=0.0`: if given, computes the PCAC mass with the improved axial current. 
   - `filename::String=""`: if given, overrides the automatically generated filename. `mkpath(dirname(filename))` is called before saving. If not given, a folder "mpcac/" is made and the correlator `μ` or `κ` are used to name the figure 
@@ -218,14 +218,23 @@ function mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca:
     fig,ax = plot_data(aux,x,wpm=wpm,ylabel =L"$m_\mathrm{PCAC}$",xlabel=L"$x_0$";plot_kw...)
     isnothing(wpm) ? uwerr(mass) : uwerr(mass,wpm)
     fig,ax = plot_func(x->mass,[plat[1]:plat[2]...],figs=(fig,ax);plot_kw...)
-    if !isnothing(kappa)
-      fig.suptitle(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
+    
+    
+    if !haskey(plot_kw,:title)
+      temp = isnothing(kappa) ? L"$\mu_1 = %$(mu[1])$ $\mu_2 = %$(mu[2])$" : L"$\kappa_1 = %$(kappa[1])$ $\kappa_2 = %$(kappa[2])$"
+      fig.suptitle(temp);
     end
-    if !isnothing(mu)
-      fig.suptitle(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2]))
-    end   
+
+    if !haskey(plot_kw,:ylabel)
+      ax.set_ylabel(L"$m_\mathrm{PCAC}$")
+    end
+    if !haskey(plot_kw,:xlabel) 
+      ax.set_xlabel(L"$x_0$")
+    end
+
     _, max_idx = findmax(value.(pp))
     ax.set_xlim(left=max_idx)
+    ax.set_ylim(v-20*e,v+20*e)
       
     if savepl
       fig.tight_layout();  
@@ -286,21 +295,41 @@ end
 
 ## Decay constants
 @doc raw"""
-    dec_const(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+    dec_const(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=false, data::Bool=false,
+    kappa::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,savepl::Bool=false,
+    filename::String="",plot_kw...)
 
-    dec_const(a0p::Corr, pp::Corr, plat::Vector{Int64}, m::uwreal; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+    dec_const(a0p::Corr, pp::Corr, plat::Vector{Int64}, m::uwreal; kwargs...) 
 
-    dec_const(a0pL::Vector{uwreal}, a0pR::Vector{uwreal}, ppL::Vector{uwreal}, ppR::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+    dec_const(a0p::Vector{Corr},pp::Vector{Corr},plt::Vector{Vector{Int64}},m::Vector{uwreal}; kwargs...)
 
-    dec_const(a0pL::Corr, a0pR::Corr, ppL::Corr, ppR::Corr, plat::Vector{Int64}, m::uwreal; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+    dec_const(a0pL::Vector{uwreal}, a0pR::Vector{uwreal}, ppL::Vector{uwreal}, ppR::Vector{uwreal}, plat::Vector{Int64}, 
+    m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=false, data::Bool=false, kappa::Union{Vector{Float64}, Nothing}=nothing, 
+    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, 
+    savepl::Bool=false, filename::String="",plot_kw...)
 
+    dec_const(a0pL::Corr, a0pR::Corr, ppL::Corr, ppR::Corr, plat::Vector{Int64}, m::uwreal; kwargs...)
 
-Computes the bare decay constant using ``A_0P`` and ``PP`` correlators . 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, effective mass `m` and source position `y0`
-must be specified.
+    dec_const(a0pL::Vector{Corr},a0pR::Vector{Corr},ppL::Vector{Corr}, ppR::Vector{Corr}, plat::Vector{Vector{Int64}}, m::Vector{uwreal}; kwargs...)
+
+    dec_const(vv::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; pl::Bool=false, data::Bool=false,
+    kappa::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
+    savepl::Bool=false,filename::String="",plot_kw...)  
+
+    dec_const(vv::Corr, plat::Vector{Int64}, m::uwreal;kwargs...) 
+    
+    dec_const(vv::Vector{Corr}, plat::Vector{Vector{Int64}},m::Vector{uwreal}; kwargs...)
+
+Computes the bare decay constant using `A_0P` and `PP` or the `VV` correlators and the associated masses `m`. 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, source position `y0` must be specified.
 
-The flags `pl` and `data` allow to show the plots and return data as an extra result. The `ca` variable allows to compute `dec_const` using
-the improved axial current.
+# Arguments
+  - `pl::Bool=false` and `savepl::Bool=false`: show and save the plots. The plots are made through the functions `plot_data` and `plot_func`. It is possible to pass keyword arguments for these two functions. 
+    When `savepl = true`, the plots are saved and closed. When `pl = true` the plot are generated, displayed and the last two returns are always `fig,ax` objects to allow more costumability 
+    If both are `true`, `pl` is ignored.
+  - `data::Bool=false`: if set to `true`, the second return is a vector with the ratio `R`(see below).
+  - `ca::Float64=0.0`: (ONLY IF `A_0P`) if given, computes the decay constant mass with the improved axial current. 
+  - `filename::String=""`: if given, overrides the automatically generated filename. `mkpath(dirname(filename))` is called before saving. If not given, a folder "dec\_cont/" or "dec\_const\_vec" is made and the correlator's `κ` is used to name the figure    
 
 **The method assumes that the source is close to the boundary.** It takes the following ratio to cancel boundary effects.
 ``R = \frac{f_A(x_0, y_0)}{\sqrt{f_P(T-y_0, y_0)}} * e^{m (x_0 - T/2)}``
@@ -308,6 +337,10 @@ the improved axial current.
 **If left and right correlators are included in the input. The result is computed with the following ratio**
 `` R = \sqrt{f_A(x_0, y_0) * f_A(x_0, T - 1 - y_0) / f_P(T - 1 - y_0, y_0)}``
 
+**If the vectorial correlator is used, ther results is computed with the following ratio**
+`` R = \sqrt{\frac{2 f_V(x_0,y_0) * e^{m*(x0-y0)}{m}}}``
+
+
 ```@example
 data_pp = read_mesons(path, "G5", "G5", legacy=true)
 data_a0p = read_mesons(path, "G5", "G0G5", legacy=true)
@@ -331,8 +364,9 @@ f = dec_const(corr_a0p[1], corr_pp[1], [50, 60], m, pl=true, ca=ca)
 f_ratio = dec_const(corr_a0pL, corr_a0pR, corr_ppL, corr_ppR, [50, 60], m, pl=true, ca=ca)
 ```
 """
-function dec_const(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false,
-    kappa::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,savepl::Bool=false)     
+function dec_const(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=false, data::Bool=false,
+    kappa::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,savepl::Bool=false,
+    filename::String="",plot_kw...)     
     
     corr_a0p = -a0p[2:end-1]
     corr_pp = pp[2:end-1]
@@ -354,26 +388,34 @@ function dec_const(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64},
     if pl || savepl
       isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm)
       isnothing(wpm) ? uwerr(f) : uwerr(f, wpm)
-      x = 1:length(R)
-      y = value.(R)
-      dy = err.(R)
-      v = value(R_av)
-      e = err(R_av)
 
-      fig,ax = subplots(1)
+      x = [1:length(R)...]
       lbl = string(L"$af = $", sprint(show, f))
-      ax.fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75, label=L"$R$")
-      ax.errorbar(x, y, dy, fmt="x", color="black", label=lbl)
-      ax.legend()
+      fig,ax = plot_data(R,x,label = lbl;plot_kw...)
+      
+      fig,ax = plot_func(x->R_av,x[plat[1]:plat[2]],figs=(fig,ax),label = L"$R$";plot_kw...)
+      v,e = value(R_av),err(R_av)
+
       ax.set_ylim(v-10*e, v+10*e)
-      ax.set_ylabel(L"$R_\mathrm{av}$")
-      ax.set_xlabel(L"$x_0$")
-      filename = "dec_const/"
-      if !isnothing(kappa)
-          fig.suptitle(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
-          filename*=string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2],".png")
+
+      if !haskey(plot_kw,:ylabel)
+        ax.set_ylabel(L"$R_\mathrm{av}$")
+      end
+      
+      if !haskey(plot_kw,:xlabel)
+        ax.set_xlabel(L"$x_0$")
       end
+
+      if !haskey(plot_kw,:title)
+        temp = isnothing(kappa) ? "" : L"$\kappa_1 = %$(kappa[1])$ $\kappa_2 = %$(kappa[2])$"
+        fig.suptitle(temp);
+      end
+
       if savepl 
+        if filename==""
+          filename = "dec_const/"
+          filename *= isnothing(kappa) ? "" : "kappa_1_$(kappa[1])_kappa_2_$(kappa[2]).png"
+        end
         fig.tight_layout();
         mkpath(dirname(filename));
         fig.savefig(filename);
@@ -382,6 +424,7 @@ function dec_const(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64},
         display(fig)
       end
     end
+
     if pl &&!savepl
       return data ? (f,R,fig,ax) : (f,fig,ax)
     else
@@ -392,7 +435,7 @@ end
 function dec_const(a0p::Corr, pp::Corr, plat::Vector{Int64}, m::uwreal; kwargs...) 
     
     if (a0p.y0 == pp.y0) && (a0p.kappa == pp.kappa)
-        return dec_const(a0p.obs, pp.obs, plat, m, a0p.y0; kwargs...)
+        return dec_const(a0p.obs, pp.obs, plat, m, a0p.y0,kappa=a0p.kappa; kwargs...)
     else 
         error("y0 or kappa values does not match")
     end
@@ -404,7 +447,7 @@ function dec_const(a0p::Vector{Corr},pp::Vector{Corr},plt::Vector{Vector{Int64}}
     return a;
   end
   # Due to the broadcast, a is a Vector{Tuple{...}} where each element of the vector is the output of the function 
-  # with one set of a0p, pp and plat. The following lines of code reshape a into a Matrix such that each column 
+  # . The following lines of code reshape a into a Matrix such that each column 
   # correspond to one Tuple and each row correspond to one component of the Tuple. The it will return each row 
   # separately.
   # 
@@ -412,15 +455,16 @@ function dec_const(a0p::Vector{Corr},pp::Vector{Corr},plt::Vector{Vector{Int64}}
   # V = [(a1,a2,a3),(b1,b2,b3)]
   # V = collect.(V) #[[a1,a2,a3],[b1,b2,b3]]
   # V stack(V) #[a1 b1; a2 b2; a3 b3]
-  # return Tuple([[r...] for e in eachrow(r)]) #([a1,b1],[a1,b2],[a3,b3])
+  # return Tuple([r...] for e in eachrow(r)) #([a1,b1],[a1,b2],[a3,b3])
 
   a = stack(collect.(a));
   return Tuple([r...] for r in eachrow(a))
 end 
 
 function dec_const(a0pL::Vector{uwreal}, a0pR::Vector{uwreal}, ppL::Vector{uwreal}, ppR::Vector{uwreal}, plat::Vector{Int64}, 
-    m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false, kappa::Union{Vector{Float64}, Nothing}=nothing, 
-    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+    m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=false, data::Bool=false, kappa::Union{Vector{Float64}, Nothing}=nothing, 
+    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, 
+    savepl::Bool=false, filename::String="",plot_kw...)
 
     corr_pp = (ppL + ppR[end:-1:1]) / 2
     T = length(corr_pp)
@@ -441,48 +485,86 @@ function dec_const(a0pL::Vector{uwreal}, a0pR::Vector{uwreal}, ppL::Vector{uwrea
     R_av = plat_av(R, plat, wpm)
 
     f = sqrt(2) * sqrt(R_av^2) / sqrt(m)
-    if pl
-        isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm)
-        isnothing(wpm) ? uwerr(f) : uwerr(f, wpm)
-        x = 1:length(R)
-        y = value.(R)
-        dy = err.(R)
-        v = value(R_av)
-        e = err(R_av)
+    if pl || savepl
+      isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm)
+      isnothing(wpm) ? uwerr(f) : uwerr(f, wpm)
 
-        figure()
-        lbl = string(L"$af = $", sprint(show, f))
-        fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75, label=L"$R$")
-        errorbar(x, y, dy, fmt="x", color="black", label=lbl)
-        legend()
-        ylabel(L"$R_\mathrm{av}$")
-        xlabel(L"$x_0$")
-    
-        if !isnothing(kappa)
-            title(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
+      x = [1:length(R)...]
+      lbl = string(L"$af = $", sprint(show, f))
+      fig,ax = plot_data(R,x,label = lbl;plot_kw...)
+      
+      fig,ax = plot_func(x->R_av,x[plat[1],plat[2]],figs=(fig,ax),label = L"$R$";plot_kw...)
+      v,e = value(R_av),err(R_av)
+
+      ax.set_ylim(v-10*e, v+10*e)
+
+      if !haskey(plot_kw,:ylabel)
+        ax.set_ylabel(L"$R_\mathrm{av}$")
+      end
+      
+      if !haskey(plot_kw,:xlabel)
+        ax.set_xlabel(L"$x_0$")
+      end
+
+      if !haskey(plot_kw,:title)
+        temp = isnothing(kappa) ? L"$\mu_1 = %$(mu[1])$ $\mu_2 = %$(mu[2])$" : L"$\kappa_1 = %$(kappa[1])$ $\kappa_2 = %$(kappa[2])$"
+        fig.suptitle(temp);
+      end
+
+      if savepl 
+        if filename==""
+          filename = "dec_const/"
+          filename *= isnothing(kappa) ? "" : "kappa_1_$(kappa[1])_kappa_2_$(kappa[2]).png"
         end
-        display(gcf())
+        fig.tight_layout();
+        mkpath(dirname(filename));
+        fig.savefig(filename);
+        close(fig)
+      elseif pl
+        display(fig)
+      end
     end
-    if !data
-        return f
+
+    if pl&&!savepl
+      return data ? (f,R,fig,ax) : (f,fig,ax)
     else
-        return (f, R)
+        return data ? (f, R) : f
     end
 end
 
-function dec_const(a0pL::Corr, a0pR::Corr, ppL::Corr, ppR::Corr, plat::Vector{Int64}, m::uwreal; ca::Float64=0.0, pl::Bool=true, data::Bool=false, 
-    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+function dec_const(a0pL::Corr, a0pR::Corr, ppL::Corr, ppR::Corr, plat::Vector{Int64}, m::uwreal; kwargs...)
 
     T = length(a0pL.obs)
     if (a0pL.y0 == ppL.y0) && (a0pR.y0 == ppR.y0) && (a0pL.kappa == ppL.kappa) && (a0pR.kappa == ppR.kappa) && (a0pL.y0 == T - 1 - a0pR.y0)
-        return dec_const(a0pL.obs, a0pR.obs, ppL.obs, ppR.obs, plat, m, a0pL.y0, ca=ca, kappa=a0pL.kappa, pl=pl, data=data, wpm=wpm)
+        return dec_const(a0pL.obs, a0pR.obs, ppL.obs, ppR.obs, plat, m, a0pL.y0; kwargs...)
     else 
         error("y0 or kappa values does not match")
     end
 end
-## ADD DOCUMENTATION FOR VECTOR DECAY 
-function dec_const(vv::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; pl::Bool=true, data::Bool=false,
-    kappa::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)     
+
+function dec_const(a0pL::Vector{Corr},a0pR::Vector{Corr},ppL::Vector{Corr}, ppR::Vector{Corr}, plat::Vector{Vector{Int64}}, m::Vector{uwreal}; kwargs...)
+  a = dec_const.(a0pL,a0pR,ppL,ppR,plat,m;kwargs...)
+  if length(a[1])==1
+    return a;
+  end
+  # Due to the broadcast, a is a Vector{Tuple{...}} where each element of the vector is the output of the function 
+  # The following lines of code reshape a into a Matrix such that each column 
+  # correspond to one Tuple and each row correspond to one component of the Tuple. The it will return each row 
+  # separately.
+  # 
+  # Example:
+  # V = [(a1,a2,a3),(b1,b2,b3)]
+  # V = collect.(V) #[[a1,a2,a3],[b1,b2,b3]]
+  # V stack(V) #[a1 b1; a2 b2; a3 b3]
+  # return Tuple([r...] for e in eachrow(r)) #([a1,b1],[a1,b2],[a3,b3])
+
+  a = stack(collect.(a));
+  return Tuple([r...] for r in eachrow(a))
+end
+
+function dec_const(vv::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; pl::Bool=false, data::Bool=false,
+    kappa::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
+    savepl::Bool=false,filename::String="",plot_kw...)     
     
     corr_vv = vv[2:end-1]
     T = length(corr_vv)
@@ -493,46 +575,78 @@ function dec_const(vv::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64
     R_av = plat_av(R, plat, wpm)
     f = sqrt(2 / m)  * R_av
     R .*= sqrt.(2 ./ [m for i in 1:length(R)])
-    if pl
+    if pl || savepl
         uwerr.(R)
         R_av *= sqrt(2/m)
         isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm)
         isnothing(wpm) ? uwerr(f) : uwerr(f, wpm)
-        x = 1:length(R)
-        y = value.(R)
-        dy = err.(R)
-        v = value(R_av)
-        e = err(R_av)
-
-        figure()
+        x = [1:length(R)...]
         lbl = string(L"$af = $", sprint(show, f))
-        fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75, label=L"$R$")
-        errorbar(x, y, dy, fmt="x", color="black", label=lbl)
-        legend()
-        xlim(left=y0)
-        ylim(v-10*e, v+10*e)
-        ylabel(L"$af_{B}$")
-        #ylabel(L"$R_\mathrm{av}$")
-        xlabel(L"$x_0$")
-        #title(L"$f_{B^*}$")
-        if !isnothing(kappa)
-            title(string( L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
+        fig,ax = plot_data(R,x,label=lbl;plot_kw...)
+        fig,ax = plot_func(x->R_av,x[plat[1]:plat[2]],label = L"$R$"; plot_kw...)
+        v,e = value(R_av),err(R_av)
+
+        ax.set_xlim(left=y0)
+        ax.set_ylim(v-10*e, v+10*e)
+        
+        if !haskey(plot_kw,:ylabel)
+          ylabel(L"$R_\mathrm{av}$")
         end
-        display(gcf())
-        #t = "ps_decay_$(kappa[1])_$(kappa[2]).pdf"
-        #savefig(joinpath("/Users/ale/Il mio Drive/phd/secondment/3pf test/analysis/plots",t))
+        
+        if !haskey(plot_kw,:xlabel)
+          ax.xlabel(L"$x_0$")
+        end
+        
+        if !haskey(plot_kw,:title)
+          if !isnothing(kappa)
+            fig.suptitle(string( L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
+          end
+          ax.title(L"$f_{B^*}$")
+        end
+
+        if savepl
+          if filename==""
+            filename = "dec_const_vec/"
+            filename *= isnothing(kappa) ? "f_B_vec.png" : string( L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2],".png")
+          end
+          mkpath(filename)
+          fig.tight_layout()
+          fig.savepl(filename)
+          close(fig)
+        else
+          display(gcf())
+        end
+
     end
-    if !data
-        return f
+    if pl&&!savepl
+      return data ? (f,R,fig,ax) : (f,fig,ax)
     else
-        return f, R
+      return data ? (f,R) : f
     end
 end
-function dec_const(vv::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false,
-    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) 
-    
-    return dec_const(vv.obs, plat, m, vv.y0, kappa=vv.kappa, pl=pl, data=data, wpm=wpm)
+function dec_const(vv::Corr, plat::Vector{Int64}, m::uwreal;kwargs...) 
+    return dec_const(vv.obs, plat, m, vv.y0, kappa=vv.kappa; kwargs...)
 end
+function dec_const(vv::Vector{Corr}, plat::Vector{Vector{Int64}},m::Vector{uwreal}; kwargs...)
+  a= dec_const.(vv,plat,m;kwargs...)
+  if length(a[1])==1
+    return a;
+  end
+  # Due to the broadcast, a is a Vector{Tuple{...}} where each element of the vector is the output of the function 
+  # The following lines of code reshape a into a Matrix such that each column 
+  # correspond to one Tuple and each row correspond to one component of the Tuple. The it will return each row 
+  # separately.
+  # 
+  # Example:
+  # V = [(a1,a2,a3),(b1,b2,b3)]
+  # V = collect.(V) #[[a1,a2,a3],[b1,b2,b3]]
+  # V stack(V) #[a1 b1; a2 b2; a3 b3]
+  # return Tuple([r...] for e in eachrow(r)) #([a1,b1],[a1,b2],[a3,b3])
+
+  a = stack(collect.(a));
+  return Tuple([r...] for r in eachrow(a))
+end
+
 # this can be used also for pseudoscalar decay constants with wilson fermions 
 function dec_const_w(vv::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false,
     kappa::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)     

src/juobs_plots.jl

@@ -86,7 +86,7 @@ end
       wpm::Union{Dict{String,Vector{Float64}},Nothing} = nothing,
       c::String ="green",
       alpha::Float64=0.25,
-      label::String = "",
+      label::AbstractString = "",
       figs::Union{Tuple{Figure, PyPlot.PyCall.PyObject},Nothing}=nothing,
 It plots the function `f` evaluated at `xdata`through broadcasting. 
 The function assumes that `f` return type is `uwreal` and is called as `f(x)`.
@@ -115,7 +115,7 @@ function plot_func(f::Function, xdata::Vector{T} where T<:Real;
   c::String ="green",
   alpha::Float64=0.25,
   figs::Union{Tuple{Figure, PyPlot.PyCall.PyObject},Nothing}=nothing,
-  label::String = "",
+  label::AbstractString = "",
   kwargs...) 
   fig,ax = isnothing(figs) ? subplots(1) : figs
   y = f.(xdata);