small bug fixed, test with t-tnot fixed in gevp, analysis added for all the ensembles

gevp tools/results.rtf

@@ -5,7 +5,23 @@
 \paperw11900\paperh16840\margl1440\margr1440\vieww10800\viewh8400\viewkind0
 \pard\tx566\tx1133\tx1700\tx2267\tx2834\tx3401\tx3968\tx4535\tx5102\tx5669\tx6236\tx6803\pardirnatural\partightenfactor0
 
-\f0\fs24 \cf0 $N200\
+\f0\fs24 \cf0 $H400\
+$Z^\{tm\}_M \\mu_c \\sqrt\{8t_0\} =  $3.030627780344985 +/- 0.07902594921440798\
+$M_D \\sqrt\{8t_0\} =  $3.981135432856687 +/- 0.05781090729513162\
+$M_Ds \\sqrt\{8t_0\} =  $3.981135432856687 +/- 0.05781090729513162\
+$M_D* \\sqrt\{8t_0\} =  $4.278807884262479 +/- 0.05002067830312485\
+$M_Ds* \\sqrt\{8t_0\} =  $4.278807884262479 +/- 0.05002067830312485\
+$M_eta* \\sqrt\{8t_0\} =  $6.096944664722402 +/- 0.09777085253831523\
+$M_jpsi* \\sqrt\{8t_0\} =  $6.3445407150455315 +/- 0.09568980579456061\
+$f_D \\sqrt\{8t_0\} =  $0.516849796654617 +/- 0.0030631272287505134\
+$f_Ds \\sqrt\{8t_0\} =  $0.516849796654617 +/- 0.0030631272287505134\
+$f_D* \\sqrt\{8t_0\} =  $0.21312484071325435 +/- 0.003968470514564074\
+$f_Ds* \\sqrt\{8t_0\} =  $0.21312484071325435 +/- 0.003968470514564074\
+$f_eta \\sqrt\{8t_0\} =  $0.9689427427736942 +/- 0.015398690942284236\
+$f_jpsi \\sqrt\{8t_0\} =  $0.4391703592723876 +/- 0.009322107275879348\
+$\
+\pard\tx566\tx1133\tx1700\tx2267\tx2834\tx3401\tx3968\tx4535\tx5102\tx5669\tx6236\tx6803\pardirnatural\partightenfactor0
+\cf0 $N200\
 $Z^\{tm\}_M \\mu_c \\sqrt\{8t_0\} =  $3.095508084190539 +/- 0.07424875770953374\
 $M_D \\sqrt\{8t_0\} =  $3.958264727738654 +/- 0.05515406625525418\
 $M_Ds \\sqrt\{8t_0\} =  $4.072538504215807 +/- 0.05467443552294891\
@@ -35,6 +51,21 @@ $f_Ds* \\sqrt\{8t_0\} =  $0.20485862913653624 +/- 0.0024571277939203565\
 $f_eta \\sqrt\{8t_0\} =  $0.9262934445841321 +/- 0.011380030137209169\
 $f_jpsi \\sqrt\{8t_0\} =  $0.4287220496155749 +/- 0.007707612363034976\
 $\
+$N300\
+$Z^\{tm\}_M \\mu_c \\sqrt\{8t_0\} =  $2.996157019535518 +/- 0.07903034846063876\
+$M_D \\sqrt\{8t_0\} =  $3.9418982012432258 +/- 0.0596919625451062\
+$M_Ds \\sqrt\{8t_0\} =  $3.9418982012432258 +/- 0.0596919625451062\
+$M_D* \\sqrt\{8t_0\} =  $4.2918918432116255 +/- 0.04985913651457226\
+$M_Ds* \\sqrt\{8t_0\} =  $4.2918918432116255 +/- 0.04985913651457226\
+$M_eta* \\sqrt\{8t_0\} =  $6.077559366401819 +/- 0.10543357205926039\
+$M_jpsi* \\sqrt\{8t_0\} =  $6.322315861178219 +/- 0.1024859602694168\
+$f_D \\sqrt\{8t_0\} =  $0.502826744878361 +/- 0.0026201800440417583\
+$f_Ds \\sqrt\{8t_0\} =  $0.502826744878361 +/- 0.0026201800440417583\
+$f_D* \\sqrt\{8t_0\} =  $0.2068623601937379 +/- 0.003546907525867658\
+$f_Ds* \\sqrt\{8t_0\} =  $0.2068623601937379 +/- 0.003546907525867658\
+$f_eta \\sqrt\{8t_0\} =  $0.8669858354378843 +/- 0.010734940135830663\
+$f_jpsi \\sqrt\{8t_0\} =  $0.402925306634495 +/- 0.008184816863154212\
+$\
 $J303\
 $Z^\{tm\}_M \\mu_c \\sqrt\{8t_0\} =  $3.106805094851814 +/- 0.07428519131808012\
 $M_D \\sqrt\{8t_0\} =  $3.951817203146415 +/- 0.05634686437630071\

gevp tools/tools.jl

@@ -95,6 +95,23 @@ mutable struct  infoMatt
         mat = get_matrix(diag, subdiag)
         return new(mat, mu, y0)
     end
+    function infoMatt(a11::Array{juobs.Corr}, a12::Array{juobs.Corr}, a13::Array{juobs.Corr}, a22::Array{juobs.Corr}, a23::Array{juobs.Corr}, a33::Array{juobs.Corr})
+        mu = getfield(a11[1], :mu)
+        y0 = Int64(getfield(a11[1], :y0))
+        elem11 = sum(a11[i].obs for i =1:length(a11))
+        elem12 = sum(a12[i].obs for i =1:length(a12))
+        elem22 = sum(a22[i].obs for i =1:length(a22))
+        elem13 = sum(a13[i].obs for i =1:length(a13))
+        elem23 = sum(a23[i].obs for i =1:length(a23))
+        elem33 = sum(a33[i].obs for i =1:length(a33))
+        diag = [elem11, elem22, elem33]
+        subdiag =[ elem12, elem13, elem23 ]
+        #diag = [a11.obs, a22[1].obs.+a22[2].obs.+a22[3].obs]
+        #subdiag = [a12[1].obs.+a12[2].obs.+a12[3].obs]
+        mat = get_matrix(diag, subdiag)
+        return new(mat, mu, y0)
+    end
+    
 end
 function comp_mat(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}})
     mu = getfield.(tot11[1],:mu)
@@ -108,6 +125,21 @@ function comp_mat(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr
     end
     return res
 end
+function comp_mat(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot13::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}}, tot33::Array{Array{juobs.Corr}}, tot23::Array{Array{juobs.Corr}})
+    mu = getfield.(tot11[1],:mu)
+    res = Array{infoMatt}(undef, length(mu))
+    for i = 1:length(mu)
+        a11 = [ tot11[j][i] for j = 1:size(tot11,1) ]
+        a12 = [ tot12[j][i] for j = 1:size(tot12,1) ]
+        a22 = [ tot22[j][i] for j = 1:size(tot22,1) ]
+        a33 = [ tot33[j][i] for j = 1:size(tot33,1) ]
+        a13 = [ tot13[j][i] for j = 1:size(tot13,1) ]
+        a23 = [ tot23[j][i] for j = 1:size(tot23,1) ]
+        res[i] = infoMatt(a11, a12, a13, a22, a23, a33)
+        #res[i]=  infoMatt(tot11[i], [tot12[1][i],tot12[2][i],tot12[3][i]], [tot22[1][i],tot22[2][i],tot22[3][i]])
+    end
+    return res
+end
 function comp_energy(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}}, evec::Bool=false)
     aux_mat = comp_mat(tot11, tot12, tot22)
     y0 = getfield(aux_mat[1], :y0)
@@ -124,6 +156,30 @@ function comp_energy(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.C
         for i = 1:length(aux_mat)
             mat_obs = getfield(aux_mat[i], :mat)[y0:end-1]
             evalues, eigvec = uwgevp_tot(mat_obs, 5, evec=true)
+            println("in comp_energy, i is = ", i)
+            res_en[i] = infoEn(energies(evalues), aux_mat[i])
+            evec[i] = eigvec
+        end
+        return res_en, evec
+    end
+end
+function comp_energy(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot13::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}}, tot33::Array{Array{juobs.Corr}}, tot23::Array{Array{juobs.Corr}}, evec::Bool=false; tnot::Int64=5)
+    aux_mat = comp_mat(tot11, tot12, tot13, tot22, tot23, tot33)
+    y0 = getfield(aux_mat[1], :y0)
+    res_en = Array{infoEn}(undef, length(aux_mat))
+    if !evec
+        for i = 1:length(aux_mat)
+            mat_obs = getfield(aux_mat[i], :mat)[y0:end-1]
+            evalues = uwgevp_tot(mat_obs, tnot)
+            res_en[i] = infoEn(energies(evalues), aux_mat[i])
+        end
+        return res_en
+    else
+        evec = Array{Array{}}(undef, length(aux_mat))
+        for i = 1:length(aux_mat)
+            mat_obs = getfield(aux_mat[i], :mat)[y0:end-1]
+            evalues, eigvec = uwgevp_tot(mat_obs, tnot, evec=true)
+            println("in comp_energy, i is = ", i)
             res_en[i] = infoEn(energies(evalues), aux_mat[i])
             evec[i] = eigvec
         end
@@ -142,15 +198,15 @@ mutable struct infoEn
     end
 end
 function pseudo_mat_elem(evec::Array{Array{T,2} where T,1}, ppcorr::Array{uwreal,1}, mass::uwreal)
-    Rn = [ (uwdot(evec[i][:,1], uwdot(ppcorr[97+i], evec[i][:,1])))^(0.5) * exp(mass*i/2)  for i=1:length(evec)-2]
+    Rn = [ (uwdot(evec[i][:,1], uwdot(ppcorr[28+i], evec[i][:,1])))^(0.5) * exp(mass*i/2)  for i=1:length(evec)-2]
     return Rn
 end
 function extract_dec_const(mat_elem::uwreal, mass::uwreal, mu::Array{Float64})
     return sqrt(2)*sum(mu)*mat_elem / mass^1.5
 end
 function vec_mat_elem(evec::Array{Array{T,2} where T,1}, akakcorr::Array{uwreal,1},gevpmat::Array{Array{T,2} where T,1},  mass::uwreal )
-    Rn = [ (uwdot(evec[i][:,2], uwdot(akakcorr[97+i]/3, evec[i][:,2]))^2)^(-0.25) * exp(mass*i/2)  for i=1:length(evec)-2]
-    aux = vec([ uwdot(evec[i][:,2], gevpmat[97+i][:,2]/3) for i =1:length(evec)-2])
+    Rn = [ (uwdot(evec[i][:,2], uwdot(akakcorr[28+i], evec[i][:,2]))^2)^(-0.25) * exp(mass*i/2)  for i=1:length(evec)-2]
+    aux = vec([ uwdot(evec[i][:,2], gevpmat[28+i][:,2]) for i =1:length(evec)-2])
     return Rn.*aux
 end
 function match_muc(muh, m_lh, m_lh_star, target)

src/juobs_linalg.jl

@@ -117,20 +117,26 @@ The columns of such matrix are the eigenvectors associated with eigenvalues eval
 The keyword iter, set by default to 30, selects the number of iterations of the qr algorithm before stopping.
 """
 function uwgevp_tot(mat_list::Vector{Matrix}, tnot::Int64; iter::Int64 = 30, evec::Bool = false)
-    c_inv = invert(mat_list[tnot])
     if !evec
+        c_inv = invert(mat_list[tnot])
         aux_evals = uwgevp.(mat_list, c_tnot_inv=c_inv, iter = iter)
         evals = get_diag.(aux_evals)
         return evals
     else
         n = length(mat_list)
-        evals = Array{Array{uwreal}}(undef, n )
-        evecs = Array{Matrix}(undef, n)
-        for i = 1:n
+        println(n)
+        evals = Array{Array{uwreal}}(undef, 0 )
+        evecs = Array{Matrix}(undef, 0)
+        for i = 6:n-1
+            println("The time t is =", i)
+            #tnot+=1
+            c_inv = invert(mat_list[tnot])
+            println("The time t0 is =", tnot)
             aux_evals, aux_evec =  uwgevp(mat_list[i], c_tnot_inv=c_inv, iter = iter, evec = evec)
-            evals[i] =  get_diag(aux_evals)
-            evecs[i] = aux_evec
+            push!(evals,  get_diag(aux_evals))
+            push!(evecs, aux_evec)
         end
+        println("In uwgevp_tot, the time extent is ", length(evals))
         return evals, evecs
         #=
         (aux_evals, aux_evec) =  uwgevp.(mat_list, c_tnot_inv=c_inv, iter = iter, evec= true)

src/juobs_obs.jl

@@ -16,7 +16,7 @@ m = meff(corr_pp[1], [50, 60], pl=false)
 """
 function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bool=false)     
     dim = length(corr)                                                                                    
-    aux = log.(corr[2:dim-2] ./ corr[3:dim-1])     
+    aux = 0.5.*log.((corr[2:dim-2] ./ corr[3:dim-1]).^2)     
     mass = plat_av(aux, plat)                                                                            
     uwerr(mass)                       
     if pl == true
@@ -69,7 +69,7 @@ function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu
     dim = length(corr_pp)
                                                                  
     aux = exp.((collect(1:dim) .- y0 ) .* [m for k in 1:dim])
-    R = (aux .* corr_pp).^0.5
+    R = ((aux .* corr_pp).^2).^0.25
     R_av = plat_av(R, plat)
     f = sqrt(2) * (mu[1] + mu[2]) *R_av / m^1.5
     uwerr(f)
@@ -88,7 +88,7 @@ function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu
     if data == false
         return f
     else
-        return (f, R)                   
+        return (f, uwdot(sqrt(2) * (mu[1] + mu[2]) / m^1.5, R))                   
     end 
 end 
 dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false) =

src/juobs_tools.jl

@@ -120,10 +120,14 @@ m2_phys = fitp[1] + fitp[2] * phi2_phys
 function lin_fit(x::Vector{Float64}, y::Vector{uwreal})
     uwerr.(y)
     par = lin_fit(x, value.(y), err.(y))
-
     chisq(p, d) = sum((d .- p[1] .- p[2].*x).^2 ./ err.(y) .^2)
     (fitp, csqexp) = fit_error(chisq, par, y)
-    
+    for i = 1:length(fitp)
+        uwerr(fitp[i])
+        print("\n Fit parameter: ", i, ": ")
+        details(fitp[i])
+    end
+    println("Chisq / chiexp: ", chisq(par, y), " / ", csqexp, " (dof: ", x[end]-length(par),")")
     return (fitp, csqexp)
 end