Skip to content

G
GEVP automated analysis
Commit 6c956856 authored by Alessandro 's avatar Alessandro
Browse files
Options

General updated of the analysis code used for the 2022 Lattice Conference

parent 838f8293
Hide whitespace changes
Inline Side-by-side
Showing with 4459 additions and 485 deletions
README.md
View file @ 6c956856
Software for automated data analysis of Wilson_tm fermions with GEVP techniques.
Software for general purpose analysis of charm physics data from a tmQCD mixed action, from mass shift to gevp extraction of the observables to continuum limit fits.
-> src
|----> aux obs computation
| |------> md.jl : main routine for the computation of delta m and t0 shifted
| |------> tools.jl, types.jl, const.jl : auxiliary files to md.jl
|
|-----> gevp analysis: obs. analysis once md.jl has been executed
|-----> gevp_tm.jl : extract observable ensemble by ensemble using gevp and store in BDIO
|-----> gevp_fit.jl : read BDIO from previous file and perform chiral-continumm extrapolations
|-----> func_comb.jl : contains all the functional forms for the chiral-continuum fits
|-----> gevp_plotter.jl : useful routines for plotting complicated chiral-continuum fits and more
|-----> read_bdio.jl : reading BDIO routines ( the order is hardcoded but consistent within the whole files)
|-----> const.jl, tools.jl, types.jl : auxiliary routines for the analysis
The GEVP is solved as outlined in 1010.0202, i.e., given a correlator f(t) the following matrix is constructed:
```math
......@@ -16,57 +32,27 @@ The effective energy spectrum is determined from:
```math
E_n^{eff}(t, t_0, \tau) = \log( \frac{\lambda_n(t, t_0, \tau)}{\lambda_n(t+a, t_0, \tau)}) = E_n + O(e^{-(E_m-e_n)t})
```
Finally the energy levels are determined by fitting the effective energy with the functional form
```math
E_n^{eff}(t, t_0, \tau) = E_n + p_1*e^{-\Delta E * t}
```
The fit is performed in the time extent
[range_t_fit[1]:range_t_fit[2], end-5]
where range_t_fit is an array of integers set by default to [2,5], but this can be changed by the user.
Final results for the energy spectrum is given by the weighthed average of the fit parameters determined by the above outlined procedure.
N.B. at the moment the code does not extract matrix element, hence decay constants, from the gevp.
To run the code in the juliaREPL type:
julia> include("/path_to/gevp.jl")
Go to SET UP VARIABLES section in the gevp.jl file to properly configure before running.
You have to up set the following variables:
- path_data where data are stored
- path_plat where a plat.txt file is stored with plateaus for all ensembles
- path_result where you want to save results.txt and plots
- ensembles to analyse
- sector you are interested in: light-light, light-heavy, heavy-heavy...
- tau value to use for matrix construction
- _t0 value to use for solving the GEVP
- range_t_fit array of integers with the fit range value for the extractions of the energies
- rwf set to true if you want to include reweighting factors
- compute_t0 set to true if you want to extract t0 from ms.dat file. Otherwise, t0 is taken from 1608.08900 in const.jl
- mass_shift set to true if you want to include the mass shift
go to line 484 to modify continuum limit + chiral extrapolatio fit model
Go to SET UP VARIABLES section in md.jl, gevp_tm.jl and gevp_fit.jl files to properly configure before running.
Several plauteau and data paths have to be set before running.
Data in path_data need to be stored as follow:
Data in the various path_data need to be stored as follow:
../data_folder:
path_data_md : contains pbp.dat files for action derivatives
------> path_data_md/ens_i/*pbp.dat
path_data_w: contains wilson light correlators for mass shift and ms.dat files for t0
------> path_data_w/ens_i/*ms.dat
------> path_data_w/ens_i/*mesons.dat
../data_folder/ens_id_1_folder:
- mesons.dat
- ms.dat
- ms1.dat
- pbp.dat
path_data_rw: contains reweighting factors ms1.dat files
-------> path_data_rw/ens_i/*ms1.dat
path_data_tm: contains twisted mass (heavy) correlators
--------> path_data_tm/ens_i/*mesons.dat
../data_folder/ens_id_2_folder:
Replicas are automatically taken into account for each ensembles.
- r0_mesons.dat
- r1_mesons.dat
- r0_ms.dat
- r1_ms.dat
- r0_ms1.dat
- r1_ms1.dat
- r0_pbp.dat
- r1_pbp.dat
N.B.
In the gamma method the number of configurations of different observables but with same ensemble HAS to match. There might be runs where i.e. ms.dat files contains more cnfgs than mesons.dat files and this generates errors. To avoid it, set in const.jl/trunc_db (database for each ensembles) the number of configurations at which you want to cut all the observables with that ensembe id.
and so on
\ No newline at end of file
src/plat_dec.txt → plat/plat_tm_dps.txt
View file @ 6c956856
#H101
ll 58 75
lh 58 70
hh 60 72
lh 63 76
hh 65 74
t0 20 70
#H102r001
ll 60 85
ls 60 85
lh 60 85
ss 60 85
sh 60 85
hh 60 85
lh 62 80
ss 60 85
sh 62 80
hh 65 81
t0 20 80
#H102r002
ll 58 73
ls 60 75
lh 62 80
lh 65 75
ss 60 80
sh 62 80
hh 63 80
sh 65 75
hh 66 78
t0 20 80
#H105
ll 58 73
ls 60 75
lh 63 72
ss 60 80
sh 65 80
hh 66 80
t0 20 80
#H400
ll 60 75
......@@ -50,6 +58,14 @@ ss 85 105
sh 83 110
hh 90 104
t0 20 110
#D200
ll 85 105
ls 82 100
lh 80 88
ss 85 100
sh 82 88
hh 87 93
t0 20 110
#J303
ll 120 140
ls 123 140
......
src/plat_dec_vec.txt → plat/plat_tm_dvec.txt
View file @ 6c956856
#H101
ll 56 70
lh 59 68
ll 60 70
lh 63 72
hh 68 75
t0 20 70
#H102r001
ll 60 80
ls 60 85
lh 57 70
ss 57 67
ss 57 67
sh 60 70
hh 63 79
t0 20 80
#H102r002
ll 60 80
ls 60 85
lh 57 70
ss 57 67
sh 60 70
ll 58 65
ls 58 65
lh 60 67
ss 57 62
sh 60 66
hh 65 79
t0 20 80
#H105
ll 56 62
ls 56 65
lh 60 70
ss 57 63
sh 61 70
hh 69 80
t0 20 80
#H400
ll 56 68
lh 60 73
hh 65 78
t0 20 80
#N300
ll 80 110
lh 83 98
hh 90 108
ll 78 82
lh 84 95
hh 92 110
t0 20 105
#N200
ll 75 90
......@@ -45,11 +53,19 @@ t0 20 110
#N203
ll 77 95
ls 76 95
lh 78 95
lh 85 98
ss 76 95
sh 80 100
hh 89 102
t0 20 110
#D200
ll 76 92
ls 78 88
lh 82 95
ss 82 90
sh 85 95
hh 88 93
t0 20 110
#J303
ll 116 130
ls 118 135
......
src/plat_meff.txt → plat/plat_tm_mps.txt
View file @ 6c956856
......@@ -12,25 +12,33 @@ sh 60 85
hh 68 85
t0 20 80
#H102r002
ll 65 72
ls 67 74
lh 67 75
ll 65 76
ls 65 74
lh 68 73
ss 67 74
sh 68 81
hh 68 85
t0 20 80
#H105
ll 60 68
ls 65 74
lh 65 75
ss 67 74
sh 68 81
hh 68 85
t0 20 80
#H400
ll 65 88
ll 65 80
lh 67 80
hh 67 85
t0 20 80
#N300
ll 80 110
lh 92 115
ll 78 105
lh 92 108
hh 92 115
t0 20 105
#N200
ll 85 105
ll 80 92
ls 85 105
lh 88 99
ss 85 105
......@@ -43,15 +51,23 @@ lh 87 100
hh 90 115
t0 20 110
#N203
ll 85 105
ll 87 103
ls 85 105
lh 88 100
ss 85 105
sh 88 110
hh 88 110
t0 20 110
#D200
ll 79 86
ls 85 105
lh 84 87
ss 85 105
sh 85 87
hh 90 93
t0 25 100
#J303
ll 120 140
ll 120 148
ls 123 140
lh 125 145
ss 123 140
......
src/plat_meff_vec.txt → plat/plat_tm_mvec.txt
View file @ 6c956856
......@@ -7,7 +7,7 @@ t0 20 80
ll 60 85
ls 60 85
lh 60 85
ss 60 85
ss 60 85
sh 60 85
hh 60 85
t0 20 80
......@@ -16,7 +16,15 @@ ll 60 69
ls 60 68
lh 63 68
ss 60 67
sh 68 73
sh 62 68
hh 70 85
t0 20 80
#H105
ll 60 69
ls 60 68
lh 63 68
ss 60 67
sh 64 68
hh 70 85
t0 20 80
#H400
......@@ -25,8 +33,8 @@ lh 68 74
hh 68 80
t0 20 80
#N300
ll 80 110
lh 88 100
ll 80 90
lh 88 95
hh 95 115
t0 20 105
#N200
......@@ -39,23 +47,31 @@ hh 92 110
t0 20 105
#N202
ll 84 95
lh 87 97
lh 82 88
hh 92 110
t0 20 110
#N203
ll 85 95
ls 85 95
lh 85 95
lh 90 100
ss 85 95
sh 86 100
sh 90 100
hh 89 110
t0 20 110
#D200
ll 78 85
ls 80 90
lh 88 95
ss 85 95
sh 88 94
hh 88 98
t0 20 110
#J303
ll 120 140
ls 123 140
lh 123 140
lh 128 140
ss 123 140
sh 123 145
sh 128 140
hh 133 165
t0 25 170
#end
results_gevp/results.txt deleted 100644 → 0
View file @ 838f8293
Data analysis results
#================================#
Setup:
Ensembles analysed: ["J303", "N300", "H400"]
Data storage: /Users/ale/Desktop/data
Plateaux file: /Users/ale/automation/plat.txt
Results storage: /Users/ale/Desktop/results_gevp
Sector analysed: ["sh", "hh", "lh"]
Reweighting factor: false
Mass shift: false
t0 values: extracted from const.jl
Phi2 = 8t0m_pi^2: m_pi extracted from ens_db in const.jl
#================================#
Results at finite lattice spacings
Ensemble J303
$Z^{tm}_M \mu_c \sqrt{8t_0}$ = 3.089038311419621 +/- 0.06368579322874132
$m_{D_s} \sqrt{8t_0}$ = 4.0767724442326205 +/- 0.006996133249689365
$m_{D_s^*} \sqrt{8t_0}$ = 4.358985941909339 +/- 0.011323225630393369
$m_etac \sqrt{8t_0}$ = 6.198961563442169 +/- 0.010631637001414646
$m_{J/ψ} \sqrt{8t_0}$ = 6.438575482495065 +/- 0.011974541830682067
$m_D \sqrt{8t_0}$ = 3.9345824989703133 +/- 0.00788650891221943
$m_{D^*} \sqrt{8t_0}$ = 4.202091508641248 +/- 0.03026831950893543
Ensemble N300
$Z^{tm}_M \mu_c \sqrt{8t_0}$ = 3.0534116918858816 +/- 0.06423563962756408
$m_{D_s} \sqrt{8t_0}$ = 3.9819783981924606 +/- 0.008771802807377663
$m_{D_s^*} \sqrt{8t_0}$ = 4.307408596257999 +/- 0.023636314382326967
$m_etac \sqrt{8t_0}$ = 6.153439293650704 +/- 0.01079032654593637
$m_{J/ψ} \sqrt{8t_0}$ = 6.395244247561371 +/- 0.015177663706514
$m_D \sqrt{8t_0}$ = 3.9819783981924606 +/- 0.008771802807377663
$m_{D^*} \sqrt{8t_0}$ = 4.307408596257999 +/- 0.023636314382326967
Ensemble H400
$Z^{tm}_M \mu_c \sqrt{8t_0}$ = 3.082565030768681 +/- 0.07238415767395359
$m_{D_s} \sqrt{8t_0}$ = 3.981978398192199 +/- 0.020054878479507945
$m_{D_s^*} \sqrt{8t_0}$ = -1048.4099785502206 +/- 2.292257001461028
$m_etac \sqrt{8t_0}$ = 6.161662152741498 +/- 0.04555911269588078
$m_{J/ψ} \sqrt{8t_0}$ = 6.407190978248098 +/- 0.0779209000164276
$m_D \sqrt{8t_0}$ = 3.981978398192199 +/- 0.020054878479507945
$m_{D^*} \sqrt{8t_0}$ = -1048.4099785502206 +/- 2.292257001461028
#================================#
Results in the continuum limit
$Z^{tm}_M \mu_c \sqrt{8t_0}$ = 3.082309275644263 +/- 0.06438222973160301
mu_c(MeV) = 1472.6943874619703 +/- 15.576870185223473
$m_{D_s} \sqrt{8t_0}$ = 4.116369939134062 +/- 0.019235286187268108
m_Ds(MeV) = 1966.7575067764074 +/- 25.70202455882141
$m_{D_s^*} \sqrt{8t_0}$ = 784.7478160631761 +/- 1.699331215580599
m_Ds_star(MeV) = 374944.10876329575 +/- 4611.4192288057875
$m_etac \sqrt{8t_0}$ = 6.211881675661275 +/- 0.0384654831143396
m_etac(MeV) = 2967.97058997659 +/- 40.44155105716508
$m_{J/ψ} \sqrt{8t_0}$ = 6.4478198895589305 +/- 0.061288707527612424
m_jpsi(MeV) = 3080.6993437523615 +/- 47.48802002562802
$m_D \sqrt{8t_0}$ = 3.914784220954249 +/- 0.01987534965036722
m_D(MeV) = 1870.442007841328 +/- 24.46605658005924
$m_{D^*} \sqrt{8t_0}$ = 784.5253834787122 +/- 1.6997966664873079
m_D_star(MeV) = 374837.8328547369 +/- 4609.877645087271
src/aux obs computation/const.jl 0 → 100644
View file @ 6c956856
#========= ENSEMBLE DATABASE ========#
ens_db = Dict(
#"ens_id"=>[L, beta, is_deg?, m_pi, dtr]
"H102r002" => [32, 3.4, false, 0.15306, 2],
#"H102r001" => [32, 3.4, false, 0.15306, 2],
"H105" => [32, 3.4, false, 0.12151, 2],
"H101" => [32, 3.4, true, 0.17979, 2],
"H400" => [32, 3.46, true, 0.16345, 1],
"N200" => [48, 3.55, false, 0.09222, 1],
"N202" => [48, 3.55, true, 0.13407, 2],
"N203" => [48, 3.55, false, 0.11224, 1],
"D200" => [64, 3.55, false, 0.06611, 2],
"N300" => [48, 3.70, true, 0.10630, 1],
"J303" => [64, 3.70, false, 0.06514, 2]
)
trunc_db = Dict(
"H102r002" => nothing,
#"H102r001" => nothing,
"H105" => nothing,
"H101" => nothing,
"H400" => nothing,
"N200" => nothing,
"N202" => nothing,
"N203" => nothing,
"D200" => 1000,
"N300" => nothing,
"J303" => nothing
)
#PDG
const hc = 197.3269804 #MeV fm
const M_values = [1869.65, 2010.26, 1968.34, 2112.2, 2980.3, 3096.916] #MD, MD*, MDs, MDs*, \eta_c, J/\psi (MeV)
const M_error = [0.05, 0.05, 0.07, 0.4, 1.2, 0.011]
#1802.05243
const b_values = [3.40, 3.46, 3.55, 3.70, 3.85]
const b_values2 = [3.40, 3.46, 3.55, 3.70]
const ZM_data = [1.9684, 1.9935, 2.0253, 2.0630, 2.0814]
const ZM_error = [35, 27, 33, 38, 45] .* 1e-4
const ZM_tm_data = [2.6047, 2.6181, 2.6312, 2.6339, 2.6127]
const ZM_tm_error = [42, 33, 42, 48, 55] .* 1e-4
#1608.08900
const t0_data = [2.86, 3.659, 5.164, 8.595]
const t0_error = [11, 16, 18, 29] .* 1e-3
const t0_ph_value = #=[0.415] =# [0.4137]
const t0_ph_error = ones(1,1) .* #=4e-3 =# 3.6e-3
# 1808.09236
const ZA_data = [0.75642, 0.76169, 0.76979, 0.78378, 0.79667]
const ZA_err = [72, 93, 43, 47, 47] .*1e-5
#Covariance matrices (Uncorrelated)
const C1 = zeros(5, 5)
const C2 = zeros(5, 5)
const C3 = zeros(4, 4)
const C4 = zeros(6,6)
const C5 = zeros(5, 5)
for i = 1:6
C4[i,i] = M_error[i] ^ 2
if i<= 5
C1[i, i] = ZM_error[i] ^ 2
C2[i, i] = ZM_tm_error[i] ^ 2
C5[i, i] = ZA_err[i]^2
if i<=4
C3[i,i] = t0_error[i] ^ 2
end
end
end
#Definitions uwreal variables
const ZM = cobs(ZM_data, C1, "ZM values")
const ZM_tm = cobs(ZM_tm_data, C2, "ZM_tm values")
const M = cobs(M_values, C4, "charmed meson masses")
const mpi_ph = uwreal([134.8, 0.3], "mpi phys")
const t0_ph = uwreal([t0_ph_value[1], t0_ph_error[1]], "sqrt(8 t0) (fm)")
const t0_ = cobs(t0_data, C3, "t0")
const a_ = fill(t0_ph, length(t0_)) ./ sqrt.(8 .* t0_)
const Za = cobs(ZA_data, C5, "ZA")
#1802.05243
C = [[0.375369e-6, 0.429197e-6, -0.186896e-5] [0.429197e-6, 0.268393e-4, 0.686776e-4] [-0.186896e-5, 0.686776e-4, 0.212386e-3]]
z = cobs([0.348629, 0.020921, 0.070613], C, "z")
ZP(b::Float64) = z[1] + z[2] * (b - 3.79) + z[3] * (b - 3.79)^2
Mrat = uwreal([0.9148, 0.0088], "M/mhad")
#Mrat = 0.9148
# Renormalization constants, t0 and a as functions (f(beta))
zm(beta::Float64) = ZM[b_values .== beta][1]
zm_tm(beta::Float64) = Mrat / ZP(beta)
t0(beta::Float64) = t0_[b_values2 .== beta][1]
a(beta::Float64) = a_[b_values2 .== beta][1]
za(beta::Float64) = Za[b_values .== beta][1]
\ No newline at end of file
src/aux obs computation/md.jl 0 → 100644
View file @ 6c956856
#=
# The code produces BDIO files for md_info (see OrderedDict md_info) and Wilson observables (see OrderedDict obs_w)
# The BDIO files are saved in path_bdio_md and path_bdio_w respectivily. Each BDIO file contains information about a given ensemble.
# Plateaux for different ensembles are defined in txt files (path_plat_w)
# The dat files must be located in the following way
path_data_w/ensemble -> *mesons.dat (Wilson)
path_data_w/ensemble -> *ms.dat
path_rw/ensemble -> *.ms1.dat
path_md/ensemble -> *.pbp.dat
# Example: path_rw/H400 contains H400r001.ms1.dat and H400r002.ms1.dat
# The target value of phi4 can be modified
# Masses and decay constants are in terms of the reference scale t0
# See read_bdio.jl for more information
=#
using Base: String
using LaTeXStrings: length
using OrderedCollections
using juobs, BDIO, DelimitedFiles, ADerrors, LaTeXStrings
#================== PATHS + INFO ==================#
const path_data_w = "/Users/alessandroconigli/Lattice/data/aux_obs_data/wilson"
const path_plat_w = "/Users/alessandroconigli/Lattice/data/aux_obs_data/plat_wilson.txt"
const path_rw = "/Users/alessandroconigli/Lattice/data/aux_obs_data/rwf"
const path_md = "/Users/alessandroconigli/Lattice/data/aux_obs_data/md"
const ens_list = ["D200"] # ["H101", "H102r002", "H400", "N202", "N200", "N203", "N300", "J303"]
const path_bdio_md = "/Users/alessandroconigli/Lattice/data/bdio_charm/md"
const path_bdio_w = "/Users/alessandroconigli/Lattice/data/bdio_charm/wilson"
const phi4_target = uwreal([1.117, 0.012], "phi4_target")
const rwf = true
#========= INCLUDES ==========#
include("./types.jl")
include("./tools.jl")
include("./const.jl")
#======== GET ENSEMBLE INFORMATION FROM DATABASE ==========#
ensinfo = Vector{EnsInfo}(undef, length(ens_list))
for i in 1:length(ens_list)
ens = ens_list[i]
try
ensinfo[i]= EnsInfo(ens, ens_db[ens], trunc_db[ens])
catch
error("The ensemble id ", ens, " was not found in the const.jl ens_db database.
Please check the ensemble id or update the database")
end
end
#t0 [2.86, 3.659, 5.164, 8.595]
wpmm = Dict{String, Vector{Float64}}()
wpmm["H101"] = [-1.0, 4.0, -1.0, 14.0*2.86]
wpmm["H102r002"] = [-1.0, 4.0, -1.0, 14.0*2.86]
wpmm["H105"] = [-1.0, 4.0, -1.0, 14.0*2.86]
wpmm["H400"] = [-1.0, 4.0, -1.0, 14.0*3.659]
wpmm["N202"] = [-1.0, 4.0, -1.0, 14.0*5.164]
wpmm["N200"] = [-1.0, 4.0, -1.0, 14.0*5.164]
wpmm["N203"] = [-1.0, 4.0, -1.0, 14.0*5.164]
wpmm["N300"] = [-1.0, 4.0, -1.0, 14.0*8.595]
wpmm["J303"] = [-1.0, 4.0, -1.0, 14.0*8.595]
##
#======== ANALYSIS ==========#
gen_mulist = Vector{Vector{Vector{Float64}}}(undef, length(ens_list))
md_info = Vector(undef, length(ensinfo))
obs_tm_shifted = Vector(undef, length(ensinfo))
#obs
obs_w = Vector(undef, length(ensinfo))
for (k, ens) in enumerate(ensinfo)
# t0
t0_ens, YW, W = get_t0(path_data_w, ens, rw=rwf, path_rw=path_rw, pl=true)
# WILSON correlators
pp_w, pp_wW, _ = get_corr(path_data_w, ens, "G5", "G5", rw=rwf, path_rw=path_rw, legacy=true)
pp_w_d1, _, _ = get_corr(path_data_w, ens, "G5_d1", "G5_d1", rw=rwf, path_rw=path_rw, legacy=true)
pp_w_d2, _, _ = get_corr(path_data_w, ens, "G5_d2", "G5_d2", rw=rwf, path_rw=path_rw, legacy=true)
# WILSON masses
m_w = get_meff_w(pp_w, ens, pl=true)
kappa_list = getfield.(pp_w, :kappa)
mll_w = get_ll_w(kappa_list, m_w, ens.deg)
mls_w = ens.deg ? mll_w : get_ls_w(kappa_list, m_w, ens.deg)
phi2 = 8 * t0_ens * mll_w^2
phi4 = ens.deg ? 8 * t0_ens * 1.5 * mll_w^2 : 8 * t0_ens * (0.5 * mll_w^2 + mls_w^2)
# MASS SHIFTS
dSdm = read_pbp(path_md, ens.id)
ppll = get_ll_w(kappa_list, pp_w, ens.deg)
ppllW = get_ll_w(kappa_list, pp_wW, ens.deg)
ppll_d1 = get_ll_w(kappa_list, pp_w_d1, ens.deg)
ppll_d2 = get_ll_w(kappa_list, pp_w_d2, ens.deg)
dphi4_s1, dphi4_s2 = md_sea(phi4, dSdm, YW, W) .+ md_sea(phi4, dSdm, ppllW, W)
dphi2_s1, dphi2_s2 = md_sea(phi2, dSdm, YW, W) .+ md_sea(phi2, dSdm, ppllW, W)
dphi4_v1, dphi4_v2 = md_val(phi4, ppll, [ppll_d1, ppll_d2])
dphi2_v1, dphi2_v2 = md_val(phi2, ppll, [ppll_d1, ppll_d2])
if !ens.deg
ppls = get_ls_w(kappa_list, pp_w, ens.deg)
pplsW = get_ls_w(kappa_list, pp_wW, ens.deg)
ppls_d1 = get_ls_w(kappa_list, pp_w_d1, ens.deg)
ppls_d2 = get_ls_w(kappa_list, pp_w_d2, ens.deg)
aux1, aux2 = md_val(phi4, ppls, [ppls_d1, ppls_d2])
dphi4_v1 = dphi4_v1 + aux1
dphi4_v2 = dphi4_v2 + aux2
dphi4_s1, dphi4_s2 = (dphi4_s1, dphi4_s2) .+ md_sea(phi4, dSdm, pplsW, W)
end
# Dict: Deltam, phi2, phi4, t0 (shifted)
deltam = (phi4_target - phi4) / (2*dphi4_s1 + dphi4_s2 + dphi4_v1 + dphi4_v2)
dt0_s1, dt0_s2 = md_sea(t0_ens, dSdm, YW, W)
t0_ens_shifted = t0_ens + deltam * (2 * dt0_s1 + dt0_s2)
phi2_shifted = phi2 + deltam * (2*dphi2_s1 + dphi2_s2 + dphi2_v1 + dphi2_v2)
phi4_shifted = phi4 + deltam * (2*dphi4_s1 + dphi4_s2 + dphi4_v1 + dphi4_v2)
md_info[k] = OrderedDict(
"deltam" => deltam,
"t0_shifted" => t0_ens_shifted,
"phi2_shifted" => phi2_shifted,
"phi4_shifted" => phi4_shifted
)
obs_w[k] = OrderedDict(
"mll" => sqrt(8 * t0_ens) * mll_w,
"mls" => sqrt(8 * t0_ens) * mls_w,
"t0" => t0_ens,
"phi2" => phi2,
"phi4" => phi4
)
end
## WRITE BDIO
# MD INFO
for (k, obs) in enumerate(md_info) #loop over ensembles
ens = ensinfo[k].id
p = joinpath(path_bdio_md, string(ens, ".bdio"))
fb = BDIO_open(p, "w", ens)
uinfo = 0
for (key, value) in obs #loop over observables
write_uwreal(value, fb, uinfo)
uinfo += 1
end
BDIO_close!(fb)
end
# WILSON
for (k, obs) in enumerate(obs_w) #loop over ensembles
ens = ensinfo[k].id
p = joinpath(path_bdio_w, string(ens, ".bdio"))
fb = BDIO_open(p, "w", ens)
uinfo = 0
for (key, value) in obs #loop over observables
write_uwreal(value, fb, uinfo)
uinfo += 1
end
BDIO_close!(fb)
end
\ No newline at end of file
src/aux obs computation/tools.jl 0 → 100644
View file @ 6c956856
function read_data(path::String, ens::String, g1::String="G5", g2::String="G5"; legacy::Bool=false)
p = joinpath(path, ens)
rep = filter(x->occursin("mesons.dat", x), readdir(p, join=true))
if length(rep)!=0
length(rep)==1 ? (return read_mesons(rep[1], g1, g2, legacy=legacy, id=ens)) : (return read_mesons(rep, g1, g2, legacy=legacy, id=ens))
else
error("mesons.dat file not found for ensemble ", ens, " in path ", p)
end
end
function read_pbp(path::String, ens::String)
p = joinpath(path, ens)
rep = filter(x->occursin("pbp.dat", x), readdir(p, join=true))
if length(rep)!=0
return read_md.(rep)
else
error("pbp.dat file not found for ensemble ", ens, " in path ", p)
end
end
function read_rw(path::String, ens::String)
p = joinpath(path, ens)
rep = filter(x->occursin("ms1.dat", x), readdir(p, join=true))
if length(rep)!=0
try
length(rep) == 1 ? (return read_ms1(rep[1])) : (return read_ms1.(rep))
catch
length(rep) == 1 && length(rep)!=0 ? (return read_ms1(rep[1], v="1.4")) : (return read_ms1.(rep, v="1.4"))
end
else
error("ms1.dat file not found for ensemble ", ens, " in path ", p)
end
end
function read_t0(path::String, ens::String, dtr::Int64)
p = joinpath(path, ens)
rep = filter(x->occursin("ms.dat", x), readdir(p, join=true))
if length(rep)!=0
length(rep) == 1 ? (return read_ms(rep[1], dtr=dtr, id=ens)) : (return read_ms.(rep, dtr=dtr, id=ens))
else
error("ms.dat file not found for ensemble ", ens, " in path ", p)
end
end
function get_corr(path::String, ens::EnsInfo, g1::String="G5", g2::String="G5"; rw::Bool=false, legacy::Bool=false, path_rw::String="")
if path_rw == ""
p_rw = path
else
p_rw = path_rw
end
aux = read_data(path, ens.id, g1, g2, legacy=legacy)
if !isnothing(ens.trunc)
truncate_data!(aux, ens.trunc)
end
if !rw
obs = corr_obs.(aux, L=ens.L, info=true)
return (getindex.(obs, 1), getindex.(obs, 2))
else
obs = corr_obs.(aux, L=ens.L, rw=read_rw(p_rw, ens.id), info=true)
return (getindex.(obs, 1), getindex.(obs, 2), getindex.(obs, 3))
end
end
function get_t0(path::String, ens::EnsInfo; rw::Bool=false, pl::Bool=false, path_rw::String="")
if path_rw == ""
p_rw = path
else
p_rw = path_rw
end
data = read_t0(path, ens.id, ens.dtr)
if !isnothing(ens.trunc)
truncate_data!(data, ens.trunc)
end
plat = select_plateau_t0(ens)
if !rw
t0 = comp_t0(data, plat, L=ens.L, pl=pl, info=true)
return (getindex(t0, 1), getindex(t0, 2))
else
t0 = comp_t0(data, plat, L=ens.L, pl=pl, rw=read_rw(p_rw, ens.id), info=true)
return (getindex(t0, 1), getindex(t0, 2), getindex(t0, 3))
end
end
function get_mu(mu_list::Vector{Vector{Float64}}, deg::Bool)
mu_sorted = unique(sort(minimum.(mu_list)))
mul = mu_sorted[1]
deg ? mus = 0.0 : mus = mu_sorted[2]
muh = unique(maximum.(mu_list))
muh = filter(x-> x > mul && x > mus, muh)
return mul, mus, muh
end
function get_kappa(kappa_list::Vector{Vector{Float64}}, deg::Bool)
kappa_sorted = unique(sort(maximum.(kappa_list), rev=true))
kappal = kappa_sorted[1]
deg ? kappas = 0.0 : kappas = kappa_sorted[2]
return kappal, kappas
end
function get_ll(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 mul in mu && mu[1] == mu[2] #l-l
return obs[k]
end
end
end
function get_ll_w(kappa_list, obs::Vector, deg::Bool)
kappal, kappas = get_kappa(kappa_list, deg)
for k = 1:length(kappa_list)
kappa = kappa_list[k]
if kappal in kappa && kappa[1] == kappa[2] #l-l
return obs[k]
end
end
end
function get_ls(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 mul in mu && mus in mu #l-l
return obs[k]
end
end
end
function get_ls_w(kappa_list, obs::Vector, deg::Bool)
kappal, kappas = get_kappa(kappa_list, deg)
for k = 1:length(kappa_list)
kappa = kappa_list[k]
if kappal in kappa && kappas in kappa #l-l
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)
for k = 1:length(mu_list)
mu = mu_list[k]
if mul in mu && mu[1] != mu[2] && !(mus in mu) #l-h
push!(obs_lh, obs[k])
end
end
return obs_lh
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_sh(mu_list, obs::Vector{uwreal}, deg::Bool)
mul, mus, muh = get_mu(mu_list, deg)
obs_sh = Vector{uwreal}(undef, 0)
for k = 1:length(mu_list)
mu = mu_list[k]
if mus in mu && mu[1] != mu[2] && !(mul in mu)#s-h
push!(obs_sh, obs[k])
end
end
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]
for i =1:length(muh)
if muh[i] in mu && mu[1] == mu[2]#h-h
push!(obs_hh, obs[k])
end
end
end
return obs_hh
end
function get_meff_tm(obs::Vector{juobs.Corr}, ens::EnsInfo; pl::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
mu_list = getfield.(obs, :mu)
plat = select_plateau_tm(ens, mu_list)
println(plat)
return meff.(obs, plat, pl=pl, wpm=wpm)
end
function get_meff_vec(obs::Vector{juobs.Corr}, ens::EnsInfo; pl::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
mu_list = getfield.(obs, :mu)
plat = select_plateau_vec(ens, mu_list)
println(plat)
return meff.(obs, plat, pl=pl, wpm=wpm)
end
function get_meff_w(obs::Vector{juobs.Corr}, ens::EnsInfo; pl::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
kappa_list = getfield.(obs, :kappa)
plat = select_plateau_w(ens, kappa_list)
println(plat)
return meff.(obs, plat, pl=pl, wpm=wpm)
end
function get_f_tm(obs::Vector{juobs.Corr}, m::Vector{uwreal}, ens::EnsInfo; pl::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
mu_list = getfield.(obs, :mu)
plat = select_plateau_tm(ens, mu_list)
return dec_const_pcvc.(obs, plat, m, pl=pl, wpm=wpm)
end
function select_plateau_t0(ensinfo::EnsInfo)
ens = ensinfo.id
f = readdlm(path_plat_w)
head = String.(f[:,1])
delim = findall(x->occursin("#", x), head)
edelim = findfirst(x->occursin(ens, x), head)
cdelim = findfirst(x-> x.== edelim, delim)
del = delim[cdelim]+1 : delim[cdelim+1]-1
head_ = String.(f[del, 1])
plat_ = Int64.(f[del, 2:3])
plat = Vector{Int64}(undef, 2)
n = findfirst(x-> occursin("t0", x), head_)
plat[1] = plat_[n,1]
plat[2] = plat_[n,2]
return plat
end
function select_plateau_tm(ensinf::EnsInfo, mu_list)
ens =ensinf.id
deg = ensinf.deg
mul, mus, muh = get_mu(mu_list, deg)
f = readdlm(path_plat_tm)
head = String.(f[:, 1])
delim = findall(x-> occursin("#", x ), head)
edelim = findfirst(x-> occursin(ens, x), head)
cdelim = findfirst(x-> x.== edelim, delim)
del = delim[cdelim]+1 : delim[cdelim+1]-1
#del = delim
plat = Vector{Vector{Int64}}(undef, 0)
#plat = Dict()
plat_ = Int64.(f[del, 2:3])
head_ = String.(f[del, 1])
for k = 1:length(mu_list)
mu = mu_list[k]
if mul in mu && mu[1] != mu[2] && !(mus in mu) #heavy-light
n = findfirst(x-> occursin("lh", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif mul in mu && mu[1] == mu[2] #light-light
n = findfirst(x-> occursin("ll", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif mul in mu && mus in mu#strange-light
n = findfirst(x-> occursin("ls", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif mus in mu && mu[1] != mu[2] && !(mul in mu)#heavy-strange
n = findfirst(x-> occursin("sh", x ), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif mus in mu && mu[1] == mu[2] && !(mul in mu)#strange-strange
n = findfirst(x-> occursin("ss", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif !(mul in mu) && !(mus in mu) #heavy-heavy
n = findfirst(x-> occursin("hh", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
end
end
return plat
end
function select_plateau_vec(ensinf::EnsInfo, mu_list)
ens =ensinf.id
deg = ensinf.deg
mul, mus, muh = get_mu(mu_list, deg)
f = readdlm(path_plat_tm)
head = String.(f[:, 1])
delim = findall(x-> occursin("#", x ), head)
edelim = findfirst(x-> occursin(ens, x), head)
cdelim = findfirst(x-> x.== edelim, delim)
del = delim[cdelim]+1 : delim[cdelim+1]-1
#del = delim
plat = Vector{Vector{Int64}}(undef, 0)
#plat = Dict()
plat_ = Int64.(f[del, 2:3])
head_ = String.(f[del, 1])
for k = 1:length(mu_list)
mu = mu_list[k]
if mul in mu && mu[1] != mu[2] && !(mus in mu) #heavy-light
n = findfirst(x-> occursin("lh_vec", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif mul in mu && mu[1] == mu[2] #light-light
n = findfirst(x-> occursin("ll_vec", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif mul in mu && mus in mu#strange-light
n = findfirst(x-> occursin("ls_vec", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif mus in mu && mu[1] != mu[2] && !(mul in mu)#heavy-strange
n = findfirst(x-> occursin("sh_vec", x ), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif mus in mu && mu[1] == mu[2] && !(mul in mu)#strange-strange
n = findfirst(x-> occursin("ss_vec", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif !(mul in mu) && !(mus in mu) #heavy-heavy
n = findfirst(x-> occursin("hh_vec", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
end
end
return plat
end
function select_plateau_w(ensinf::EnsInfo, kappa_list)
ens = ensinf.id
deg = ensinf.deg
kappal, kappas = get_kappa(kappa_list, deg)
f = readdlm(path_plat_w)
head = String.(f[:, 1])
delim = findall(x-> occursin("#", x ), head)
edelim = findfirst(x-> occursin(ens, x), head)
cdelim = findfirst(x-> x.== edelim, delim)
del = delim[cdelim]+1 : delim[cdelim+1]-1
#del = delim
plat = Vector{Vector{Int64}}(undef, 0)
#plat = Dict()
plat_ = Int64.(f[del, 2:3])
head_ = String.(f[del, 1])
for k = 1:length(kappa_list)
kappa = kappa_list[k]
if kappal in kappa && kappa[1] == kappa[2] #light-light
n = findfirst(x-> occursin("ll", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif kappal in kappa && kappas in kappa#strange-light
n = findfirst(x-> occursin("ls", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif kappas in kappa && kappa[1] == kappa[2] && !(kappal in kappa)#strange-strange
n = findfirst(x-> occursin("ss", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
end
end
return plat
end
function match_muc(muh, m_lh, m_sh, target)
M = (2/3 .* m_lh .+ 1/3 .* m_sh)
par, chi2exp = lin_fit(muh, M)
muh_target = x_lin_fit(par, target)
return muh_target
end
function aic_weight(aic)
w_aux = exp.(-1/2 * aic /minimum(aic)) #/minimum(aic) -> +0.5*mininum(aic)
w_aux1=filter(x->x<10000, w_aux)
index = findall(x->x<10000, w_aux)
norm = sum(w_aux1)
w = w_aux1 ./ norm
return w, index
end
function aic_weight(c::Cat)
aic = c.aic
w_aux = exp.(-1/2 * aic /minimum(aic)) #/minimum(aic) -> +0.5*mininum(aic)
w_aux1=filter(x->x<10000, w_aux)
index = findall(x->x<10000, w_aux)
norm = sum(w_aux1)
w = w_aux1 ./ norm
return w, index
end
function aic_syst(aic, fit_res)
w,index = aic_weight(aic)
aux1 = sum( w .* fit_res[index].^2)
aux2 = sum( w .* fit_res[index])^2
return sqrt(abs(value(aux1 - aux2)))
end
function aic_syst(c::Cat)
w,index = aic_weight(c.aic)
aux1 = sum( w .* c.fit_res[index].^2)
aux2 = sum( w .* c.fit_res[index])^2
return sqrt(abs(value(aux1 - aux2)))
end
function aic_model_ave(aic, fit_res, wpm::Union{Nothing, Dict}=nothing)
w, index = aic_weight(aic)
res = sum( w.* fit_res[index])
isnothing(wpm) ? uwerr(res) : uwerr(res, wpm)
return res
end
function aic_model_ave(c::Cat, wpm::Union{Nothing, Dict}=nothing)
w, index = aic_weight(c)
res = sum( w.* getfield(c, :fit_res)[index])
isnothing(wpm) ? uwerr(res) : uwerr(res, wpm)
return res
end
function category_av(cat_arr::Vector{Cat})
best_res = aic_model_ave.(cat_arr)
syst = aic_syst.(cat_arr)
tot_err = err.(best_res).^2 .+ syst.^2
w_aux = 1 ./ (tot_err)
w = w_aux ./ sum(w_aux)
println(findmax(w)[2])
aux1 = sum( w .* best_res.^2)
aux2 = sum( w .* best_res)^2
sys1 = sqrt(abs(value(aux1 - aux2)))
av = sum( w .* best_res)
return av, sys1
end
function category_av_tot(cat_arr::Vector{Cat})
fit_res = vcat(getfield.(cat_arr, :fit_res)...)
aic = vcat(getfield.(cat_arr, :aic)...)
best = aic_model_ave(aic, fit_res)
syst = aic_syst(aic, fit_res)
return best, syst
end
function syst_av(syst_err::Vector{Float64})
aux = sum(syst_err.^2)
return sqrt(aux)
end
function cat_obs(obs::Vector{Vector{uwreal}})
N = length(obs[1])
if !all(length.(obs) .== N)
error("Dimension mismatch")
end
r = getindex.(obs, 1)
for i = 2:N
r = vcat(r, getindex.(obs, i))
end
return r
end
function Base.:*(x::uwreal, y::Array{Any})
N = size(y, 1)
return fill(x, N) .* y
end
function Base.:*(x::uwreal, y::Array{Float64})
N = size(y, 1)
return fill(x, N) .* y
end
function Base.:*(x::uwreal, y::Array{uwreal})
N = size(y, 1)
return fill(x, N) .* y
end
function Base.:+(x::uwreal, y::Vector{uwreal})
N = size(y, 1)
return fill(x, N) .+ y
end
function Base.:+(x::Float64, y::Vector{Float64})
N = size(y, 1)
return fill(x, N) .+ y
end
function Base.:+(x::Float64, y::Vector{uwreal})
N = size(y, 1)
return fill(x, N) .+ y
end
function Base.:-(x::Float64, y::Vector{Any})
N = size(y, 1)
return fill(x, N) .- y
end
function Base.:-(x::Float64, y::Vector{Float64})
N = size(y, 1)
return fill(x, N) .- y
end
function Base.:-(x::Float64, y::Vector{uwreal})
N = size(y, 1)
return fill(x, N) .- y
end
\ No newline at end of file
src/aux obs computation/types.jl 0 → 100644
View file @ 6c956856
mutable struct EnsInfo
id::String
L::Int64
beta::Float64
deg::Bool
mpi::Float64
dtr::Int64
trunc::Union{Int64, Vector{Int64}, Nothing}
function EnsInfo(ens_id::String, info::Vector{Float64})
id = ens_id
L = info[1]
beta = info[2]
deg = info[3]
mpi = info[4]
dtr = info[5]
return new(id, L, beta, deg, mpi, dtr, nothing)
end
function EnsInfo(ens_id::String, info::Vector{Float64}, trunc::Union{Nothing, Int64, Vector{Int64}})
id = ens_id
L = info[1]
beta = info[2]
deg = info[3]
mpi = info[4]
dtr = info[5]
return new(id, L, beta, deg, mpi, dtr, trunc)
end
end
mutable struct EnsObs
ensinfo::EnsInfo
mu_list:: Vector{Vector{Float64}}
is_pseudo::Bool
m_ll::Union{Nothing,uwreal}
m_ls::Union{Nothing,uwreal}
m_lh::Union{Nothing,Vector{uwreal}}
m_ss::Union{Nothing,uwreal}
m_sh::Union{Nothing,Vector{uwreal}}
m_hh::Union{Nothing,Vector{uwreal}}
f_ll::Union{Nothing,uwreal}
f_ls::Union{Nothing,uwreal}
f_lh::Union{Nothing,Vector{uwreal}}
f_ss::Union{Nothing,uwreal}
f_sh::Union{Nothing,Vector{uwreal}}
f_hh::Union{Nothing,Vector{uwreal}}
m_ll_vec::Union{Nothing,uwreal}
m_ls_vec::Union{Nothing,uwreal}
m_lh_vec::Union{Nothing,Vector{uwreal}}
m_ss_vec::Union{Nothing,uwreal}
m_sh_vec::Union{Nothing,Vector{uwreal}}
m_hh_vec::Union{Nothing,Vector{uwreal}}
f_ll_vec::Union{Nothing,uwreal}
f_ls_vec::Union{Nothing,uwreal}
f_lh_vec::Union{Nothing,Vector{uwreal}}
f_ss_vec::Union{Nothing,uwreal}
f_sh_vec::Union{Nothing,Vector{uwreal}}
f_hh_vec::Union{Nothing,Vector{uwreal}}
m_lh_match::Union{Nothing,uwreal}
m_sh_match::Union{Nothing,uwreal}
m_hh_match::Union{Nothing,uwreal}
f_lh_match::Union{Nothing,uwreal}
f_sh_match::Union{Nothing,uwreal}
f_hh_match::Union{Nothing,uwreal}
m_lh_vec_match::Union{Nothing,uwreal}
m_sh_vec_match::Union{Nothing,uwreal}
m_hh_vec_match::Union{Nothing,uwreal}
f_lh_vec_match::Union{Nothing,uwreal}
f_sh_vec_match::Union{Nothing,uwreal}
f_hh_vec_match::Union{Nothing,uwreal}
muh_target::Union{Nothing,uwreal}
a::Union{Nothing,uwreal}
t0::Union{Nothing,uwreal}
deltam::Union{Nothing, uwreal}
function EnsObs(ens::EnsInfo, mu::Vector{Vector{Float64}}, m_ps::Vector{uwreal}, f_ps::Vector{uwreal}, m_vec::Vector{uwreal}, f_vec::Vector{uwreal})
a = new()
a.ensinfo = ens
a.mu_list = mu
a.is_pseudo = true
a.m_ll = get_ll(a.mu_list, m_ps,a.ensinfo.deg)
a.ensinfo.deg ? a.m_ls = a.m_ll : a.m_ls = get_ls(a.mu_list, m_ps, a.ensinfo.deg)
a.ensinfo.deg ? a.m_ss = a.m_ll : a.m_ss = get_ss(a.mu_list, m_ps, a.ensinfo.deg)
a.m_lh = get_lh(a.mu_list, m_ps,a.ensinfo.deg)
a.ensinfo.deg ? a.m_sh = a.m_lh : a.m_sh = get_sh(a.mu_list, m_ps, a.ensinfo.deg)
a.m_hh = get_hh(a.mu_list, m_ps,a.ensinfo.deg)
a.f_ll = get_ll(a.mu_list, f_ps, a.ensinfo.deg)
a.ensinfo.deg ? a.f_ls = a.f_ll : a.f_ls = get_ls(a.mu_list, f_ps, a.ensinfo.deg)
a.ensinfo.deg ? a.f_ss = a.f_ll : a.f_ss = get_ss(a.mu_list, f_ps, a.ensinfo.deg)
a.f_lh = get_lh(a.mu_list, f_ps,a.ensinfo.deg)
a.ensinfo.deg ? a.f_sh = a.f_lh : a.f_sh = get_sh(a.mu_list, f_ps, a.ensinfo.deg)
a.f_hh = get_hh(a.mu_list, f_ps,a.ensinfo.deg)
a.m_ll_vec = get_ll(a.mu_list, m_vec,a.ensinfo.deg)
a.ensinfo.deg ? a.m_ls_vec = a.m_ll_vec : a.m_ls_vec = get_ls(a.mu_list, m_vec, a.ensinfo.deg)
a.ensinfo.deg ? a.m_ss_vec = a.m_ll_vec : a.m_ss_vec = get_ss(a.mu_list, m_vec, a.ensinfo.deg)
a.m_lh_vec = get_lh(a.mu_list, m_vec, a.ensinfo.deg)
a.ensinfo.deg ? a.m_sh_vec = a.m_lh_vec : a.m_sh_vec = get_sh(a.mu_list, m_vec, a.ensinfo.deg)
a.m_hh_vec = get_hh(a.mu_list, m_vec, a.ensinfo.deg)
a.f_ll_vec = get_ll(a.mu_list, f_vec, a.ensinfo.deg)
a.ensinfo.deg ? a.f_ls_vec = a.f_ll_vec : a.f_ls_vec = get_ls(a.mu_list, f_vec, a.ensinfo.deg)
a.ensinfo.deg ? a.f_ss_vec = a.f_ll_vec : a.f_ss_vec = get_ss(a.mu_list, f_vec, a.ensinfo.deg)
a.f_lh_vec = get_lh(a.mu_list, f_vec,a.ensinfo.deg)
a.ensinfo.deg ? a.f_sh_vec = a.f_lh_vec : a.f_sh_vec = get_sh(a.mu_list, f_vec, a.ensinfo.deg)
a.f_hh_vec = get_hh(a.mu_list, f_vec, a.ensinfo.deg)
a.m_lh_match = nothing
a.m_sh_match = nothing
a.m_hh_match = nothing
a.f_lh_match = nothing
a.f_sh_match = nothing
a.f_hh_match = nothing
a.m_lh_vec_match = nothing
a.m_sh_vec_match = nothing
a.m_hh_vec_match = nothing
a.f_lh_vec_match = nothing
a.f_sh_vec_match = nothing
a.f_hh_vec_match = nothing
a.muh_target = nothing
a.a = nothing
a.t0 = nothing
a.deltam = nothing
return a
end
function EnsObs(ens::EnsInfo, mu::Vector{Vector{Float64}}, m_ps::Vector{uwreal}, m_vec::Vector{uwreal})
a = new()
a.ensinfo = ens
a.mu_list = mu
a.is_pseudo = true
a.m_ll = get_ll(a.mu_list, m_ps,a.ensinfo.deg)
a.ensinfo.deg ? a.m_ls = a.m_ll : a.m_ls = get_ls(a.mu_list, m_ps, a.ensinfo.deg)
a.ensinfo.deg ? a.m_ss = a.m_ll : a.m_ss = get_ss(a.mu_list, m_ps, a.ensinfo.deg)
a.m_lh = get_lh(a.mu_list, m_ps,a.ensinfo.deg)
a.ensinfo.deg ? a.m_sh = a.m_lh : a.m_sh = get_sh(a.mu_list, m_ps, a.ensinfo.deg)
a.m_hh = get_hh(a.mu_list, m_ps,a.ensinfo.deg)
a.m_ll_vec = get_ll(a.mu_list, m_vec,a.ensinfo.deg)
a.ensinfo.deg ? a.m_ls_vec = a.m_ll_vec : a.m_ls_vec = get_ls(a.mu_list, m_vec, a.ensinfo.deg)
a.ensinfo.deg ? a.m_ss_vec = a.m_ll_vec : a.m_ss_vec = get_ss(a.mu_list, m_vec, a.ensinfo.deg)
a.m_lh_vec = get_lh(a.mu_list, m_vec, a.ensinfo.deg)
a.ensinfo.deg ? a.m_sh_vec = a.m_lh_vec : a.m_sh_vec = get_sh(a.mu_list, m_vec, a.ensinfo.deg)
a.m_hh_vec = get_hh(a.mu_list, m_vec, a.ensinfo.deg)
a.m_lh_match = nothing
a.m_sh_match = nothing
a.m_hh_match = nothing
a.m_lh_vec_match = nothing
a.m_sh_vec_match = nothing
a.m_hh_vec_match = nothing
a.muh_target = nothing
a.a = nothing
a.t0 = nothing
return a
end
end
mutable struct MatInfo
ensinfo::EnsInfo
mu::Vector{Float64}
mat_list::Vector{Matrix{uwreal}}
y0::Int64
function MatInfo(_ensinfo::EnsInfo, _mat_list::Array{Array{T,2} where T,1}, _mu::Vector{Float64}, _y0::Int64)
a = new()
a.ensinfo = _ensinfo
a.mu = _mu
a.mat_list = _mat_list
a.y0 = _y0
return a
end
end
mutable struct Cat
obs::Vector{uwreal}
phih::Vector{uwreal}
phih_ph::uwreal
#models::Vector{Vector{Function}}
x_tofit::Array{uwreal}
info::String
fit_param::Vector{Vector{uwreal}} #store all fit parameters
fit_res::Vector{uwreal} # store all fit results
chi2_corr::Vector{Float64} # store all chi2 corrected
n_param::Vector{Int64} # store the number of parameters. Maybe this is redundant
aic::Vector{Float64} # store the AIC value
function Cat(_obs::Vector{uwreal}, _xtofit::Array{uwreal}, _phi_ph::uwreal, _info::String)
a = new()
a.obs = _obs
#a.models = _models
a.info = _info
a.x_tofit = _xtofit
a.phih_ph = _phi_ph
a.fit_param = Vector{Vector{uwreal}}(undef, 0)
a.fit_res = Vector{uwreal}(undef, 0)
a.chi2_corr = Vector{Float64}(undef, 0)
a.n_param = Vector{Int64}(undef, 0)
a.aic = Vector{Float64}(undef, 0)
return a
end
end
\ No newline at end of file
src/gevp analysis/const.jl 0 → 100644
View file @ 6c956856
#========= ENSEMBLE DATABASE ========#
ens_db = Dict(
#"ens_id"=>[L, beta, is_deg?, m_pi, dtr]
"H102r002" => [32, 3.4, false, 0.15306, 2],
#"H102r001" => [32, 3.4, false, 0.15306, 2],
"H105" => [32, 3.4, false, 0.12151, 2],
"H101" => [32, 3.4, true, 0.17979, 2],
"H400" => [32, 3.46, true, 0.16345, 1],
"N200" => [48, 3.55, false, 0.09222, 1],
"N202" => [48, 3.55, true, 0.13407, 2],
"N203" => [48, 3.55, false, 0.11224, 1],
"D200" => [64, 3.55, false, 0.06611, 2],
"N300" => [48, 3.70, true, 0.10630, 1],
"J303" => [64, 3.70, false, 0.06514, 2]
)
trunc_db = Dict(
"H102r002" => nothing,
#"H102r001" => nothing,
"H105" => nothing,
"H101" => nothing,
"H400" => nothing,
"N200" => nothing,
"N202" => nothing,
"N203" => nothing,
"D200" => 1000,
"N300" => nothing,
"J303" => nothing
)
#PDG
const hc = 197.3269804 #MeV fm
const M_values = [1869.65, 2010.26, 1968.34, 2112.2, 2980.3, 3096.916] #MD, MD*, MDs, MDs*, \eta_c, J/\psi (MeV)
const M_error = [0.05, 0.05, 0.07, 0.4, 1.2, 0.011]
#1802.05243
const b_values = [3.40, 3.46, 3.55, 3.70, 3.85]
const b_values2 = [3.40, 3.46, 3.55, 3.70]
const ZM_data = [1.9684, 1.9935, 2.0253, 2.0630, 2.0814]
const ZM_error = [35, 27, 33, 38, 45] .* 1e-4
const ZM_tm_data = [2.6047, 2.6181, 2.6312, 2.6339, 2.6127]
const ZM_tm_error = [42, 33, 42, 48, 55] .* 1e-4
#1608.08900
const t0_data = [2.86, 3.659, 5.164, 8.595]
const t0_error = [11, 16, 18, 29] .* 1e-3
const t0_ph_value = #=[0.415] =# [0.4137]
const t0_ph_error = ones(1,1) .* #=4e-3 =# 3.6e-3
# 1808.09236
const ZA_data = [0.75642, 0.76169, 0.76979, 0.78378, 0.79667]
const ZA_err = [72, 93, 43, 47, 47] .*1e-5
#Covariance matrices (Uncorrelated)
const C1 = zeros(5, 5)
const C2 = zeros(5, 5)
const C3 = zeros(4, 4)
const C4 = zeros(6,6)
const C5 = zeros(5, 5)
for i = 1:6
C4[i,i] = M_error[i] ^ 2
if i<= 5
C1[i, i] = ZM_error[i] ^ 2
C2[i, i] = ZM_tm_error[i] ^ 2
C5[i, i] = ZA_err[i]^2
if i<=4
C3[i,i] = t0_error[i] ^ 2
end
end
end
#Definitions uwreal variables
const ZM = cobs(ZM_data, C1, "ZM values")
const ZM_tm = cobs(ZM_tm_data, C2, "ZM_tm values")
const M = cobs(M_values, C4, "charmed meson masses")
const mpi_ph = uwreal([134.8, 0.3], "mpi phys")
const t0_ph = uwreal([t0_ph_value[1], t0_ph_error[1]], "sqrt(8 t0) (fm)")
const t0_ = cobs(t0_data, C3, "t0")
const a_ = fill(t0_ph, length(t0_)) ./ sqrt.(8 .* t0_)
const Za = cobs(ZA_data, C5, "ZA")
#1802.05243
C = [[0.375369e-6, 0.429197e-6, -0.186896e-5] [0.429197e-6, 0.268393e-4, 0.686776e-4] [-0.186896e-5, 0.686776e-4, 0.212386e-3]]
z = cobs([0.348629, 0.020921, 0.070613], C, "z")
ZP(b::Float64) = z[1] + z[2] * (b - 3.79) + z[3] * (b - 3.79)^2
Mrat = uwreal([0.9148, 0.0088], "M/mhad")
#Mrat = 0.9148
# Renormalization constants, t0 and a as functions (f(beta))
zm(beta::Float64) = ZM[b_values .== beta][1]
zm_tm(beta::Float64) = Mrat / ZP(beta)
t0(beta::Float64) = t0_[b_values2 .== beta][1]
a(beta::Float64) = a_[b_values2 .== beta][1]
za(beta::Float64) = Za[b_values .== beta][1]
\ No newline at end of file
src/gevp analysis/func_comb.jl 0 → 100644
View file @ 6c956856
###################################
########## MASSS ##########
###################################
##FUNCTIONS
basemodel(x,p) = p[1] .+ p[2] .* x[:,2] .+ p[3] .* x[:,3] .+ p[4] .* x[:,1]
a2phi2(x) = x[:,1] .* x[:,2] # phi2*a^2/8t0
a2phih2(x) = x[:,1] .* x[:,3].^2 # phih^2*a^2/8t0
a4phih4(x) = x[:,1].^(2) .* x[:,3].^4 # phih^4*(a^2/8t0)^2
a4phih4_2(x) = x[:,1].^(2) .* x[:,3].^2 # (a^2/8t0)^2*phih^2
a4cutoff(x) = x[:,1].^(2) # (a^2/8t0)^2
a2loga(x) = x[:,1] .* log.(x[:,1]) # (a^2/8t0)^2 * log(a^2/8t0)
phih2loga(x) = x[:,3].^2 .* log.(x[:,1]) # phih^2*a^2/8t0
# full list
func_list = [a2phi2, a2phih2, a4phih4, a4phih4_2, a4cutoff]
label_list = ["a2l", "a2h2", "a4h4", "a4h2", "a4"]
func_map = [Bool.([i,j,k,m,n]) for i=0:1 for j=0:1 for k=0:1 for m=0:1 for n=0:1]
# a^2log(a) and phih^2log(a)
# func_list = [ a2phih2, a2loga, phih2loga]
# label_list = [ "a2h2", "a2loga", "h2loga"]
# func_map = [Bool.([i,j,k]) for i=0:1 for j=0:1 for k=0:1 ]
# no a2l
# func_list = [ a2phih2, a4phih4, a4phih4_2, a4cutoff]
# label_list = [ "a2h2", "a4h4", "a4h2", "a4"]
# func_map = [Bool.([i,j,k,m]) for i=0:1 for j=0:1 for k=0:1 for m=0:1]
##COMBINATIONS (LINEAR)
npar = length(func_list) # number of extra parameters
n_param = Vector(4:4)
label_cutoff = Vector{Vector{String}}(undef, 0)
push!(label_cutoff, ["a2"])
#f_lin definition
f_lin = Vector{Vector{Function}}(undef, npar+1) #f[i][j]-> i: number of extra parameters, j: combinations
f_lin[1] = Vector{Function}(undef, 1)
f_lin[1][1] = (x,p) -> basemodel(x,p)
#f_aux definition -> similar to f_lin but without vectorized opeartors -> useful for plots
f_aux = Vector{Vector{Function}}(undef, npar+1)
f_aux[1] = Vector{Function}(undef, 1)
f_aux[1][1] = (x,p) -> p[1] + p[2] * x[:,2] + p[3] * x[:,3] + p[4] * x[:,1]
#loop over combinations
for n = 2:npar+1
aux = filter(x->sum(x) == n-1, func_map)
f_lin[n] = Vector{Function}(undef, length(aux))
f_aux[n] = Vector{Function}(undef, length(aux))
k = 1
for a in aux
#vectorized function
push!(n_param, 4+n-1)
push!(label_cutoff, label_list[a])
f_lin[n][k] = (x,p) -> basemodel(x,p) .+ sum([p[i+4] for i=1:(n-1)] .* (fill(x,n-1) .|> func_list[a]))
f_aux[n][k] = (x,p) -> f_aux[1][1](x, p) + sum([p[i+4] for i=1:(n-1)] .* (fill(x,n-1) .|> func_list[a]))
k = k + 1
end
end
##COMBINATIONS (NONLINEAR)
f_non_lin = Vector{Vector{Function}}(undef, npar) #f[i][j]-> i: number of extra parameters, j: combinations
f_non_lin_aux = Vector{Vector{Function}}(undef, npar)
n_non_lin_param = Vector(5:npar+4) #param for each non-linear functions
for n = 1:npar
aux = filter(x->sum(x) == n, func_map)
f_non_lin[n] = Vector{Function}(undef, length(aux))
f_non_lin_aux[n] = Vector{Function}(undef, length(aux))
k = 1
for a in aux
#vectorized function
push!(n_param, 4+n)
f_non_lin[n][k] = (x,p) -> basemodel(x,p) .* (1 .+ sum([p[i+4] for i=1:(n)] .* (fill(x,n) .|> func_list[a])))
f_non_lin_aux[n][k] = (x,p) -> p[1] + p[2] * x[:,2] + p[3] * x[:,3] + p[4] * x[:,1] +
(p[1] + p[2] * x[:,2] + p[3] * x[:,3] + p[4] * x[:,1]) .* sum([p[i+4] for i=1:(n)] .* (fill(x,n) .|> func_list[a]))
k = k + 1
end
end
##APPEND LIN + NONLIN
f = f_lin
append!(f, f_non_lin)
append!(f_aux, f_non_lin_aux)
#n_param2 = vcat(n_param, n_non_lin_param)
########################################
########## DEC. CONST. ##########
########################################
phi2_sym = 0.744666
##FUNCTIONS
f_a2phi2(x) = x[:,1] .* x[:,2] # phi2*a^2/8t0
f_a2phih2(x) = x[:,1] .* (x[:,3]).^2 # phih^2*a^2/8t0
f_a4phih4(x) = x[:,1].^(2) .* x[:,3].^4 # phih^4*(a^2/8t0)^2
f_a4phih2(x) = x[:,1].^(2) .* x[:,3].^2 # phih^2*(a^2/8t0)^2
f_a4cutoff(x) = x[:,1].^(2) # (a^2/8t0)^2
f_a2loga(x) = x[:,1] .* log.(x[:,1]) # (a^2/8t0)^2 * log(a^2/8t0)
f_phih2loga(x) = x[:,3].^2 .* log.(x[:,1]) # phih^2*a^2/8t0
##COMBINATIONS (LINEAR)
func_map = [Bool.([i,j,k, m, n]) for i=0:1 for j=0:1 for k=0:1 for m=0:1 for n=0:1]
f_basemodel1(x,p) = p[1] .+ p[2] .* x[:,2] .+ p[3] .* (1.0 ./ sqrt.(x[:,3])) .+ p[4] .* x[:,1]
f_basemodel2(x,p) = p[1] .+ p[2] .* (3 .* phi2_sym - 2 .* x[:,2]) .+ p[3] .* (1.0 ./ sqrt.(x[:,3])) .+ p[5] .* x[:,1]
#CHIRAL LOGS
f_basemodel1_clog(x,p) = p[1] .+ p[2] .* x[:,2] .+ p[3] .* (1.0 ./ sqrt.(x[:,3])) .+ p[5] .* x[:,1] .+ p[4] .* (3 .* x[:,2] .*log.(x[:,2]) .+
(3 .*phi2_sym .- x[:,2]) .* log.(1.5 .*phi2_sym .- 0.5*x[:,2]) .+ (2 .* phi2_sym .- x[:,2]) ./ 3 .* log.(2 .*phi2_sym .- x[:,2]))
f_basemodel2_clog(x,p) = p[1] .+ p[2] .* (3 .* phi2_sym - 2 .* x[:,2]) .+ p[3] .* (1.0 ./ sqrt.(x[:,3])) .+ p[6] .* x[:,1] .+ 4 .* p[4] .*(
(1.5 .* phi2_sym .- 0.5.*x[:, 2]) .* log.(1.5 .* phi2_sym .- 0.5.*x[:, 2]) .+ (2 .* phi2_sym .- x[:,2]) ./ 3 .* log.((2 .* phi2_sym .- x[:,2])))
f_func_list = [f_a2phi2, f_a2phih2, f_a4phih4, f_a4phih2, f_a4cutoff]
#f_func_list = [f_a2loga, f_a2phih2, f_a4phih4, f_a4phih2, f_a4cutoff]
f_npar = length(f_func_list) # number of extra parameters
## FLIN + FAUX
f_f_lin1 = Vector{Vector{Function}}(undef, f_npar+1) #f[i][j]-> i: number of extra parameters, j: combinations
f_f_lin2 = Vector{Vector{Function}}(undef, f_npar+1) #f[i][j]-> i: number of extra parameters, j: combinations
f_f_aux1 = Vector{Vector{Function}}(undef, f_npar+1)
f_f_aux2 = Vector{Vector{Function}}(undef, f_npar+1)
f_f_lin1[1] = Vector{Function}(undef, 1)
f_f_lin2[1] = Vector{Function}(undef, 1)
f_f_aux1[1] = Vector{Function}(undef, 1)
f_f_aux2[1] = Vector{Function}(undef, 1)
f_f_lin1[1][1] = (x,p) -> f_basemodel1(x, p)
f_f_lin2[1][1] = (x,p) -> f_basemodel2(x, p)
f_f_aux1[1][1] = (x,p) -> p[1] + p[2] * x[:,2] + p[3] * (1.0 ./ sqrt.(x[:,3])) + p[4] * x[:,1]
f_f_aux2[1][1] = (x,p) -> p[1] + p[2] * (3*phi2_sym - 2*x[:,2]) + p[3] * (1.0 ./ sqrt.(x[:,3])) + p[4] * x[:,1]
# FLIN + FAUX (CLOGS) [hep-ph/9206230]
f_f_lin1_clog = Vector{Vector{Function}}(undef, f_npar+1) #f[i][j]-> i: number of extra parameters, j: combinations
f_f_lin2_clog = Vector{Vector{Function}}(undef, f_npar+1) #f[i][j]-> i: number of extra parameters, j: combinations
f_f_aux1_clog = Vector{Vector{Function}}(undef, f_npar+1)
f_f_aux2_clog = Vector{Vector{Function}}(undef, f_npar+1)
f_f_lin1_clog[1] = Vector{Function}(undef, 1)
f_f_lin2_clog[1] = Vector{Function}(undef, 1)
f_f_aux1_clog[1] = Vector{Function}(undef, 1)
f_f_aux2_clog[1] = Vector{Function}(undef, 1)
f_f_lin1_clog[1][1] = (x,p) -> f_basemodel1_clog(x, p)
f_f_lin2_clog[1][1] = (x,p) -> f_basemodel2_clog(x, p)
f_f_aux1_clog[1][1] = (x,p) -> p[1] + p[2] * x[:,2] + p[3] * (1.0 ./ sqrt.(x[:,3])) + p[5] * x[:,1] + p[4] * (3 * x[:,2] .*log.(x[:,2]) +
(3 *phi2_sym - x[:,2]) .* log.(1.5 *phi2_sym - 0.5*x[:,2]) + (2 * phi2_sym - x[:,2]) / 3 .* log.(2 *phi2_sym - x[:,2]))
f_f_aux2_clog[1][1] = (x,p) -> p[1] + p[2] * (3 * phi2_sym - 2 * x[:,2]) + p[3] * (1.0 ./ sqrt.(x[:,3])) + p[6] * x[:,1] + 4 * p[4] *(
(1.5 * phi2_sym - 0.5*x[:, 2]) .* log.(1.5 * phi2_sym - 0.5*x[:, 2]) + (2 * phi2_sym - x[:,2]) / 3 .* log.((2 * phi2_sym - x[:,2])))
## loop over combinations
for n = 2:f_npar+1
aux = filter(x->sum(x) == n-1, func_map)
f_f_lin1[n] = Vector{Function}(undef, length(aux))
f_f_lin2[n] = Vector{Function}(undef, length(aux))
f_f_aux1[n] = Vector{Function}(undef, length(aux))
f_f_aux2[n] = Vector{Function}(undef, length(aux))
f_f_lin1_clog[n] = Vector{Function}(undef, length(aux))
f_f_lin2_clog[n] = Vector{Function}(undef, length(aux))
f_f_aux1_clog[n] = Vector{Function}(undef, length(aux))
f_f_aux2_clog[n] = Vector{Function}(undef, length(aux))
k=1
for a in aux
#vectorized function
f_f_lin1[n][k] = (x,p) -> f_basemodel1(x,p) .+ sum([p[i+5] for i=1:2:2*(n-1)] .* (fill(x,n-1) .|> f_func_list[a]))
f_f_aux1[n][k] = (x,p) -> f_f_aux1[1][1](x, p) +
sum([p[i+5] for i=1:2:2*(n-1)] .* (fill(x,n-1) .|> f_func_list[a]))
f_f_lin2[n][k] = (x,p) -> f_basemodel2(x,p) .+ sum([p[i+6] for i=1:2:2*(n-1)] .* (fill(x,n-1) .|> f_func_list[a]))
f_f_aux2[n][k] = (x,p) -> f_f_aux2[1][1](x, p) +
sum([p[i+6] for i=1:2:2*(n-1)] .* (fill(x,n-1) .|> f_func_list[a]))
f_f_lin1_clog[n][k] = (x,p) -> f_basemodel1_clog(x,p) .+ sum([p[i+6] for i=1:2:2*(n-1)] .* (fill(x,n-1) .|> f_func_list[a]))
f_f_aux1_clog[n][k] = (x,p) -> f_f_aux1_clog[1][1](x,p) +
sum([p[i+6] for i=1:2:2*(n-1)] .* (fill(x,n-1) .|> f_func_list[a]))
f_f_lin2_clog[n][k] = (x,p) -> f_basemodel2_clog(x,p) .+ sum([p[i+7] for i=1:2:2*(n-1)] .* (fill(x,n-1) .|> f_func_list[a]))
f_f_aux2_clog[n][k] = (x,p) -> f_f_aux2_clog[1][1](x,p) +
sum([p[i+7] for i=1:2:2*(n-1)] .* (fill(x,n-1) .|> f_func_list[a]))
k = k+1
end
end
f_npar_comb = 2 * f_npar
f_n_param_comb = collect(5:2:f_npar_comb+5)
f_n_param_comb_clog = collect(6:2:f_npar_comb+6)
\ No newline at end of file
src/gevp analysis/gevp_fit.jl 0 → 100644
View file @ 6c956856
#################################
# This file was created by Alessandro Conigli - 2022
# The pourpose of this code is to read BDIO files
# and perform chiral+continuum extrapolation.
#
##################################
using Base: String
using LaTeXStrings: length
using OrderedCollections
using juobs, BDIO, DelimitedFiles, ADerrors, LaTeXStrings, PyPlot, Statistics
#============= SET UP VARIABLES ===========#
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["text.usetex"] = false
rcParams["mathtext.fontset"] = "cm"
rcParams["font.size"] =12
rcParams["axes.labelsize"] =20
rcParams["axes.titlesize"] = 18
plt.rc("text", usetex=false) # set to true once you install latex
#rcParams.update(PyPlot.matplotlib.mpl.rcParamsDefault)
#================== PATHS ==================#
# data by ale (carfull with the read_BDIO that you include)
const path_bdio = "/Users/alessandroconigli/Lattice/data/bdio_charm"
#const path_bdio = "/Users/alessandroconigli/Lattice/data/bdio_charm/from_gevp/tm_shifted/ale"
const path_bdio_tm = joinpath(path_bdio, "tm_shifted")
const path_bdio_md = joinpath(path_bdio, "md")
# data by javier (carful with the read_BDIO that you include. do not use the one in this folder)
# const path_bdio = "/Users/alessandroconigli/MyDrive/phd/analysis/charm_gevp_2022/analysis_with_javier_data/bdio" #PATH to BDIO folder
# const path_bdio_tm = joinpath(path_bdio, "tm_shifted")
# const path_bdio_md = joinpath(path_bdio, "md")
# by reading javier data you have deactivated the reading of vector decay constants
const path_plot = "/Users/alessandroconigli/MyDrive/phd/analysis/charm_gevp_2022/plots"
#========= INCLUDES ==========#
include("./const.jl")
include("./types.jl")
include("./tools.jl")
include("./func_comb.jl")
# include("/Users/alessandroconigli/MyDrive/phd/analysis/charm_gevp_2022/analysis_with_javier_data/read_bdio.jl")
include("./read_bdio.jl") # this is the one to include with data by ale
include("./gevp_plotter.jl")
#======== GET ENSEMBLE INFORMATION FROM DATABASE ==========#
enslist = filter(x -> occursin(".bdio", x), readdir(path_bdio_tm))
Nens = length(enslist)
ensname = [enslist[k][1:end-5] for k = 1:Nens]
ensinfo = Vector{EnsInfo}(undef, length(ensname))
for (i, ens) in enumerate(ensname)
try
ensinfo[i]= EnsInfo(ens, ens_db[ens], trunc_db[ens])
catch
error("The ensemble id ", ens, " was not found in the const.jl ens_db database.
Please check the ensemble id or update the database")
end
end
wpmm = Dict{String, Vector{Float64}}()
wpmm["H101"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["H102r002"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["H400"] = [-1.0, 1.5, -1.0, -1.0]
wpmm["N202"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["N200"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["N203"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["N300"] = [-1.0, 1.5, -1.0, -1.0]
wpmm["J303"] = [-1.0, 2.0, -1.0, -1.0]
# READ MASSES + DECAY CONSTANTS
muh = Vector{Vector{uwreal}}(undef, Nens)
mlh = Vector{Vector{uwreal}}(undef, Nens)
msh = Vector{Vector{uwreal}}(undef, Nens)
mhh = Vector{Vector{uwreal}}(undef, Nens)
mlh_v = Vector{Vector{uwreal}}(undef, Nens)
msh_v = Vector{Vector{uwreal}}(undef, Nens)
mhh_v = Vector{Vector{uwreal}}(undef, Nens)
flh = Vector{Vector{uwreal}}(undef, Nens)
fsh = Vector{Vector{uwreal}}(undef, Nens)
fhh = Vector{Vector{uwreal}}(undef, Nens)
flh_v = Vector{Vector{uwreal}}(undef, Nens)
fsh_v = Vector{Vector{uwreal}}(undef, Nens)
fhh_v = Vector{Vector{uwreal}}(undef, Nens)
phih = Vector{Dict{String, Vector{uwreal}}}(undef, Nens)
a28t0 = Vector{uwreal}(undef, Nens)
phi2 = Vector{uwreal}(undef, Nens)
phi4 = Vector{uwreal}(undef, Nens)
for (k, ens) in enumerate(enslist)
p = joinpath(path_bdio_tm, ens)
mlh[k] = read_BDIO(p, "tm", "mlh")
msh[k] = read_BDIO(p, "tm", "msh")
mhh[k] = read_BDIO(p, "tm", "mhh")
mll = read_BDIO(p, "tm", "mll")[1]
mls = read_BDIO(p, "tm", "mls")[1]
mlh_v[k] = read_BDIO(p, "tm", "mlh_star")
msh_v[k] = read_BDIO(p, "tm", "msh_star")
mhh_v[k] = read_BDIO(p, "tm", "mhh_star")
flh[k] = read_BDIO(p, "tm", "flh")
fsh[k] = read_BDIO(p, "tm", "fsh")
fhh[k] = read_BDIO(p, "tm", "fhh")
flh_v[k] = read_BDIO(p, "tm", "flh_star") .* za.(ensinfo[k].beta)
fsh_v[k] = read_BDIO(p, "tm", "fsh_star") .* za.(ensinfo[k].beta)
fhh_v[k] = read_BDIO(p, "tm", "fhh_star") .* za.(ensinfo[k].beta)
phih_fa = (2 * mlh[k] + msh[k]) / 3
phih_sfa = (3 * phih_fa + 6 * mlh_v[k] + 3 * msh_v[k]) / 12
phih_eta = mhh[k]
phih[k] = Dict(
"phih_fa" => phih_fa,
"phih_sfa" => phih_sfa,
"phih_eta" => phih_eta
)
muh[k] = read_BDIO(p, "tm", "muh")
zp_aux = fill(ZP(ensinfo[k].beta), length(muh[k]))
muh[k] ./= zp_aux
p = joinpath(path_bdio_md, ens)
a28t0[k] = 1 / ( 8 * read_BDIO(p, "md", "t0_shifted")[1])
phi2[k] = mll^2
phi4[k] = mls^2 + 0.5 * mll^2
end
## X VARIABLE [a^2/(8t0) phi_2 phi_h]
x_fa = [a28t0 phi2 getindex.(getindex.(phih, "phih_fa"), 1)]
x_sfa = [a28t0 phi2 getindex.(getindex.(phih, "phih_sfa"), 1)]
x_eta = [a28t0 phi2 getindex.(getindex.(phih, "phih_eta"), 1)]
for i = 2:length(muh[1])
x_fa = vcat(x_fa, [a28t0 phi2 getindex.(getindex.(phih, "phih_fa"), i)])
x_sfa = vcat(x_sfa, [a28t0 phi2 getindex.(getindex.(phih, "phih_sfa"), i)])
x_eta = vcat(x_eta, [a28t0 phi2 getindex.(getindex.(phih, "phih_eta"), i)])
end
# PHYSICAL VALUES
phih_eta_ph = M[5] * t0_ph / hc
phih_fa_ph = (2 * M[1] + M[3]) * t0_ph / (3 * hc)
phih_sfa_ph = (2 * M[1] + M[3] + 6 * M[2] + 3 * M[4]) * t0_ph / (12 * hc)
phi2_ph = (t0_ph * mpi_ph / hc)^2 #isospin symmetric mpi = 134.8 MeV
#============= CREATE MASS_CAT =================#
# these are categories for runs with 3 value of the charm mass - hard coded
# s_obs_mass = ["muc", "m_etac", "m_jpsi", "mD", "mDs", "mD_star", "mDs_star"]
# obs = cat_obs.([muh, mhh, mhh_v, mlh, msh, mlh_v, msh_v])
# s_match = ["fl_ave", "etac", "sp_fl_ave"]
# match = [x_fa, x_eta, x_sfa]
# phih_ph = [phih_fa_ph, phih_eta_ph, phih_sfa_ph]
# this is for general number of charm masses
# keep in mind that the order of the data is differente: above it was for each mh all the ensembles,
# here it is for each ensembles all the mh values. No differences in fits but xdata and ydata need to be correctly ordered
s_obs_mass = ["muc", "m_etac", "m_jpsi", "mD", "mDs", "mD_star", "mDs_star"]
obs = cat_obs_vcat.([muh, mhh, mhh_v, mlh, msh, mlh_v, msh_v])
s_match = ["fl_ave", "etac", "sp_fl_ave"]
[uwerr.(obs[i]) for i in 1:length(obs)]
ens_id = getindex.(ensembles.(obs[1]),2)
a28t0_tot = Vector{Vector{uwreal}}(undef, length(ensinfo))
phi2_tot = Vector{Vector{uwreal}}(undef, length(ensinfo))
for (k, ens) in enumerate(ensinfo)
idx = findall(x->x==ens.id, ens_id)
a28t0_tot[k] = fill(a28t0[k], length(idx))
phi2_tot[k] = fill(phi2[k], length(idx))
end
x_fa = [vcat(a28t0_tot...) vcat(phi2_tot...) vcat(getindex.(phih, "phih_fa")...) ]
x_eta = [vcat(a28t0_tot...) vcat(phi2_tot...) vcat(getindex.(phih, "phih_eta")...) ]
x_sfa = [vcat(a28t0_tot...) vcat(phi2_tot...) vcat(getindex.(phih, "phih_sfa")...) ]
match = [x_fa, x_eta, x_sfa]
phih_ph = [phih_fa_ph, phih_eta_ph, phih_sfa_ph]
# creating mass categories
mass_cat = Vector{Vector{Cat}}(undef, 0)
for (i, o) in enumerate(obs)
aux = Vector{Cat}(undef, 0)
for (j, m) in enumerate(match)
push!(aux, Cat(o, m, phih_ph[j], string(s_obs_mass[i], "_", s_match[j])))
end
push!(mass_cat, aux)
end
##
#============= FIT MESON MASSES =============#
@info("Fitting meson masses for charm quark matching")
# mD
fit_over_cat!(mass_cat[4], f, phi2_ph, n_param, weight_dof=2, correlated_fit=false, Mrat=uwreal(1.0), small_sample=true)
# mDs
fit_over_cat!(mass_cat[5], f, phi2_ph, n_param, weight_dof=2, correlated_fit=false, Mrat=uwreal(1.0), small_sample=true)
# mD_star
fit_over_cat!(mass_cat[6], f, phi2_ph, n_param, weight_dof=2, correlated_fit=false, Mrat=uwreal(1.0), small_sample=true)
# mDs_star
fit_over_cat!(mass_cat[7], f, phi2_ph, n_param, weight_dof=2, correlated_fit=false, Mrat=uwreal(1.0), small_sample=true)
##
s = [L"$\sqrt{8t_0} Z^{tm}_M \mu_c$", L"$\sqrt{8t_0}M_{\eta^{(c)}_c}$", L"$\sqrt{8t_0}M_{J/\psi}$",
L"$\sqrt{8t_0}M_D$", L"$\sqrt{8t_0}M_{D_s}$", L"$\sqrt{8t_0}M_{D^{*}}$", L"$\sqrt{8t_0}M_{D^{*}_s}$"]
idx_cat = 4 # 4,5
for i in 2:2
#plot_cat(mass_cat[idx_cat][i], s[idx_cat] , save=false)
#plot_chi2_vs_exp(mass_cat[idx_cat][i], label_cutoff ) #, p=path_plot)
plot_chi_charm(mass_cat[idx_cat][i], vcat(f_aux...), s[idx_cat], phi2_ph, ensinfo, mrat=uwreal(1.0), p=path_plot)
#plot_cl_charm(mass_cat[idx_cat][i], vcat(f_aux...), s[idx_cat], phi2_ph, ensinfo, mrat=uwreal(1.0))#, p=path_plot)
end
##
#============= FIT CHARM MASS ==============#
@info("Fitting charm quark mass")
fit_over_cat!(mass_cat[1], f, phi2_ph, n_param, weight_dof=2, correlated_fit=false, Mrat=Mrat, small_sample=true )
## plot charm
for i in 1:3
plot_cat(mass_cat[1][i], L"$M_c^{\mathrm{RGI}}(N_f=3)\mathrm{MeV}$" , save=false)
plot_chi2_vs_exp(mass_cat[1][i], label_cutoff)#, p=path_plot)
plot_chi_charm(mass_cat[1][i], vcat(f_aux...), L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, mrat=Mrat, p=path_plot)
plot_cl_charm(mass_cat[1][i], vcat(f_aux...), L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, mrat=Mrat)#, p=path_plot)
end
close("all")
## plot for multiple categories
plot_hist_tot([mass_cat[1][1], mass_cat[1][2]], L"$M^\mathrm{RGI}_c\,(N_f=3)\,(\mathrm{MeV})$" , p=path_plot)
plot_all_cl_charm([mass_cat[1][2], mass_cat[1][1], mass_cat[1][3]], vcat(f_aux...),L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, mrat=Mrat, n=[0,0,1] )#, p=path_plot)
plot_cat_tot(mass_cat[1][1:2], L"$M^\mathrm{RGI}_c\,(N_f=3)\,(\mathrm{MeV})$" , save=true)
##
#============== CREATE DECAY CAT =============#
s_obs_f = ["fD", "fDs", "f_etac"]
s_match = ["fl_ave", "etac", "sp_fl_ave"]
phih_ph = [phih_fa_ph, phih_eta_ph, phih_sfa_ph]
## for uncorrelated fit the orginal order ens1...ensN,ens1...ensN,ens1...ensN is fine
obs = cat_obs.([flh, fsh, fhh])
match = [x_fa, x_eta, x_sfa]
## for correlated fits the order must be different: ens11..ens13, ens21...ens23,...,ensN1..ensN3
obs = [vcat(flh...), vcat(fsh...), vcat(fhh...)]
[uwerr.(obs[i]) for i in 1:length(obs)]
# preparing xdata for fitting in correct order
ens_id = getindex.(ensembles.(obs[1]),2)
a28t0_tot = Vector{Vector{uwreal}}(undef, length(ensinfo))
phi2_tot = Vector{Vector{uwreal}}(undef, length(ensinfo))
for (k, ens) in enumerate(ensinfo)
idx = findall(x->x==ens.id, ens_id)
a28t0_tot[k] = fill(a28t0[k], length(idx))
phi2_tot[k] = fill(phi2[k], length(idx))
end
x_fa = [vcat(a28t0_tot...) vcat(phi2_tot...) vcat(getindex.(phih, "phih_fa")...) ]
x_eta = [vcat(a28t0_tot...) vcat(phi2_tot...) vcat(getindex.(phih, "phih_eta")...) ]
x_sfa = [vcat(a28t0_tot...) vcat(phi2_tot...) vcat(getindex.(phih, "phih_sfa")...) ]
match = [x_fa, x_eta, x_sfa]
## create categories
f_cat = Vector{Vector{Cat}}(undef, 0)
for (i, o) in enumerate(obs)
aux = Vector{Cat}(undef, 0)
for (j, m) in enumerate(match)
push!(aux, Cat(o, m, phih_ph[j], string(s_obs_f[i], "_", s_match[j])))
end
push!(f_cat, aux)
end
##########################
## TEST ON fD/fDs
#f_ratio_cat = Vector{Cat}(undef, 0)
#for (j,m) in enumerate(match)
# push!(f_ratio_cat, Cat(obs[2]./obs[1], m, phih_ph[j], string("fds_over_fd", "_", s_match[j]) ))
#end
obs_ratio = obs[2]./obs[1]
idx_no_sym = findall(x->x!=1.0, value.(obs_ratio))
obs_ratio = obs_ratio[idx_no_sym]
#obs_ratio[1:3] .+= 0.015
xratio = x_fa[idx_no_sym, 1:3]
f_ratio_cat = Cat(obs_ratio, xratio, phih_ph[1], string("fds_over_fd") )
ratio_model(x,p) = 1 .+ p[1].*(3. .*phi2_sym .- 3. .* x[:,2]) .+ p[2].*x[:,1]
ratio_model_aux = (x,p) -> 1. + p[1]*(3. *phi2_sym - 3. * x[:,2]) + p[2]*x[:,1]
par, chi_vs_chi2exp = fit_routine(ratio_model, value.(f_ratio_cat.x_tofit), f_ratio_cat.obs, 2)
xx = [0.0 phi2_ph]
ress = ratio_model_aux(xx, par)[1]
push!(f_ratio_cat.fit_res, ress )
push!(f_ratio_cat.fit_param, par)
## plot ratio
idx_ens_nosim = findall( x-> x==false, getfield.(ensinfo, :deg))
obs_ratio .+= 0.01
plot_chi_fds_ratio(f_ratio_cat, ratio_model_aux, L"$f_{D_s}/f_D$", phi2_ph, ensinfo[idx_ens_nosim], p=path_plot)
##
function plot_chi_fds_ratio(c::Cat, ff::Function, ylab::LaTeXString, phi2_ph::uwreal, ensinfo::Vector{EnsInfo}; p::Union{String, Nothing}=nothing )
color = ["orange", "darkgreen", "magenta", "navy"]
fmt = ["^", "v", "<", ">"]
a_sp = [0.087, 0.077, 0.065, 0.05]
phi2_max = phi2_sym#maximum(value.(phi2))
upar = getfield(c, :fit_param)[1]
phih_2_use = value(c.phih_ph)
res = getfield(c, :fit_res)[1]
xx = [fill(1e-8, 100) Float64.(range(0.01, phi2_max, length=100)) ]
yy = ff(xx, upar)
#Project one value of muc_i (i=1,2,3)
n_unique = unique(i-> c.x_tofit[i,2], 1:length(c.x_tofit[:,2]))
p2 = c.x_tofit[n_unique, 2]
a2 = c.x_tofit[n_unique, 1]
pc = c.x_tofit[n_unique, 3]
x = [a2 p2 pc]
x2 = [a2 p2 fill(phih_2_use, length(p2))] #physical value phih
y = ( c.obs[n_unique] - ff(x, upar) + ff(x2, upar)) #Projection
uwerr.(yy)
uwerr(res)
uwerr.(y)
#Plotting
figure(figsize=(10,6))
plot(phi2_sym, 1., marker="X", ms=10, color="black", label= L"$M_\pi=M_K$")
fill_between(xx[:,2], value.(yy)+ err.(yy), value.(yy)-err.(yy), alpha=0.3, color="royalblue") #CL extrapolation (band)
errorbar(value(phi2_ph), value(res), err(res), fmt="s", color="red", label=L"$\phi_2 = \phi^\mathrm{phys}_2$", capsize=5, ms=12)#CL extrapolation
for (k,b) in enumerate(sort(unique(getfield.(ensinfo, :beta))))
n_ = findall(x->x.beta == b, ensinfo)
a2_aux = mean(value.(a2[n_]))
errorbar(value.(x[n_,2]), value.(y[n_]), err.(y[n_]), fmt=fmt[k], label=string(L"$a \approx$", a_sp[k], " fm"), color=color[k], capsize=5, ms=12, mfc="none")#projected datapoints
#dashed lines
#xxx = [fill(a2_aux, 100) Float64.(range(0.0, phi2_max, length=100)) fill(value(c.phih_ph), 100)]
#yyy = ff(xxx, upar)
#plot(xxx[:,2], value.(yyy), ls="--", color=color[k])#dashed lines
end
ylabel(ylab, size=22)
xlabel(L"$\phi_2$", size=22)
legend(ncol=2, fontsize=16)
#ylim(3.02, 3.18)
tight_layout()
grid(ls="dashed")
if !isnothing(p)
t = "fit_fds_ratio_" * c.info * ".pdf"
savefig(joinpath(p, t))
close()
else
display(gcf())
close("all")
end
end
##
#=========== FIT DECAY CONSTANT ==============#
## FIT DEC. CONSTS fD + fDs (LIN)
for i = 1:3
f1 = f_cat[1][i]
f2 = f_cat[2][i]
for j = 1:length(f_f_lin1)
for k = 1:length(f_f_lin1[j])
println(" j is :", j, " k is equal: ", k)
par, chi_vs_chi2exp = fit_routine([f_f_lin1[j][k], f_f_lin2[j][k]], [value.(f1.x_tofit), value.(f2.x_tofit)],
[f1.obs, f2.obs], f_n_param_comb[j], correlated_fit=false)
#aux1 = par[1] + par[2] * phi2_ph + par[3] / sqrt(f1.phih_ph)
#aux2 = par[1] + par[2] * (3*phi2_sym - 2*phi2_ph) + par[3] / sqrt(f1.phih_ph)
x = [0.0 phi2_ph f1.phih_ph]
aux1 = f_f_aux1[1][1](x, par)[1]
aux2 = f_f_aux2[1][1](x, par)[1]
doff = sum(length.([f1.obs, f2.obs])) - f_n_param_comb[j]
push!(f1.fit_res, aux1)
push!(f2.fit_res, aux2)
push!(f1.fit_param, par)
push!(f2.fit_param, par)
push!(f1.chi2_vs_exp, chi_vs_chi2exp)
push!(f2.chi2_vs_exp, chi_vs_chi2exp)
push!(f1.n_param, f_n_param_comb[j])
push!(f2.n_param, f_n_param_comb[j])
push!(f1.dof, doff)
push!(f2.dof, doff)
# use first 2 lines for uncorrelated fits
AIC = chi_vs_chi2exp * doff + 2 * f_n_param_comb[j] + ( 2*f_n_param_comb[j]^2 + 2*f_n_param_comb[j]) / (length(f1.obs)- f_n_param_comb[j] -1)
push!(f1.aic, AIC)
push!(f2.aic, AIC)
# push!(f2.aic, chi_vs_chi2exp - 2 * doff)
# push!(f1.aic, chi_vs_chi2exp - 2 * doff)
end
end
end
## FIT DEC. CONSTS fD + fDs (LIN + CLOGS)
for i = 1:3
f1 = f_cat[1][i]
f2 = f_cat[2][i]
for j = 1:length(f_f_lin1_clog)
for k = 1:length(f_f_lin1_clog[j])
println(" j is :", j, " k is equal: ", k)
par, chi2_vs_chi2exp = fit_routine([f_f_lin1_clog[j][k], f_f_lin2_clog[j][k]], [value.(f1.x_tofit), value.(f2.x_tofit)],
[f1.obs, f2.obs], f_n_param_comb_clog[j], correlated_fit=false)
x = [0.0 phi2_ph f1.phih_ph]
aux1 = f_f_aux1_clog[1][1](x, par)[1]
aux2 = f_f_aux2_clog[1][1](x, par)[1]
doff = sum(length.([f1.obs, f2.obs])) - f_n_param_comb_clog[j]
push!(f1.fit_res, aux1)
push!(f2.fit_res, aux2)
push!(f1.fit_param, par)
push!(f2.fit_param, par)
push!(f1.chi2_vs_exp, chi2_vs_chi2exp)
push!(f2.chi2_vs_exp, chi2_vs_chi2exp)
push!(f1.n_param, f_n_param_comb_clog[j])
push!(f2.n_param, f_n_param_comb_clog[j])
push!(f1.dof, doff)
push!(f2.dof, doff)
# use first 2 lines for uncorrelated fits
push!(f1.aic, chi2_vs_chi2exp * doff + 2 * f_n_param_comb_clog[j])
push!(f2.aic, chi2_vs_chi2exp * doff + 2 * f_n_param_comb_clog[j])
# push!(f2.aic, chi2_vs_chi2exp - 2 * doff)
# push!(f1.aic, chi2_vs_chi2exp - 2 * doff)
end
end
end
##
# ffD = vcat(f_f_aux1...)
# ffDs = vcat(f_f_aux2...)
ffD = vcat(f_f_aux1_clog...)
ffDs = vcat(f_f_aux2_clog...)
# ffD = vcat([vcat(f_f_aux1...), vcat(f_f_aux1_clog...)]...) # lin+log
# ffDs = vcat([vcat(f_f_aux2...), vcat(f_f_aux2_clog...)]...) # lin+log
s = [L"$\sqrt{8t_0}f_D$", L"$\sqrt{8t_0}f_{D_s}$", L"$\sqrt{8t_0}f_{\eta^{(c)}_c}$"]
for i=1:3
plot_cl_charm( f_cat[1][i], ffD, s[1], phi2_ph, ensinfo, dec=true)#, p=path_plot)
plot_cl_charm( f_cat[2][i], ffDs, s[2], phi2_ph, ensinfo, dec=true)#, p=path_plot)
plot_chi_charm(f_cat[1][i], ffD, s[1], phi2_ph, ensinfo, dec=true)#, p=path_plot)
plot_chi_charm(f_cat[2][i], ffDs, s[2], phi2_ph, ensinfo, dec=true)#, p=path_plot)
plot_cat(f_cat[1][i], s[1] , save=false)
plot_cat(f_cat[2][i], s[2] , save=false)
#plot_chi2_vs_exp(f_cat[1][i], label_cutoff)# , p=path_plot)
#plot_chi2_vs_exp(f_cat[2][i], label_cutoff)#, p=path_plot)
# plot_chi2_vs_exp_clog(f_cat[1][i], label_cutoff , p=path_plot)
# plot_chi2_vs_exp_clog(f_cat[2][i], label_cutoff , p=path_plot)
end
close("all")
##
plot_cl_charm(f_cat[1][1], vcat(f_f_aux1...),L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, dec=true, n=1 ) # 35 nice
plot_chi_charm(f_cat[1][2], vcat(f_f_aux1...), L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, dec=true, n=1) #, p=path_plot)
#############
## RESULTS ##
#############
#t0 [2.86, 3.659, 5.164, 8.595]
wpmm2 = Dict{String, Vector{Float64}}()
wpmm2["H101"] = [-1.0, 1.5, -1.0, -1.0]
wpmm2["H102r002"] = [-1.0, 2.0, -1.0, -1.0]
wpmm2["H400"] = [-1.0, 4.0, -1.0, -1.0]
wpmm2["N202"] = [-1.0, 1.0, -1.0, 14.0*5.146]
wpmm2["N200"] = [-1.0, 2.0, -1.0, -1.0]
wpmm2["N203"] = [-1.0, 2.0, -1.0, -1.0]
wpmm2["N300"] = [-1.0, 1.5, -1.0, 14.0*8.595]
wpmm2["J303"] = [-1.0, 2.0, -1.0, -1.0]
#============ CHARM MASS =============#
# category av fl_ave + etac
for (k, c) in enumerate([mass_cat[1]])
o, sys = category_av_tot([c[1],c[2]])
o_p = o * hc / t0_ph
sys_p = sys * hc / value(t0_ph)
uwerr(o, wpmm2)
uwerr(o_p, wpmm2)
println(s_obs_mass[k], "\t $o +/- $sys")
println(s_obs_mass[k], "\t $o_p +- $sys_p")
println(s_obs_mass[k], details(o_p))
end
# fl-av / sp-fl-av / etac
for i = 1:3
for (k, c) in enumerate([mass_cat[1]])
o = aic_model_ave(c[i])
sys = aic_syst(c[i])
uwerr(o, wpmm2)
println(s_obs_mass[k], "\t $o +/- $sys")
end
println()
end
##
#============ DECAY CONSTANTS ==============#
for (k, c) in enumerate(f_cat[1:2])
o, sys = category_av_tot([c[1], c[2]])
o_p = o * hc / value(t0_ph)
sys_p = sys * hc / value(t0_ph)
uwerr(o, wpmm2)
uwerr(o_p, wpmm2)
#uwerr(o)
#uwerr(o_p)
println(s_obs_f[k], "\t$o +/- $sys")
println(s_obs_f[k], "\t $o_p +- $sys_p")
println(s_obs_f[k], details(o_p))
push!(store, o_p)
end
# fl-av / sp-fl-av / etac
##
for i = 1:3
for (k, c) in enumerate(f_cat[1:2])
o = aic_model_ave(c[i])
sys = aic_syst(c[i])
uwerr(o, wpmm2)
println(s_obs_f[k], "\t $o +/- $sys")
end
println()
end
##
#######################
### EXTRA PLOTS
#######################
ens_tot = getindex.(ensembles.(obs[1]),2)
color = ["orange", "darkgreen", "magenta", "navy", "yellow", "gray", "black", "lightgreen"]
for (k,ens) in enumerate(ensname)
idx = findall(x->x==ens, ens_tot)
y_muc = obs[1][idx]
x_muc = x_sfa[idx,3]
lin_model(x,p) = p[1] .+ x.*p[2]
par, chi2 = juobs.fit_routine(lin_model, value.(x_muc), y_muc, 2)
uwerr.(x_muc)
plot(value.(x_muc), lin_model(value.(x_muc), value.(par)), color=color[k])
errorbar(value.(x_muc), xerr=err.(x_muc), value.(y_muc), yerr=err.(y_muc), color=color[k], fmt="^", label="$(ens)")
end
xlabel(L"$\phi_H^{sp-fl}$")
ylabel(L"$\sqrt{8t_0}\mu_c^R$")
legend()
tight_layout()
savefig(joinpath(path_plot,"mc_heavy_dep.pdf"))
display(gcf())
close("all")
################
###############
## TEST WITH CORRELATED FITS
ydat = f_cat[1][1].obs
uwerr.(ydat)
aux_id = getindex.(ensembles.(ydat),2)
val_map = tuple.(ydat,aux_id )
##
covar_set = Vector{Matrix{Float64}}(undef, length(unique(aux_id)) )
for (k,id) in enumerate(unique(aux_id))
idx = findall(x-> x==id, getindex.(val_map,2) )
fluc_aux = mchist.(ydat[idx], fill(id, length(idx)))
fluc = mapreduce(permutedims, vcat, fluc_aux)'
println(size(fluc)[1])
covar_set[k] = fluc' * fluc ./ size(fluc)[1]^2
end
##
covar_set[1]
err(ydat[1])^2
function covar_multi_id(obs::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
isnothing(wpm) ? uwerr.(obs) : [uwerr(yaux, wpm) for yaux in obs]
id_set = getindex.(ensembles.(obs), 2) # hardcoded assuming that ens label is always the second
#covar_set = Vector{Matrix{Float64}}(undef, length(unique(id_set)) )
cov_tot = zeros(length(id_set), length(id_set))
count = 1
for (k,id) in enumerate(unique(id_set))
idx = findall(x-> x==id, id_set )
fluc_aux = mchist.(obs[idx], fill(id, length(idx)))
fluc = mapreduce(permutedims, vcat, fluc_aux)'
covar_set = fluc' * fluc ./ size(fluc)[1]^2
dim_tmp = size(covar_set,1)
cov_tot[count:dim_tmp+count-1, count:dim_tmp+count-1] = covar_set
count +=dim_tmp
end
return cov_tot
end
##
ydat= mass_cat[1][2].obs
test = covar_multi_id(ydat)
aa = inv(test)
aa = (test + test')/2
issymmetric(aa)
isposdef(aa)
evals = Diagonal(test)^(-0.5)
bb = evals * test * evals
##
imshow(bb, cmap="hot")
title("covariance")
colorbar()
display(gcf())
close("all")
##
dim = sum(size(test[i],1) for i in 1:length(test) )
cov_tot = zeros(dim, dim)
count = 1
for (i, covar) in enumerate(test)
dim_tmp = size(test[i],1)
println(dim_tmp)
cov_tot[count:dim_tmp+count-1, count:dim_tmp+count-1] = covar
count+=dim_tmp
end
##
src/gevp analysis/gevp_plotter.jl 0 → 100644
View file @ 6c956856
# WORKING PLOTS
# plot category average
function plot_cat(cat::Cat, ylab::LaTeXString ; cut::Float64=2.5, save::Bool=false)
aic = getfield(cat, :aic)
best = aic_model_ave(cat) * hc /t0_ph
syst = aic_syst(cat) * hc /t0_ph
weigths, _ = aic_weight(cat)
all_res_aux = getfield(cat, :fit_res) .* hc
all_res = [all_res_aux[i] /t0_ph for i =1:length(all_res_aux)]
uwerr.(all_res)
uwerr(best)
uwerr(syst)
ratio_w = weigths #aic ./ minimum(aic)
fig = figure(figsize=(10,5))
#cm = get_cmap("YlOrRd")
cm = get_cmap("autumn")
fill_between(collect(1:length(all_res)), value(best-syst), value(best+syst), alpha=0.5, color="royalblue", label="Systematics")
sc = scatter(collect(1:length(all_res)), value.(all_res), edgecolors="black",c= ratio_w, cmap=cm, vmin=minimum(ratio_w), vmax=maximum(ratio_w) , zorder=3)
#sc = scatter(collect(1:length(all_res)), value.(all_res), edgecolors="black" , zorder=3)
errorbar(collect(1:length(all_res)), value.(all_res), yerr=err.(all_res), fmt="none",zorder=0, capsize=2,marker="s", lw=0.8, mec="black", ecolor="black")
errorbar(-2, value(best), yerr=err(best), fmt="s",mfc="navy",ms=8, mec="navy",ecolor="navy", capsize=2, zorder=0, label="Model average")
xlabel("Models", size=20)
ylabel(ylab, size=20)
cb = colorbar(sc)
#cb[:set_label](string("AIC", L"$(M_i)/$", "min(AIC)"))
cb[:set_label]( "AIC weights")
legend()
#ylim(0.46,0.58)
display(fig)
if save
t = string("cat_ave_", cat.info, ".pdf")
savefig(joinpath(path_plot,t))
end
end
function plot_cat_tot(cat::Vector{Cat}, ylab::LaTeXString ; cut::Float64=2.5, save::Bool=false)
best, syst = category_av_tot(cat)
println(best)
best *=hc/t0_ph
best -= 2.
syst *=hc/t0_ph
println(syst)
weigths = aic_weight.(cat)
weigths = vcat(weigths[1][1], weigths[2][1])
all_res_aux = vcat(getfield.(cat, :fit_res) .* hc...)
all_res = [all_res_aux[i] /t0_ph for i =1:length(all_res_aux)]
all_res .-=2.
uwerr.(all_res)
uwerr(best)
uwerr(syst)
ratio_w = weigths #aic ./ minimum(aic)
fig = figure(figsize=(10,5))
cm = get_cmap("YlGn")
#cm = get_cmap("autumn")
fill_between(collect(1:length(all_res)), value(best-syst), value(best+syst), alpha=0.5, color="royalblue", label="Systematics")
sc = scatter(collect(1:length(all_res)), value.(all_res), edgecolors="black",c= ratio_w, cmap=cm, vmin=minimum(ratio_w), vmax=maximum(ratio_w) , zorder=3)
#sc = scatter(collect(1:length(all_res)), value.(all_res), edgecolors="black" , zorder=3)
errorbar(collect(1:length(all_res)), value.(all_res), yerr=err.(all_res), fmt="none",zorder=0, capsize=2,marker="s", lw=0.8, mec="black", ecolor="black")
errorbar(-2, value(best), yerr=err(best), fmt="s",mfc="navy",ms=8, mec="navy",ecolor="navy", capsize=2, zorder=0, label="Model average")
xlabel("Models", size=20)
ylabel(ylab, size=20)
axvline(64, ls="--", color="red")
cb = colorbar(sc)
#cb[:set_label](string("AIC", L"$(M_i)/$", "min(AIC)"))
cb[:set_label]( "IC weights")
legend()
#ylim(0.46,0.58)
display(fig)
if save
t = string("cat_ave_tot", ".pdf")
savefig(joinpath(path_plot,t))
end
end
# plot continuum limit for a given category
function plot_cl_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph::uwreal, ensinfo::Vector{EnsInfo}; mrat::uwreal=uwreal(1.), dec::Bool=false, p::Union{String, Nothing}=nothing,
wpm::Union{Dict{String,Vector{Float64}}, Nothing}=nothing, n::Union{Nothing, Int64}=nothing)
#Format
a2_max = 0.05
color = ["orange", "darkgreen", "magenta", "navy"]
fmt = ["^", "v", "<", ">"]
a_sp = [0.087, 0.077, 0.065, 0.05]
#Choose model
if isnothing(n) #n is nothing -> model with lowest AIC
aic, n = findmin(getfield(c, :aic))
println(n)
else
aic = 0.0
end
#Model info
phih_2_use = value(c.phih_ph)
par_n = getfield(c, :n_param)[n]
upar = getfield(c, :fit_param)[n]
if !dec
res = value(mrat) * (upar[1] + upar[2]*phi2_ph + upar[3]*phih_2_use)
else
res = getfield(c, :fit_res)[n]
end
func = ff[n]
println("Results corresponding to the best fit: ", res)
println("AIC corresponding to the best fit ", aic)
println("chi2/chiexp corresponding to the best fit ", c.chi2_vs_exp[n] )
##PLOTS
#xx, yy -> Continuum limit band
xx = [Float64.(range(1e-8, a2_max, length=100)) fill(phi2_ph, 100) fill(phih_2_use,100)]
#xx = [Float64.(range(0.0, a2_max, length=100)) fill(phi2[1], 100) fill(phih_2_use,100)]
yy = func(xx, upar) * value(mrat)
#Project one value of muc_i (i=1,2,3)
n_unique = unique(i-> c.x_tofit[i,2], 1:length(c.x_tofit[:,2]))
p2 = c.x_tofit[n_unique, 2]
a2 = c.x_tofit[n_unique, 1]
pc = c.x_tofit[n_unique, 3]
x = [a2 p2 pc]
x2 = [a2 fill(phi2_ph, length(a2)) fill(phih_2_use, length(p2))] #physical value phih
#x2 = [a2 fill(phi2[1], length(a2)) fill(phih_2_use, length(p2))] #project to symmetric line only
y = ( c.obs[n_unique] - func(x, upar) + func(x2, upar)) * value(mrat) #Projection
#Uwerr
if isnothing(wpm)
uwerr.(yy)
uwerr(res)
uwerr.(y)
else
[uwerr(yy_aux, wpm) for yy_aux in yy]
uwerr(res, wpm)
[uwerr(y_aux, wpm) for y_aux in y]
end
#Plotting
figure(figsize=(10,6))
fill_between(xx[:,1], value.(yy)+ err.(yy), value.(yy)-err.(yy), alpha=0.3, color="royalblue") #CL extrapolation (band)
errorbar(0.0, value(res), err(res), fmt="s", color="red", label="CL", capsize=5, ms=12)#CL extrapolation
for (k,b) in enumerate(sort(unique(getfield.(ensinfo, :beta))))
n_ = findall(x->x.beta == b, ensinfo)
a2_aux = mean(value.(a2[n_]))
errorbar(value.(x[n_,1]), value.(y[n_]), err.(y[n_]), fmt=fmt[k], label=string(L"$a \approx$", a_sp[k], " fm"), color=color[k], capsize=5, ms=12, mfc="none")#projected datapoints
end
ylabel(ylab, size=20)
xlabel(L"$a^2/8t_0$", size=20)
legend(ncol=2, fontsize=16)
grid(ls="dashed")
tight_layout()
#ylim(0.46, 0.58)
if !isnothing(p)
t = "fit_cl_" * c.info * ".pdf"
savefig(joinpath(p, t))
close()
else
display(gcf())
close("all")
end
end
# plot continuum limit for all categories
function plot_all_cl_charm(cvec::Vector{Cat}, ff::Vector{Function}, ylab::LaTeXString, phi2_ph::uwreal, ensinfo::Vector{EnsInfo}; mrat::uwreal=uwreal(1.), dec::Bool=false, p::Union{String, Nothing}=nothing,wpm::Union{Dict{String,Vector{Float64}}, Nothing}=nothing, n::Union{Nothing, Vector{Int64}}=nothing)
#Format
a2_max = 0.05
color = ["orange", "navy", "darkgreen", "magenta", "navy"]
fmt = ["^", "v", "<", ">"]
a_sp = [0.087, 0.077, 0.065, 0.05]
labels = [L"$\phi_H = \sqrt{8t_0}m_{\eta_c}$", L"$\phi_H = \sqrt{8t_0}m_{\bar D}$", L"$\phi_H = \sqrt{8t_0}m_{\mathrm{sp-fl}}$" ]
figure(figsize=(10,6))
for (idx,c) in enumerate(cvec)
#Choose model
if n[idx] == 0 #n is nothing -> model with lowest AIC
aic, n_best = findmin(getfield(c, :aic))
else
n_best = n[idx]
aic = 0.0
end
phih_2_use = value(c.phih_ph)
par_n = getfield(c, :n_param)[n_best]
upar = getfield(c, :fit_param)[n_best]
if !dec
res = value(mrat) * (upar[1] + upar[2]*phi2_ph + upar[3]*phih_2_use)
else
@info("Here for dec const the best res is not consistent wiht chiral logs")
res = upar[1] + upar[2]*phi2_ph + upar[3] / sqrt(phih_2_use)
end
func = ff[n_best]
println("Category: ", c.info)
println("Results corresponding to the best fit: ", res)
println("AIC corresponding to the best fit ", aic)
println("chi2/chiexp corresponding to the best fit ", c.chi2_vs_exp[n_best] )
##PLOTS
#xx, yy -> Continuum limit band
xx = [Float64.(range(1e-8, a2_max, length=100)) fill(phi2_ph, 100) fill(phih_2_use,100)]
yy = func(xx, upar) * value(mrat)
#Project one value of muc_i (i=1,2,3)
n_unique = unique(i-> c.x_tofit[i,2], 1:length(c.x_tofit[:,2]))
p2 = c.x_tofit[n_unique, 2]
a2 = c.x_tofit[n_unique, 1]
pc = c.x_tofit[n_unique, 3]
x = [a2 p2 pc]
# before it was
# x2 = [a2 p2 fill(phih_2_use, length(p2))] #physical value phih
# but to me should be
x2 = [a2 fill(phi2_ph, length(p2)) fill(phih_2_use, length(p2))] #physical value phih
y = ( c.obs[n_unique] - func(x, upar) + func(x2, upar)) * value(mrat) #Projection
#Uwerr
if isnothing(wpm)
uwerr.(yy)
uwerr(res)
uwerr.(y)
else
[uwerr(yy_aux, wpm) for yy_aux in yy]
uwerr(res, wpm)
[uwerr(y_aux, wpm) for y_aux in y]
end
#Plotting
fill_between(xx[:,1], value.(yy)+ err.(yy), value.(yy)-err.(yy), alpha=0.2, color=color[idx]) #CL extrapolation (band)
#errorbar(0.0 +(idx-1)/800, value(res), err(res), fmt="s", color=color[idx], capsize=5, ms=12)#CL extrapolation
errorbar(value.(x[:,1]), value.(y[:]), err.(y[:]), fmt=fmt[idx], color=color[idx], capsize=5, ms=12, label = labels[idx], mfc="none")#projected datapoints
end
#tight_layout()
xlim(-0.001, 0.05)
#ylim(2.3,3.5)
grid(ls="dashed")
xlabel(L"$a^2/8t_0$", size=20)
ylabel(ylab, size=20)
legend(fontsize=16, loc="upper left")
if !isnothing(p)
t = "mc_cl_all_cat.pdf"
savefig(joinpath(p, t))
close()
else
display(gcf())
close("all")
end
end
# plot chiral extrapolation for a given category
function plot_chi_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph::uwreal, ensinfo::Vector{EnsInfo}; p::Union{String, Nothing}=nothing,
mrat::uwreal=uwreal(1.), dec::Bool=false,wpm::Union{Dict{String,Vector{Float64}}, Nothing}=nothing, n::Union{Nothing, Int64}=nothing)
#Format
phi2_max = maximum(value.(phi2))
color = ["orange", "darkgreen", "magenta", "navy"]
fmt = ["^", "v", "<", ">"]
a_sp = [0.087, 0.077, 0.065, 0.05]
#Choose model
if isnothing(n) #n is nothing -> model with lowest AIC
aic, n = findmin(getfield(c, :aic))
println(n)
else
aic = 0.0
end
#Model info
@info("error in t0_ph not included")
phih_2_use = c.phih_ph #value(c.phih_ph)
par_n = getfield(c, :n_param)[n]
upar = getfield(c, :fit_param)[n]
if !dec
res = value(mrat) * (upar[1] + upar[2]*phi2_ph + upar[3]*phih_2_use)
else
res = getfield(c, :fit_res)[n]
end
func = ff[n]
println("Results corresponding to the best fit: ", res)
println("AIC corresponding to the best fit ", aic)
println("chi2/chiexp corresponding to the best fit ", c.chi2_vs_exp[n] )
##PLOTS
#xx, yy -> Continuum limit band
xx = [fill(1e-8, 100) Float64.(range(0.01, phi2_max, length=100)) fill(phih_2_use,100)]
yy = func(xx, upar) * value(mrat)
#Project one value of muc_i (i=1,2,3)
n_unique = unique(i-> c.x_tofit[i,2], 1:length(c.x_tofit[:,2]))
p2 = c.x_tofit[n_unique, 2]
a2 = c.x_tofit[n_unique, 1]
pc = c.x_tofit[n_unique, 3]
x = [a2 p2 pc]
x2 = [a2 p2 fill(phih_2_use, length(p2))] #physical value phih
y = ( c.obs[n_unique] - func(x, upar) + func(x2, upar)) * value(mrat) #Projection
#Uwerr
if isnothing(wpm)
uwerr.(yy)
uwerr(res)
uwerr.(y)
else
[uwerr(yy_aux, wpm) for yy_aux in yy]
uwerr(res, wpm)
[uwerr(y_aux, wpm) for y_aux in y]
end
#Plotting
figure(figsize=(10,6))
fill_between(xx[:,2], value.(yy)+ err.(yy), value.(yy)-err.(yy), alpha=0.3, color="royalblue") #CL extrapolation (band)
errorbar(value(phi2_ph), value(res), err(res), fmt="s", color="red", label=L"$\phi_2 = \phi^\mathrm{phys}_2$", capsize=5, ms=12)#CL extrapolation
for (k,b) in enumerate(sort(unique(getfield.(ensinfo, :beta))))
n_ = findall(x->x.beta == b, ensinfo)
a2_aux = mean(value.(a2[n_]))
errorbar(value.(x[n_,2]), value.(y[n_]), err.(y[n_]), fmt=fmt[k], label=string(L"$a \approx$", a_sp[k], " fm"), color=color[k], capsize=5, ms=12, mfc="none")#projected datapoints
#dashed lines
xxx = [fill(a2_aux, 100) Float64.(range(0.0, phi2_max, length=100)) fill(value(c.phih_ph), 100)]
yyy = func(xxx, upar) * value(mrat)
plot(xxx[:,2], value.(yyy), ls="--", color=color[k])#dashed lines
end
# fill between physical values
am = M[1] * t0_ph / hc
uwerr(am)
fill_between(xx[:,2], value(am)+err(am), value(am)-err(am), color="forestgreen", alpha=0.2, zorder=1, label="PDG")
ylabel(ylab, size=22)
xlabel(L"$\phi_2$", size=22)
legend(ncol=2, fontsize=16)
ylim(3.85, 4.2)
tight_layout()
grid(ls="dashed")
if !isnothing(p)
t = "fit_phi2_" * c.info * ".pdf"
savefig(joinpath(p, t))
close()
else
display(gcf())
close("all")
end
end
# plot chi2_vs_exp adn weights for a given category
function plot_chi2_vs_exp(c::Cat, label::Vector{Vector{String}}; p::Union{String, Nothing}=nothing)
chi2_exp = getfield(c, :chi2_vs_exp)
label = vcat(label, label[2:end])
figure(figsize=(12,8))
label_plot = []
for i in 1:length(label)
aux = ""
[aux *= label[i][k] * " " for k in 1:length(label[i])]
push!(label_plot, aux)
end
xx = collect(1:length(chi2_exp))
subplots_adjust(hspace=0.0)
subplot(311)
ax1 = gca()
setp(ax1.get_xticklabels(),visible=false) # Disable x tick labels
plot(xx, aic_weight(c)[1], marker="^", ls="dashed", color="forestgreen")
ylabel(L"$w^\mathrm{AIC}$")
subplot(312)
plot(xx, chi2_exp, marker="^", ls="dashed", color="royalblue")
xticks(xx, label_plot, rotation=90)
ylabel(L"$\chi^2 / \chi^2_{\mathrm{exp}}$")
if !isnothing(p)
t = "chi2_vs_exp_" * c.info * ".pdf"
savefig(joinpath(p, t))
close()
else
display(gcf())
close("all")
end
end
function plot_chi2_vs_exp_clog(c::Cat,label::Vector{Vector{String}}; p::Union{String, Nothing}=nothing)
chi2_exp = getfield(c, :chi2_vs_exp)
label_log = deepcopy(label)
label_log = vcat.("LOG", label_log)
#label = vcat(label, label_log)
figure(figsize=(12,8))
label_plot = []
for i in 1:length(label)
aux = ""
[aux *= label[i][k] * " " for k in 1:length(label[i])]
push!(label_plot, aux)
end
xx = collect(1:length(chi2_exp))
subplots_adjust(hspace=0.0)
subplot(311)
ax1 = gca()
setp(ax1.get_xticklabels(),visible=false) # Disable x tick labels
plot(xx, aic_weight(c)[1], marker="^", ls="dashed", color="forestgreen")
ylabel(L"$w^\mathrm{AIC}$")
subplot(312)
plot(xx, chi2_exp, marker="^", ls="dashed", color="royalblue")
xticks(xx, label_plot, rotation=90)
ylabel(L"$\chi^2 / \chi^2_{\mathrm{exp}}$")
if !isnothing(p)
t = "chi2_vs_exp_" * c.info * ".pdf"
savefig(joinpath(p, t))
close()
else
display(gcf())
close("all")
end
end
##
# plot histogram for a vector of categories
function plot_hist_tot(c::Vector{Cat}, lbl::LaTeXString; p::Union{Nothing, String}=nothing)
aux_aic = getfield.(c, :aic)
aux_fit_res = getfield.(c, :fit_res)
aic = aux_aic[1]
fit_res = aux_fit_res[1]
for i = 2:length(aux_aic)
aic = vcat(aic, aux_aic[i])
fit_res = vcat(fit_res, aux_fit_res[i])
end
best = aic_model_ave(aic, fit_res) *hc /t0_ph
sys = aic_syst(aic, fit_res) * hc / t0_ph
sys = value(sys)
println(best)
w, _ = aic_weight(aic)
all_res_aux = vcat(getfield.(c, :fit_res) .* hc... )
all_res = [all_res_aux[i] /t0_ph for i =1:length(all_res_aux)]
println(length(all_res))
uwerr.(all_res)
uwerr(best)
fig = figure()
all_res .-=1.9
hist(value.(all_res), bins=15, histtype="stepfilled",alpha=0.5, ec="k", color="navy", weights=w)
axvline(value(best) -1.9, 0, 1, color="orange", label="model average")
ylim = plt.ylim()
fill_betweenx(ylim, value(best)-sys-1.9,value(best)+sys-1.9, color="orange", alpha=0.4)
ylabel("Frequency")
xlabel(lbl)
legend()
tight_layout()
plt.ylim(ylim)
if !isnothing(p)
o = split(c[1].info, ", ")[1]
t = string("hist_", o, ".pdf")
savefig(joinpath(p, t))
close()
else
display(gcf())
end
end
####
##########################
#PLOTTER FOR MESON MASSES
##########################
function plot_chi_masses(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph::uwreal, ensinfo::Vector{EnsInfo}; p::Union{String, Nothing}=nothing,
wpm::Union{Dict{String,Vector{Float64}}, Nothing}=nothing, n::Union{Nothing, Int64}=nothing)
#Format
phi2_max = 0.744666
color = ["orange", "darkgreen", "magenta", "navy"]
fmt = ["^", "v", "<", ">"]
a_sp = [0.087, 0.077, 0.065, 0.05]
#Choose model
if isnothing(n) #n is nothing -> model with lowest AIC
aic, n = findmin(getfield(c, :aic))
end
#Model info
phih_2_use = c.phih_ph
par_n = getfield(c, :n_param)[n]
upar = getfield(c, :fit_param)[n]
res = upar[1] + upar[2]*phi2_ph + upar[3]*phih_2_use
func = ff[n]
##PLOTS
#xx, yy -> Continuum limit band
xx = [fill(0.0, 100) Float64.(range(0.01, phi2_max, length=100)) fill(phih_2_use,100)]
yy = func(xx, upar)
#Project one value of muc_i (i=1,2,3)
n_unique = unique(i-> c.x_tofit[i,2], 1:length(c.x_tofit[:,2]))
p2 = c.x_tofit[n_unique, 2]
a2 = c.x_tofit[n_unique, 1]
pc = c.x_tofit[n_unique, 3]
x = [a2 p2 pc]
x2 = [a2 p2 fill(phih_2_use, length(p2))] #physical value phih
y = c.obs[n_unique] - func(x, upar) + func(x2, upar) #Projection
#Uwerr
if isnothing(wpm)
uwerr.(yy)
uwerr(res)
uwerr.(y)
else
[uwerr(yy_aux, wpm) for yy_aux in yy]
uwerr(res, wpm)
[uwerr(y_aux, wpm) for y_aux in y]
end
#Plotting
figure(figsize=(10,6))
fill_between(xx[:,2], value.(yy)+ err.(yy), value.(yy)-err.(yy), alpha=0.3, color="royalblue") #CL extrapolation (band)
errorbar(value(phi2_ph), value(res), err(res), fmt="s", color="red", label=L"$\phi_2 = \phi^\mathrm{phys}_2$", capsize=5, ms=12)#CL extrapolation
for (k,b) in enumerate(sort(unique(getfield.(ensinfo, :beta))))
n_ = findall(x->x.beta == b, ensinfo)
a2_aux = mean(value.(a2[n_]))
errorbar(value.(x[n_,2]), value.(y[n_]), err.(y[n_]), fmt=fmt[k], label=string(L"$a \approx$", a_sp[k], " fm"), color=color[k], capsize=5, ms=12, mfc="none")#projected datapoints
#dashed lines
xxx = [fill(a2_aux, 100) Float64.(range(0.0, phi2_max, length=100)) fill(value(c.phih_ph), 100)]
yyy = func(xxx, value.(upar))
plot(xxx[:,2], yyy, ls="--", color=color[k])#dashed lines
end
ylabel(ylab, size=20)
xlabel(L"$\phi_2$", size=20)
#ylim(3.9, 4.2)
grid(ls="dashed")
legend(ncol=2, fontsize=16)
tight_layout()
if !isnothing(p)
t = "fit_phi2_" * c.info * ".pdf"
savefig(joinpath(p, t))
close()
else
display(gcf())
close("all")
end
end
function plot_hist(cat::Vector{Cat} ; save::Bool=false)
best = category_av(cat)[1] *hc /t0_ph[1]
println(best)
all_res_aux = vcat(getfield.(cat, :fit_res) .* hc... )
all_res = [all_res_aux[i] /t0_ph[1] for i =1:length(all_res_aux)]
println(length(all_res))
uwerr.(all_res)
uwerr(best)
weigths = vcat(aic_weigth.(cat)...)
println(size(weigths))
fig = figure()
hist(value.(all_res), bins=20, density=true, weights=weigths, histtype="stepfilled", facecolor="royalblue",ec="k", alpha=0.3)
fill_betweenx(collect(0:5), value(best)-err(best)+4,value(best)+err(best)-4, color="tomato", alpha=0.4)
axvline(value(best),0, 5, color="tomato", label="best estimate" )
xlim(1460, 1520)
ylim(0,0.1)
ylabel("Probability")
xlabel(L"$M_c^{\rm{RGI}}\ [\rm{MeV}]$")
legend()
display(fig)
if save
path_plot= "/Users/ale/Desktop/plots"
t = string("hist_muc.pdf")
savefig(joinpath(path_plot,t))
end
close()
end
####
## PLOT FROM JAVIER TO TEST decays
###
src/gevp analysis/gevp_tm.jl 0 → 100644
View file @ 6c956856
#################################
# This file was created by Alessandro Conigli - 2022
# The pourpose of this code is to compute twisted mass
# observables from raw correlators stored in mesons.dat files.
# The analysis is performed at the level of a single ensemble
# and results are then stored in BDIO files
# It extracts:
# - charm quark mass
# - ll hl hh pseudo and vector masses
# - ll hl hh pseudo and vector leptonic decays
#
# The analysis is performed by solving a GEVP
##################################
using Base: String
using LaTeXStrings: length
using OrderedCollections
using juobs, BDIO, DelimitedFiles, ADerrors, LaTeXStrings, PyPlot
#============= SET UP VARIABLES ===========#
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["text.usetex"] = false
rcParams["mathtext.fontset"] = "cm"
rcParams["font.size"] =10
rcParams["axes.labelsize"] =22
rcParams["axes.titlesize"] = 18
plt.rc("text", usetex=false) # set to true once you install latex
#============ PATHS & INFO =================#
# path to data
const path_data_w = "/Users/alessandroconigli/Lattice/data/aux_obs_data/wilson"
const path_data_tm = "/Users/alessandroconigli/Lattice/data/charm_full_Dirac"
# path to plateau
const path_plat_w = "/Users/alessandroconigli/Lattice/data/aux_obs_data/plat_wilson.txt"
const path_plat_tm_mps = "/Users/alessandroconigli/Lattice/gevp-automated-analysis/plat/plat_tm_mps.txt"
const path_plat_tm_mvec = "/Users/alessandroconigli/Lattice/gevp-automated-analysis/plat/plat_tm_mvec.txt"
const path_plat_tm_dps = "/Users/alessandroconigli/Lattice/gevp-automated-analysis/plat/plat_tm_dps.txt"
const path_plat_tm_dvec = "/Users/alessandroconigli/Lattice/gevp-automated-analysis/plat/plat_tm_dvec.txt"
# path to aux obs
const path_rw = "/Users/alessandroconigli/Lattice/data/aux_obs_data/rwf"
const path_md = "/Users/alessandroconigli/Lattice/data/aux_obs_data/md"
# path to bdio
const path_bdio_md = "/Users/alessandroconigli/Lattice/data/bdio_charm/md"
const path_bdio_tm = "/Users/alessandroconigli/Lattice/data/bdio_charm/tm"
const path_bdio_tm_shifted = "/Users/alessandroconigli/Lattice/data/bdio_charm/tm_shifted"
const path_plot = "/Users/alessandroconigli/MyDrive/phd/analysis/charm_gevp_2022/plots"
# ensembles to analyse
const ens_list = ["H102r002"]
# reweighting and mass shift flags
const rwf = true
const mass_shift = true
const tau = 3 # gevp shift parameter
const TSM = false # set whether TSM is used or not
@warning("\nTSM FLAG MODE: FALSE\n")
#=============== INCLUDES ===============#
include("types.jl")
include("tools.jl")
include("const.jl")
include("read_bdio.jl")
#=============== ENSEMBLE INFO FROM DATABASE ==================#
ensinfo = Vector{EnsInfo}(undef, length(ens_list))
for i in 1:length(ens_list)
ens = ens_list[i]
try
ensinfo[i]= EnsInfo(ens, ens_db[ens], trunc_db[ens])
catch
error("The ensemble id ", ens, " was not found in the const.jl ens_db database.
Please check the ensemble id or update the database")
end
end
wpmm = Dict{String, Vector{Float64}}()
wpmm["H101"] = [-1.0, 4.0, -1.0, 14.0*2.86]
wpmm["H102r002"] = [-1.0, 4.0, -1.0, 14.0*2.86]
wpmm["H105"] = [-1.0, 4.0, -1.0, 14.0*2.86]
wpmm["H400"] = [-1.0, 4.0, -1.0, 14.0*3.659]
wpmm["N202"] = [-1.0, 4.0, -1.0, 14.0*5.164]
wpmm["N200"] = [-1.0, 4.0, -1.0, 14.0*5.164]
wpmm["N203"] = [-1.0, 4.0, -1.0, 14.0*5.164]
wpmm["D200"] = [-1.0, 4.0, -1.0, 14.0*5.164]
wpmm["N300"] = [-1.0, 4.0, -1.0, 14.0*8.595]
wpmm["J303"] = [-1.0, 4.0, -1.0, 14.0*8.595]
##
#=============== ANALYSIS ===============#
gen_mulist = Vector{Vector{Vector{Float64}}}(undef, length(ens_list))
# obs storage
obs_tm = Vector(undef, length(ensinfo))
if mass_shift
deltam = Vector(undef, length(ensinfo))
obs_tm_shifted = Vector(undef, length(ensinfo))
end
# compute observables
for (k,ens) in enumerate(ensinfo)
# t0
t0_ens, YW, W = get_t0(path_data_w, ens, rw=rwf, path_rw=path_rw)
# tm correlators
pp_tm, pp_tmW, _ = get_corr(path_data_tm, ens,"G5", "G5", rw=rwf, path_rw=path_rw, tsm=TSM)
a1a1_tm, _, _ = get_corr(path_data_tm, ens, "G1G5", "G1G5", rw=rwf, path_rw=path_rw, tsm=TSM)
a2a2_tm, _, _ = get_corr(path_data_tm, ens, "G2G5", "G2G5", rw=rwf, path_rw=path_rw, tsm=TSM)
a3a3_tm, _, _ = get_corr(path_data_tm, ens, "G3G5", "G3G5", rw=rwf, path_rw=path_rw, tsm=TSM)
pa1_tm, _, _ = get_corr(path_data_tm, ens, "G5", "G1G5", rw=rwf, path_rw=path_rw, tsm=TSM)
pa2_tm, _, _ = get_corr(path_data_tm, ens, "G5", "G2G5", rw=rwf, path_rw=path_rw, tsm=TSM)
pa3_tm, _, _ = get_corr(path_data_tm, ens, "G5", "G3G5", rw=rwf, path_rw=path_rw, tsm=TSM)
mu_list = gen_mulist[k] = getfield.(pp_tm, :mu)
# gevp matrices for ps and vec masses
mat_tm1 = comp_mat_multigamma(ens, pp_tm, pa1_tm, a1a1_tm)
mat_tm2 = comp_mat_multigamma(ens, pp_tm, pa2_tm, a2a2_tm)
mat_tm3 = comp_mat_multigamma(ens, pp_tm, pa3_tm, a3a3_tm)
# gevp matrices for ps and vec decays
mat_dec_ps = comp_mat(ens, pp_tm, tau)
mat_dec_vec1 = comp_mat(ens, a1a1_tm, tau)
mat_dec_vec2 = comp_mat(ens, a2a2_tm, tau)
mat_dec_vec3 = comp_mat(ens, a3a3_tm, tau)
# gevp meson masses
m_tm_ps, m_tm_vec1 = gevp_mass(mat_tm1, path_plat_tm_mps, path_plat_tm_mvec, t0=3, pl=true)
_, m_tm_vec2 = gevp_mass(mat_tm2, path_plat_tm_mps, path_plat_tm_mvec, t0=3, pl=false)
_, m_tm_vec3 = gevp_mass(mat_tm3, path_plat_tm_mps, path_plat_tm_mvec, t0=3, pl=false)
m_tm_vec = (m_tm_vec1 .+ m_tm_vec2 .+ m_tm_vec3 ) ./ 3.
# pion mass
pion_plat = select_plateau(ens, mu_list, path_plat_tm_mps)[1]
#mll_g = m_tm_ps[1]
#m_tm_ps[1] = mll_f = juobs.meff(pp_tm[1], pion_plat, pl=true)
#test = meff.(pp_tm,select_plateau(ens, mu_list, path_plat_tm_mps), pl=true)
# decay constants from gevp
# dec_ps, data = gevp_dec(mat_dec_ps, m_tm_ps, path_plat_tm_dps, t0=3, n=1, pl=true)
# dec_v1, datav1 = gevp_dec(mat_dec_vec1, m_tm_vec1, path_plat_tm_dvec, t0=3, n=1, pl=true, pseudo=false)
# dec_v2, datav2 = gevp_dec(mat_dec_vec2, m_tm_vec2, path_plat_tm_dvec, t0=3, n=1, pl=true, pseudo=false)
# dec_v3, datav3 = gevp_dec(mat_dec_vec3, m_tm_vec3, path_plat_tm_dvec, t0=3, n=1, pl=true, pseudo=false)
# dec_vec = (dec_v1 .+ dec_v2 .+ dec_v3) ./ 3.
# decay constant from ratio
dec_ps = get_f_tm(pp_tm, m_tm_ps, ens, path_plat_tm_dps, pl=true )
dec_v1 = get_fstar_tm(a1a1_tm, m_tm_vec1, ens, path_plat_tm_dvec, pl=true)
dec_v2 = get_fstar_tm(a2a2_tm, m_tm_vec2, ens, path_plat_tm_dvec, pl=false)
dec_v3 = get_fstar_tm(a3a3_tm, m_tm_vec3, ens, path_plat_tm_dvec, pl=false)
dec_vec = (dec_v1 .+ dec_v2 .+ dec_v3) ./ 3.
close("all")
#======== DECAY CONSTANT PLOTTER ===========#
#=
for i in 1:length(data)
y0s = getfield(mat_dec_ps[i], :y0)
aux_data_ps = data[i] .* [ sqrt(2 / m_tm_ps[i]^3) * sum(mu_list[i]) for k in 1:length(data[i]) ]
aux_data_vec = datav1[i] .* [ za(ensinfo[k].beta) * sqrt(2 / m_tm_vec1[i]) for j in 1:length(data[i]) ]
uwerr.(aux_data_ps)
uwerr.(aux_data_vec)
dec_v1[i] *= za(ensinfo[k].beta)
uwerr(dec_ps[i])
uwerr(dec_v1[i])
xx = collect(64:length(aux_data_ps)+63)
v , e = value(dec_ps[i]), err(dec_ps[i])
v_vec, e_vec = value(dec_v1[i]), err(dec_v1[i])
plat_ps = select_plateau(ensinfo[k], mu_list, path_plat_tm_dps)[i] .- y0s
plat_vec = select_plateau(ensinfo[k], mu_list, path_plat_tm_dvec)[i] .- y0s
fig = figure(figsize=(10,6))
errorbar(xx, value.(aux_data_ps), err.(aux_data_ps), fmt="s", mfc="none", label=L"$f_{D_{(s)}}$")
errorbar(xx[1:38], value.(aux_data_vec)[1:38], err.(aux_data_vec)[1:38], fmt="s", mfc="none", label=L"$f_{D_{(s)}}^*$")
fill_between(plat_ps[1]+62:plat_ps[2]+62, v-e, v+e, alpha=0.6 )
fill_between(plat_vec[1]+62:plat_vec[2]+62, v_vec-e_vec, v_vec+e_vec, alpha=0.6 )
xlabel(L"$x_0/a$")
ylabel(L"$af$")
xlim(65, 110)
title("$(ensembles(dec_ps[i])[1]) " * L"\mu_1="*"$(mu_list[i][1]) " *L" \mu_2="*"$(mu_list[i][2]) " )
ylim(v*0.95, v_vec*1.1)
legend(fontsize=20)
display(fig)
t = "$(ensembles(dec_ps[i])[1])_dec_"*"mu_1=$(mu_list[i][1]) " *"_mu_2=$(mu_list[i][2]) "*".pdf"
savefig(joinpath(path_plot,t))
close()
end
=#
# divide different mass sectors for storage
# ps masses
mll_tm = get_ll(mu_list, m_tm_ps, ens.deg) # ll
mls_tm = ens.deg ? mll_tm : get_ls(mu_list, m_tm_ps, ens.deg) # ls
mlh_tm = get_lh(mu_list, m_tm_ps, ens.deg) # lh
msh_tm = ens.deg ? mlh_tm : get_sh(mu_list, m_tm_ps, ens.deg) # sh
mhh_tm = get_hh(mu_list, m_tm_ps, ens.deg) # hh
# vec masses
mlh_star_tm = get_lh(mu_list, m_tm_vec, ens.deg)
msh_star_tm = ens.deg ? mlh_star_tm : get_sh(mu_list, m_tm_vec, ens.deg)
mhh_star_tm = get_hh(mu_list, m_tm_vec, ens.deg)
# ps decays
flh_tm = get_lh(mu_list, dec_ps, ens.deg)
fsh_tm = ens.deg ? flh_tm : get_sh(mu_list, dec_ps, ens.deg)
fhh_tm = get_hh(mu_list, dec_ps, ens.deg)
# vec decays
flh_star_tm = get_lh(mu_list, dec_vec, ens.deg)
fsh_star_tm = ens.deg ? flh_tm : get_sh(mu_list, dec_vec, ens.deg)
fhh_star_tm = get_hh(mu_list, dec_vec, ens.deg)
mul, mus, muh = get_mu(mu_list, ens.deg)
t0_aux = sqrt(8 * t0_ens)
t0vec_aux = fill(t0_aux, length(muh))
#Dict: Masses and decay constants
obs_tm[k] = OrderedDict(
"mll" => t0_aux * mll_tm,
"mls" => t0_aux * mls_tm,
"mlh" => t0vec_aux .* mlh_tm,
"msh" => t0vec_aux .* msh_tm,
"mhh" => t0vec_aux .* mhh_tm,
"flh" => t0vec_aux .* flh_tm,
"fsh" => t0vec_aux .* fsh_tm,
"fhh" => t0vec_aux .* fhh_tm,
"mlh_star" => t0vec_aux .* mlh_star_tm,
"msh_star" => t0vec_aux .* msh_star_tm,
"mhh_star" => t0vec_aux .* mhh_star_tm,
"flh_star" => t0vec_aux .* flh_star_tm,
"fsh_star" => t0vec_aux .* fsh_star_tm,
"fhh_star" => t0vec_aux .* fhh_star_tm,
"muh" => t0vec_aux .* muh
)
# MASS SHIFT TM OBSERVABLES
if mass_shift
# read deltam from bdio
fname = joinpath(path_bdio_md, string(ens.id,".bdio"))
deltam = read_BDIO(fname, "md", "deltam")[1]
# read dS/dm
dSdm = read_pbp(path_md, ens.id)
# definitions
prim_obs = vcat(pp_tmW...)
obs_tm_shifted[k] = deepcopy(obs_tm[k])
## LOOOP OVER OBSERVABLES + MASS SHIFT
for (key, value) in obs_tm[k]
if isa(value, uwreal) # if 1 uwreal
mds1, mds2 = md_sea(value, dSdm, prim_obs, W) .+ md_sea(value, dSdm, YW, W)
obs_tm_shifted[k][key] = obs_tm[k][key] + deltam * (2 * mds1 + mds2)
elseif isa(value, Vector{uwreal}) # if vector uwreal
for (i, v) in enumerate(value)
mds1, mds2 = md_sea(v, dSdm, prim_obs, W) .+ md_sea(v, dSdm, YW, W)
obs_tm_shifted[k][key][i] = obs_tm[k][key][i] + deltam * (2 * mds1 + mds2)
end
end
end
end
end
##
#====================== WRITE BDIO =====================#
if mass_shift
#TM SHIFTED
for (k, obs) in enumerate(obs_tm_shifted) #loop over ensembles
ens = ensinfo[k].id
p = joinpath(path_bdio_tm_shifted, string(ens, ".bdio"))
fb = BDIO_open(p, "w", ens)
uinfo = 0
for (key, value) in obs #loop over observables
if isa(value, uwreal)
write_uwreal(value, fb, uinfo)
elseif isa(value, Vector{uwreal})
[write_uwreal(v, fb, uinfo) for v in value]
end
uinfo += 1
end
BDIO_close!(fb)
end
end
# TM
for (k, obs) in enumerate(obs_tm) #loop over ensembles
ens = ensinfo[k].id
p = joinpath(path_bdio_tm, string(ens, ".bdio"))
fb = BDIO_open(p, "w", ens)
uinfo = 0
for (key, value) in obs #loop over observables
if isa(value, uwreal)
write_uwreal(value, fb, uinfo)
elseif isa(value, Vector{uwreal})
[write_uwreal(v, fb, uinfo) for v in value]
end
uinfo += 1
end
BDIO_close!(fb)
end
src/gevp analysis/read_bdio.jl 0 → 100644
View file @ 6c956856
using BDIO, ADerrors
function read_BDIO(path::String, uinfo::Int64)
r = Vector{uwreal}(undef, 0)
fb = BDIO_open(path, "r")
BDIO_seek!(fb)
if BDIO_get_uinfo(fb) == uinfo
push!(r, read_uwreal(fb))
end
while BDIO_seek!(fb, 2)
if BDIO_get_uinfo(fb) == uinfo
push!(r, read_uwreal(fb))
end
end
BDIO_close!(fb)
return r
end
function read_BDIO(path::String, type::String, obs::String)
dict_w = Dict(
"mll" => 0,
"mls" => 1,
"t0" => 2,
"phi2" => 3,
"phi4" => 4
)
dict_tm = Dict(
"mll" => 0,
"mls" => 1,
"mlh" => 2,
"msh" => 3,
"mhh" => 4,
"flh" => 5,
"fsh" => 6,
"fhh" => 7,
"mlh_star" => 8,
"msh_star" => 9,
"mhh_star" => 10,
"flh_star" => 11,
"fsh_star" => 12,
"fhh_star" => 13,
"muh" => 14
)
dict_md = Dict(
"deltam" => 0,
"t0_shifted" => 1,
"phi2_shifted" => 2,
"phi4_shifted" => 3
)
dict2dict = Dict("w" => dict_w, "tm" => dict_tm, "md" => dict_md)
if !(type in keys(dict2dict))
error("Incorrect type.\ntype = $(keys(dict2dict))")
end
if !(obs in keys(dict2dict[type]))
error("Incorrect obs.\nobs = $(keys(dict2dict[type]))")
end
return read_BDIO(path, dict2dict[type][obs])
end
src/gevp analysis/tools.jl 0 → 100644
View file @ 6c956856
# read standard data no TSM
function read_data(path::String, ens::String, g1::String="G5", g2::String="G5"; legacy::Bool=false)
p = joinpath(path, ens)
rep = filter(x->occursin("mesons.dat", x), readdir(p, join=true))
if length(rep)!=0
length(rep)==1 ? (return read_mesons(rep[1], g1, g2, legacy=legacy, id=ens)) : (return read_mesons(rep, g1, g2, legacy=legacy, id=ens))
else
error("mesons.dat file not found for ensemble ", ens, " in path ", p)
end
end
# read sloppy TSM
function read_data_sloppy(path::String, ens::String, g1::String="G5", g2::String="G5"; legacy::Bool=false)
p = joinpath(path, ens)
rep = filter(x->occursin("sloppy", x), readdir(p, join=true))
if length(rep)!=0
length(rep)==1 ? (return read_mesons(rep[1], g1, g2, legacy=legacy, id=ens)) : (return read_mesons(rep, g1, g2, legacy=legacy, id=ens))
else
error("mesons.dat file not found for ensemble ", ens, " in path ", p)
end
end
# read correction TSM
function read_data_exact(path::String, ens::String, g1::String="G5", g2::String="G5"; legacy::Bool=false)
p = joinpath(path, ens)
rep = filter(x->occursin("exact", x), readdir(p, join=true))
if length(rep)!=0
length(rep)==1 ? (return read_mesons_correction(rep[1], g1, g2, legacy=legacy, id=ens)) : (return read_mesons_correction(rep, g1, g2, legacy=legacy, id=ens))
else
error("mesons.dat file not found for ensemble ", ens, " in path ", p)
end
end
function read_pbp(path::String, ens::String)
p = joinpath(path, ens)
rep = filter(x->occursin("pbp.dat", x), readdir(p, join=true))
if length(rep)!=0
return read_md.(rep)
else
error("pbp.dat file not found for ensemble ", ens, " in path ", p)
end
end
function read_rw(path::String, ens::String)
p = joinpath(path, ens)
rep = filter(x->occursin("ms1.dat", x), readdir(p, join=true))
if length(rep)!=0
try
length(rep) == 1 ? (return read_ms1(rep[1])) : (return read_ms1.(rep))
catch
length(rep) == 1 && length(rep)!=0 ? (return read_ms1(rep[1], v="1.4")) : (return read_ms1.(rep, v="1.4"))
end
else
error("ms1.dat file not found for ensemble ", ens, " in path ", p)
end
end
function read_t0(path::String, ens::String, dtr::Int64)
p = joinpath(path, ens)
rep = filter(x->occursin("ms.dat", x), readdir(p, join=true))
if length(rep)!=0
length(rep) == 1 ? (return read_ms(rep[1], dtr=dtr, id=ens)) : (return read_ms.(rep, dtr=dtr, id=ens))
else
error("ms.dat file not found for ensemble ", ens, " in path ", p)
end
end
# in this function you have to implement the tsm reader
function get_corr(path::String, ens::EnsInfo, g1::String="G5", g2::String="G5"; rw::Bool=false, legacy::Bool=false, path_rw::String="", tsm::Bool=false)
if path_rw == ""
p_rw = path
else
p_rw = path_rw
end
if !tsm
aux = read_data(path, ens.id, g1, g2, legacy=legacy)
if !isnothing(ens.trunc)
truncate_data!(aux, ens.trunc)
end
if !rw
obs = corr_obs.(aux, L=ens.L, info=true)
return (getindex.(obs, 1), getindex.(obs, 2))
else
obs = corr_obs.(aux, L=ens.L, rw=read_rw(p_rw, ens.id), info=true)
return (getindex.(obs, 1), getindex.(obs, 2), getindex.(obs, 3))
end
else
aux1 = read_data_sloppy(path, ens.id, g1, g2, legacy=legacy)
aux2 = read_data_exact(path, ens.id, g1, g2, legacy=legacy)
if !isnothing(ens.trunc)
truncate_data!(aux1, ens.trunc)
truncate_data!(aux2, ens.trunc)
end
if !rw
obs = corr_obs_TSM.(aux1, aux2, L=ens.L, info=true)
return (getindex.(obs, 1), getindex.(obs, 2))
else
obs = corr_obs_TSM.(aux1, aux2, L=ens.L, rw=read_rw(p_rw, ens.id), info=true)
return (getindex.(obs, 1), getindex.(obs, 2), getindex.(obs, 3))
end
end
end
function get_t0(path::String, ens::EnsInfo; rw::Bool=false, pl::Bool=false, path_rw::String="")
if path_rw == ""
p_rw = path
else
p_rw = path_rw
end
data = read_t0(path, ens.id, ens.dtr)
if !isnothing(ens.trunc)
truncate_data!(data, ens.trunc)
end
plat = select_plateau_t0(ens)
if !rw
t0 = comp_t0(data, plat, L=ens.L, pl=pl, info=true)
return (getindex(t0, 1), getindex(t0, 2))
else
t0 = comp_t0(data, plat, L=ens.L, pl=pl, rw=read_rw(p_rw, ens.id), info=true)
return (getindex(t0, 1), getindex(t0, 2), getindex(t0, 3))
end
end
function get_mu(mu_list::Vector{Vector{Float64}}, deg::Bool)
mu_sorted = unique(sort(minimum.(mu_list)))
mul = mu_sorted[1]
deg ? mus = 0.0 : mus = mu_sorted[2]
muh = unique(maximum.(mu_list))
muh = filter(x-> x > mul && x > mus, muh)
return mul, mus, muh
end
function get_kappa(kappa_list::Vector{Vector{Float64}}, deg::Bool)
kappa_sorted = unique(sort(maximum.(kappa_list), rev=true))
kappal = kappa_sorted[1]
deg ? kappas = 0.0 : kappas = kappa_sorted[2]
return kappal, kappas
end
function get_ll(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 mul in mu && mu[1] == mu[2] #l-l
return obs[k]
end
end
end
function get_ll_w(kappa_list, obs::Vector, deg::Bool)
kappal, kappas = get_kappa(kappa_list, deg)
for k = 1:length(kappa_list)
kappa = kappa_list[k]
if kappal in kappa && kappa[1] == kappa[2] #l-l
return obs[k]
end
end
end
function get_ls(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 mul in mu && mus in mu #l-l
return obs[k]
end
end
end
function get_ls_w(kappa_list, obs::Vector, deg::Bool)
kappal, kappas = get_kappa(kappa_list, deg)
for k = 1:length(kappa_list)
kappa = kappa_list[k]
if kappal in kappa && kappas in kappa #l-l
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)
for k = 1:length(mu_list)
mu = mu_list[k]
if mul in mu && mu[1] != mu[2] && !(mus in mu) #l-h
push!(obs_lh, obs[k])
end
end
return obs_lh
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_sh(mu_list, obs::Vector{uwreal}, deg::Bool)
mul, mus, muh = get_mu(mu_list, deg)
obs_sh = Vector{uwreal}(undef, 0)
for k = 1:length(mu_list)
mu = mu_list[k]
if mus in mu && mu[1] != mu[2] && !(mul in mu)#s-h
push!(obs_sh, obs[k])
end
end
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]
for i =1:length(muh)
if muh[i] in mu && mu[1] == mu[2]#h-h
push!(obs_hh, obs[k])
end
end
end
return obs_hh
end
function comp_mat(ensinfo::EnsInfo, tot_obs::Vector{juobs.Corr}, tau::Int64)
mat_info = Vector{MatInfo}(undef, length(tot_obs))
for i = 1:length(tot_obs)
a11 = tot_obs[i].obs[1:end-1-2*tau]
a12 = tot_obs[i].obs[1+tau:end-1-tau]
a22 = tot_obs[i].obs[1+2*tau:end-1]
aux = juobs.get_matrix([a11, a22], [a12])
mat_info[i] = MatInfo(ensinfo, aux, tot_obs[i].mu, tot_obs[i].y0)
end
return mat_info
end
function comp_mat_multigamma(ensinfo::EnsInfo,tot11::Vector{juobs.Corr}, tot12::Vector{juobs.Corr}, tot22::Vector{juobs.Corr})
mu = getfield.(tot11,:mu)
y0 = getfield.(tot11,:y0) .+1
mat_info = Vector{MatInfo}(undef, length(mu))
for i = 1:length(mu)
a11 = tot11[i].obs[y0[i]:end-1]
a12 = tot12[i].obs[y0[i]:end-1]
a22 = tot22[i].obs[y0[i]:end-1]
#res[i] = infoMatt(a11, a12, a22)
aux = juobs.get_matrix([a11, a22], [a12])
mat_info[i] = MatInfo(ensinfo, aux, mu[i], y0[i])
end
return mat_info
end
function gevp_mass(mat_obs::Vector{MatInfo}, path_plat::String; t0::Int64=2,pl::Bool=true, wpm::Union{Nothing, Dict{Int64,Vector{Float64}}}=nothing)
y0s = getfield.(mat_obs, :y0)
mu_list = getfield.(mat_obs, :mu)
plat_en0 = Vector{uwreal}(undef, length(mat_obs))
plat_en1 = Vector{uwreal}(undef, length(mat_obs))
for i in 1:length(mat_obs)
mat = getfield(mat_obs[i], :mat_list)[y0s[i]+1:end-1] #matrices to diagonalise for a given sector [i]
evals = getall_eigvals(mat, t0)
en = energies(evals, wpm=wpm) #ground and excited state energy
plat = select_plateau(mat_obs[i].ensinfo, mu_list, path_plat)[i] .- y0s[i]
println(plat)
plat_en0[i] = plat_av(en[1], plat)
if pl
uwerr(plat_en0[i])
xx = collect(1:length(en[1]))
v = value(plat_en0[i])
e = err(plat_en0[i])
errorbar(xx, value.(en[1]), err.(en[1]), fmt="s")
fill_between(plat[1]:plat[2], v-e, v+e )
ylim(e*0.9, e*1.1)
display(gcf())
close()
end
end
return plat_en0
end
function gevp_mass(mat_obs::Vector{MatInfo}, path_m_ps::String, path_m_vec::String; t0::Int64=2,pl::Bool=true, wpm::Union{Nothing, Dict{Int64,Vector{Float64}}}=nothing)
y0s = getfield.(mat_obs, :y0)
mu_list = getfield.(mat_obs, :mu)
plat_en0 = Vector{uwreal}(undef, length(mat_obs))
plat_en1 = Vector{uwreal}(undef, length(mat_obs))
for i in 1:length(mat_obs)
mat = getfield(mat_obs[i], :mat_list) #matrices to diagonalise for a given sector [i]
evals = getall_eigvals(mat, t0)
en = energies(evals, wpm=wpm) #ground and excited state energy
plat_ps = select_plateau(mat_obs[i].ensinfo, mu_list, path_m_ps)[i] .- y0s[i]
println(plat_ps)
plat_vec = select_plateau(mat_obs[i].ensinfo, mu_list, path_m_vec)[i] .- y0s[i]
println(plat_vec)
plat_en0[i] = plat_av(en[1], plat_ps)
plat_en1[i] = plat_av(en[2], plat_vec)
if pl
fig = figure(figsize=(10,6))
uwerr(plat_en0[i])
uwerr(plat_en1[i])
xx = collect(y0s[1]:length(en[1])+y0s[1]-1)
v = value(plat_en0[i])
e = err(plat_en0[i])
v_vec = value(plat_en1[i])
e_vec = err(plat_en1[i])
errorbar(xx, value.(en[1]), err.(en[1]), fmt=">", mfc="none", capsize=2, ms=12, label=L"$M_{D_{(s)}}$")
fill_between(plat_ps[1]+y0s[1]:plat_ps[2]+y0s[1], v-e, v+e, alpha=0.6 )
errorbar(xx[1:end-10], value.(en[2])[1:end-10], err.(en[2][1:end-10]), fmt="^",mfc="none", ms=12, capsize=2, label=L"$M_{D_{(s)}}^*$")
fill_between(plat_vec[1]+y0s[1]:plat_vec[2]+y0s[1], v_vec-e_vec, v_vec+e_vec, alpha=0.6 )
xlabel(L"$x_0/a$")
ylabel(L"$am_{\mathrm{eff}}$")
title("$(ensembles(plat_en0[i])[1]) " * L"\mu_1="*"$(mu_list[i][1]) " *L" \mu_2="*"$(mu_list[i][2]) " )
# xlim(65, 110)
ylim(v-25*e, v+25*e)
legend(fontsize=20)
display(gcf())
t = "$(ensembles(plat_en0[i])[1])_m_eff_"*"mu_1=$(mu_list[i][1]) " *"_mu_2=$(mu_list[i][2]) "*".pdf"
#savefig(joinpath(path_plot,t))
close()
end
end
return plat_en0, plat_en1
end
function gevp_dec(mat_obs::Vector{MatInfo}, mass::Vector{uwreal}, path_plat::String; t0::Int64=2, wpm::Union{Nothing, Dict{Int64,Vector{Float64}}}=nothing, pl::Bool=true, n::Int64=1, pseudo::Bool=true, wilson::Bool=false)
y0s = getfield.(mat_obs, :y0)
mu_list = getfield.(mat_obs, :mu)
plat_dec0 = Vector{uwreal}(undef, length(mat_obs))
data = Vector{Vector{uwreal}}(undef,0)
for i in 1:length(mat_obs)
mat = getfield(mat_obs[i], :mat_list)[y0s[i]+2:end-1] #matrices to diagonalise for a given sector [i]
evecs = getall_eig(mat, t0)
elem = mat_elem(evecs, mat, mass[i], n) #ground and excited state energy
push!(data, elem)
plat = select_plateau(mat_obs[i].ensinfo, mu_list, path_plat)[i] .- y0s[i]
try
if pseudo
plat_dec0[i] = dec(elem, plat, mass[i], mu_list[i], pl=pl, wilson=wilson)
else
plat_dec0[i] = vec_dec(elem, plat, mass[i], pl=pl)
end
catch
println("Decay constant for the ensemble ", mat_obs[i].ensinfo.id, " failed in the sector ", mu_list[i], " pseudo sector?", pseudo)
println("The associated effective mass is ", mass[i])
plat_dec0[i] = uwreal(0)
end
end
return plat_dec0, data
end
function mat_elem(evec::Vector{Matrix{uwreal}}, ct_mat::Vector{Matrix{uwreal}}, mass::uwreal, n::Int64)
t = length(evec)
res = Vector{uwreal}(undef, t)
for i in 1:t
aux = uwdot(evec[i][:,n], uwdot(ct_mat[i], evec[i][:,n]))[1]
aux2 = 1 / (aux^2)^(0.25) * exp(mass * (i)/2)
res[i] = aux2 * uwdot(evec[i][:,n], ct_mat[i][:,n])
end
return res
end
function dec(mat_elem::Vector{uwreal}, plat::Array{Int64}, mass::uwreal, mu::Vector{Float64}; pl::Bool=true, wilson::Bool=false)
aux = plat_av(mat_elem, plat)
if pl
uwerr(aux)
x = 1:length(mat_elem)
y = value.(mat_elem)
dy = err.(mat_elem)
v = value(aux)
e = err(aux)
figure()
errorbar(x, y, yerr=dy, fmt="*", color="red")
fill_between(plat[1]:plat[2], v-e, v+e, color="green")
ylim(v - 10*e, v + 10*e)
display(gcf())
close()
end
if !wilson
return sqrt(2/ (mass)^3) * sum(mu)* aux
else
return sqrt(2/mass) * aux
end
end
function vec_dec(mat_elem::Vector{uwreal}, plat::Array{Int64}, mass::uwreal; pl::Bool=true)
aux = plat_av(mat_elem, plat)
if pl
uwerr(aux)
x = 1:length(mat_elem)
y = value.(mat_elem)
dy = err.(mat_elem)
v = value(aux)
e = err(aux)
figure()
ylim(v - 10*e, v + 10*e)
errorbar(x, y, yerr=dy, fmt="*", color="red")
fill_between(plat[1]:plat[2], v-e, v+e, color="green")
display(gcf())
end
return sqrt(2 / mass) * aux #aux / mass
end
function getall_eig(a::Vector{Matrix{uwreal}}, t0; iter=10)
n = length(a)
res = Vector{Matrix{uwreal}}(undef, n)
[res[i] = uweigvecs(a[i], a[t0]) for i=1:n]
return res
end
function get_meff_tm(obs::Vector{juobs.Corr}, ens::EnsInfo; pl::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
mu_list = getfield.(obs, :mu)
plat = select_plateau_tm(ens, mu_list)
println(plat)
return meff.(obs, plat, pl=pl, wpm=wpm)
end
function get_meff_vec(obs::Vector{juobs.Corr}, ens::EnsInfo; pl::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
mu_list = getfield.(obs, :mu)
plat = select_plateau_vec(ens, mu_list)
println(plat)
return meff.(obs, plat, pl=pl, wpm=wpm)
end
function get_meff_w(obs::Vector{juobs.Corr}, ens::EnsInfo; pl::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
kappa_list = getfield.(obs, :kappa)
plat = select_plateau_w(ens, kappa_list)
println(plat)
return meff.(obs, plat, pl=pl, wpm=wpm)
end
function get_f_tm(obs::Vector{juobs.Corr}, m::Vector{uwreal}, ens::EnsInfo, path_plat::String; pl::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
mu_list = getfield.(obs, :mu)
plat = select_plateau(ens, mu_list, path_plat)
return dec_const_pcvc.(obs, plat, m, pl=pl, wpm=wpm)
end
function get_fstar_tm(obs::Vector{juobs.Corr}, m::Vector{uwreal}, ens::EnsInfo, path_plat::String; pl::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
mu_list = getfield.(obs, :mu)
plat = select_plateau(ens, mu_list, path_plat)
return dec_const_w.(obs, obs, plat, m, pl=pl, wpm=wpm)
end
function select_plateau(ensinf::EnsInfo, mu_list, path_plat)
ens =ensinf.id
deg = ensinf.deg
mul, mus, muh = get_mu(mu_list, deg)
f = readdlm(path_plat)
head = String.(f[:, 1])
delim = findall(x-> occursin("#", x ), head)
edelim = findfirst(x-> occursin(ens, x), head)
cdelim = findfirst(x-> x.== edelim, delim)
del = delim[cdelim]+1 : delim[cdelim+1]-1
#del = delim
plat = Vector{Vector{Int64}}(undef, 0)
#plat = Dict()
plat_ = Int64.(f[del, 2:3])
head_ = String.(f[del, 1])
for k = 1:length(mu_list)
mu = mu_list[k]
if mul in mu && mu[1] != mu[2] && !(mus in mu) #heavy-light
n = findfirst(x-> occursin("lh", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif mul in mu && mu[1] == mu[2] #light-light
n = findfirst(x-> occursin("ll", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif mul in mu && mus in mu#strange-light
n = findfirst(x-> occursin("ls", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif mus in mu && mu[1] != mu[2] && !(mul in mu)#heavy-strange
n = findfirst(x-> occursin("sh", x ), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif mus in mu && mu[1] == mu[2] && !(mul in mu)#strange-strange
n = findfirst(x-> occursin("ss", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif !(mul in mu) && !(mus in mu) #heavy-heavy
n = findfirst(x-> occursin("hh", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
end
end
return plat
end
function select_plateau_t0(ensinfo::EnsInfo)
ens = ensinfo.id
f = readdlm(path_plat_w)
head = String.(f[:,1])
delim = findall(x->occursin("#", x), head)
edelim = findfirst(x->occursin(ens, x), head)
cdelim = findfirst(x-> x.== edelim, delim)
del = delim[cdelim]+1 : delim[cdelim+1]-1
head_ = String.(f[del, 1])
plat_ = Int64.(f[del, 2:3])
plat = Vector{Int64}(undef, 2)
n = findfirst(x-> occursin("t0", x), head_)
plat[1] = plat_[n,1]
plat[2] = plat_[n,2]
return plat
end
function select_plateau_tm(ensinf::EnsInfo, mu_list)
ens =ensinf.id
deg = ensinf.deg
mul, mus, muh = get_mu(mu_list, deg)
f = readdlm(path_plat_tm)
head = String.(f[:, 1])
delim = findall(x-> occursin("#", x ), head)
edelim = findfirst(x-> occursin(ens, x), head)
cdelim = findfirst(x-> x.== edelim, delim)
del = delim[cdelim]+1 : delim[cdelim+1]-1
#del = delim
plat = Vector{Vector{Int64}}(undef, 0)
#plat = Dict()
plat_ = Int64.(f[del, 2:3])
head_ = String.(f[del, 1])
for k = 1:length(mu_list)
mu = mu_list[k]
if mul in mu && mu[1] != mu[2] && !(mus in mu) #heavy-light
n = findfirst(x-> occursin("lh", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif mul in mu && mu[1] == mu[2] #light-light
n = findfirst(x-> occursin("ll", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif mul in mu && mus in mu#strange-light
n = findfirst(x-> occursin("ls", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif mus in mu && mu[1] != mu[2] && !(mul in mu)#heavy-strange
n = findfirst(x-> occursin("sh", x ), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif mus in mu && mu[1] == mu[2] && !(mul in mu)#strange-strange
n = findfirst(x-> occursin("ss", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif !(mul in mu) && !(mus in mu) #heavy-heavy
n = findfirst(x-> occursin("hh", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
end
end
return plat
end
function (ensinf::EnsInfo, mu_list)
ens =ensinf.id
deg = ensinf.deg
mul, mus, muh = get_mu(mu_list, deg)
f = readdlm(path_plat_tm)
head = String.(f[:, 1])
delim = findall(x-> occursin("#", x ), head)
edelim = findfirst(x-> occursin(ens, x), head)
cdelim = findfirst(x-> x.== edelim, delim)
del = delim[cdelim]+1 : delim[cdelim+1]-1
#del = delim
plat = Vector{Vector{Int64}}(undef, 0)
#plat = Dict()
plat_ = Int64.(f[del, 2:3])
head_ = String.(f[del, 1])
for k = 1:length(mu_list)
mu = mu_list[k]
if mul in mu && mu[1] != mu[2] && !(mus in mu) #heavy-light
n = findfirst(x-> occursin("lh_vec", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif mul in mu && mu[1] == mu[2] #light-light
n = findfirst(x-> occursin("ll_vec", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif mul in mu && mus in mu#strange-light
n = findfirst(x-> occursin("ls_vec", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif mus in mu && mu[1] != mu[2] && !(mul in mu)#heavy-strange
n = findfirst(x-> occursin("sh_vec", x ), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif mus in mu && mu[1] == mu[2] && !(mul in mu)#strange-strange
n = findfirst(x-> occursin("ss_vec", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif !(mul in mu) && !(mus in mu) #heavy-heavy
n = findfirst(x-> occursin("hh_vec", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
end
end
return plat
end
function select_plateau_w(ensinf::EnsInfo, kappa_list)
ens = ensinf.id
deg = ensinf.deg
kappal, kappas = get_kappa(kappa_list, deg)
f = readdlm(path_plat_w)
head = String.(f[:, 1])
delim = findall(x-> occursin("#", x ), head)
edelim = findfirst(x-> occursin(ens, x), head)
cdelim = findfirst(x-> x.== edelim, delim)
del = delim[cdelim]+1 : delim[cdelim+1]-1
#del = delim
plat = Vector{Vector{Int64}}(undef, 0)
#plat = Dict()
plat_ = Int64.(f[del, 2:3])
head_ = String.(f[del, 1])
for k = 1:length(kappa_list)
kappa = kappa_list[k]
if kappal in kappa && kappa[1] == kappa[2] #light-light
n = findfirst(x-> occursin("ll", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif kappal in kappa && kappas in kappa#strange-light
n = findfirst(x-> occursin("ls", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
elseif kappas in kappa && kappa[1] == kappa[2] && !(kappal in kappa)#strange-strange
n = findfirst(x-> occursin("ss", x), head_)
push!(plat, [plat_[n, 1], plat_[n, 2]] )
end
end
return plat
end
function match_muc(muh, m_lh, m_sh, target)
M = (2/3 .* m_lh .+ 1/3 .* m_sh)
par, chi2exp = lin_fit(muh, M)
muh_target = x_lin_fit(par, target)
return muh_target
end
function aic_weight(aic)
w_aux = exp.(-1/2 * (aic / minimum(aic) )) #/minimum(aic) -> +0.5*mininum(aic)
w_aux1=filter(x->x<10000, w_aux)
index = findall(x->x<10000, w_aux)
norm = sum(w_aux1)
w = w_aux1 ./ norm
return w, index
end
function aic_weight(c::Cat)
aic = c.aic
w_aux = exp.(-1/2 * (aic / minimum(aic) ) ) #/minimum(aic) -> +0.5*mininum(aic)
w_aux1=filter(x->x<10000, w_aux)
index = findall(x->x<10000, w_aux)
norm = sum(w_aux1)
w = w_aux1 ./ norm
return w, index
end
function aic_syst(aic, fit_res)
w,index = aic_weight(aic)
aux1 = sum( w .* fit_res[index].^2)
aux2 = sum( w .* fit_res[index])^2
return sqrt(abs(value(aux1 - aux2)))
end
function aic_syst(c::Cat)
w,index = aic_weight(c.aic)
aux1 = sum( w .* c.fit_res[index].^2)
aux2 = sum( w .* c.fit_res[index])^2
return sqrt(abs(value(aux1 - aux2)))
end
function aic_model_ave(aic, fit_res, wpm::Union{Nothing, Dict}=nothing)
w, index = aic_weight(aic)
res = sum( w.* fit_res[index])
isnothing(wpm) ? uwerr(res) : uwerr(res, wpm)
return res
end
function aic_model_ave(c::Cat, wpm::Union{Nothing, Dict}=nothing)
w, index = aic_weight(c)
res = sum( w.* getfield(c, :fit_res)[index])
isnothing(wpm) ? uwerr(res) : uwerr(res, wpm)
return res
end
function category_av(cat_arr::Vector{Cat})
best_res = aic_model_ave.(cat_arr)
syst = aic_syst.(cat_arr)
tot_err = err.(best_res).^2 .+ syst.^2
w_aux = 1 ./ (tot_err)
w = w_aux ./ sum(w_aux)
println(findmax(w)[2])
aux1 = sum( w .* best_res.^2)
aux2 = sum( w .* best_res)^2
sys1 = sqrt(abs(value(aux1 - aux2)))
av = sum( w .* best_res)
return av, sys1
end
function category_av_tot(cat_arr::Vector{Cat})
fit_res = vcat(getfield.(cat_arr, :fit_res)...)
aic = vcat(getfield.(cat_arr, :aic)...)
best = aic_model_ave(aic, fit_res)
syst = aic_syst(aic, fit_res)
return best, syst
end
function syst_av(syst_err::Vector{Float64})
aux = sum(syst_err.^2)
return sqrt(aux)
end
function cat_obs(obs::Vector{Vector{uwreal}})
N = length(obs)
if !all(length.(obs) .== N)
error("Dimension mismatch")
end
r = getindex.(obs, 1)
for i = 2:N
r = vcat(r, getindex.(obs, i))
end
return r
end
function cat_obs_vcat(obs::Vector{Vector{uwreal}})
return vcat(obs...)
end
function fit_over_cat!(cat_tot::Vector{Cat}, models::Vector{Vector{Function}}, phi2_phys::uwreal, model_param::Vector{Int64} ; decay::Bool=false,
weight_dof::Int64=2, xycovar::Bool=false, correlated_fit::Bool=false, Mrat::uwreal=uwreal(1.0),
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, small_sample::Bool=false)
models = vcat(models...)
function continuum_dependece(params, phi2, phih; decay_phih_dep::Bool=decay)
if !decay_phih_dep
return Mrat * (params[1] + params[2]*phi2 + params[3]*phih)
else
return params[1] + params[2]*phi2 + params[3]/sqrt(phih)
end
end
for cat in cat_tot
println("\n", cat.info, "\n")
for (i, mod) in enumerate(models)
if !xycovar # correlation between x and y
if !correlated_fit
fit_params, chi2_vs_chi2exp = juobs.fit_routine(mod, value.(cat.x_tofit), cat.obs, model_param[i], wpm=wpm)
else
cov_inv = juobs.inv_covar_multi_id(cat.obs)
fit_params, chi2_vs_chi2exp = juobs.fit_routine(mod, value.(cat.x_tofit), cat.obs, model_param[i], inv_cov=cov_inv, wpm=wpm)
end
else
fit_params, chi2_vs_chi2exp = juobs.fit_routine(mod, cat.x_tofit, cat.obs, model_param[i], covar=true, wpm=wpm)
end
val_mod_i = continuum_dependece(fit_params, phi2_phys, cat.phih_ph)
push!(cat.dof, length(cat.obs)-model_param[i] )
if !small_sample
AIC = chi2_vs_chi2exp * cat.dof[i] + weight_dof * model_param[i]
else
AIC = chi2_vs_chi2exp * cat.dof[i] + weight_dof * model_param[i] + ( weight_dof*model_param[i]^2 + weight_dof*model_param[i]) / (length(cat.obs)- model_param[i] -1)
end
tmp = [val_mod_i, fit_params, chi2_vs_chi2exp, model_param[i], AIC]
symb = [:fit_res, :fit_param, :chi2_vs_exp, :n_param, :aic]
for (k,ss) in enumerate(symb)
push!(getfield(cat, ss), tmp[k] )
end
end
end
end
function Base.:*(x::uwreal, y::Array{Any})
N = size(y, 1)
return fill(x, N) .* y
end
function Base.:*(x::uwreal, y::Array{Float64})
N = size(y, 1)
return fill(x, N) .* y
end
function Base.:*(x::uwreal, y::Array{uwreal})
N = size(y, 1)
return fill(x, N) .* y
end
function Base.:+(x::uwreal, y::Vector{uwreal})
N = size(y, 1)
return fill(x, N) .+ y
end
function Base.:+(x::Float64, y::Vector{Float64})
N = size(y, 1)
return fill(x, N) .+ y
end
function Base.:+(x::Float64, y::Vector{uwreal})
N = size(y, 1)
return fill(x, N) .+ y
end
function Base.:-(x::Float64, y::Vector{Any})
N = size(y, 1)
return fill(x, N) .- y
end
function Base.:-(x::Float64, y::Vector{Float64})
N = size(y, 1)
return fill(x, N) .- y
end
function Base.:-(x::Float64, y::Vector{uwreal})
N = size(y, 1)
return fill(x, N) .- y
end
\ No newline at end of file
src/gevp analysis/types.jl 0 → 100644
View file @ 6c956856
mutable struct EnsInfo
id::String
L::Int64
beta::Float64
deg::Bool
mpi::Float64
dtr::Int64
trunc::Union{Int64, Vector{Int64}, Nothing}
function EnsInfo(ens_id::String, info::Vector{Float64})
id = ens_id
L = info[1]
beta = info[2]
deg = info[3]
mpi = info[4]
dtr = info[5]
return new(id, L, beta, deg, mpi, dtr, nothing)
end
function EnsInfo(ens_id::String, info::Vector{Float64}, trunc::Union{Nothing, Int64, Vector{Int64}})
id = ens_id
L = info[1]
beta = info[2]
deg = info[3]
mpi = info[4]
dtr = info[5]
return new(id, L, beta, deg, mpi, dtr, trunc)
end
end
mutable struct EnsObs
ensinfo::EnsInfo
mu_list:: Vector{Vector{Float64}}
is_pseudo::Bool
m_ll::Union{Nothing,uwreal}
m_ls::Union{Nothing,uwreal}
m_lh::Union{Nothing,Vector{uwreal}}
m_ss::Union{Nothing,uwreal}
m_sh::Union{Nothing,Vector{uwreal}}
m_hh::Union{Nothing,Vector{uwreal}}
f_ll::Union{Nothing,uwreal}
f_ls::Union{Nothing,uwreal}
f_lh::Union{Nothing,Vector{uwreal}}
f_ss::Union{Nothing,uwreal}
f_sh::Union{Nothing,Vector{uwreal}}
f_hh::Union{Nothing,Vector{uwreal}}
m_ll_vec::Union{Nothing,uwreal}
m_ls_vec::Union{Nothing,uwreal}
m_lh_vec::Union{Nothing,Vector{uwreal}}
m_ss_vec::Union{Nothing,uwreal}
m_sh_vec::Union{Nothing,Vector{uwreal}}
m_hh_vec::Union{Nothing,Vector{uwreal}}
f_ll_vec::Union{Nothing,uwreal}
f_ls_vec::Union{Nothing,uwreal}
f_lh_vec::Union{Nothing,Vector{uwreal}}
f_ss_vec::Union{Nothing,uwreal}
f_sh_vec::Union{Nothing,Vector{uwreal}}
f_hh_vec::Union{Nothing,Vector{uwreal}}
m_lh_match::Union{Nothing,uwreal}
m_sh_match::Union{Nothing,uwreal}
m_hh_match::Union{Nothing,uwreal}
f_lh_match::Union{Nothing,uwreal}
f_sh_match::Union{Nothing,uwreal}
f_hh_match::Union{Nothing,uwreal}
m_lh_vec_match::Union{Nothing,uwreal}
m_sh_vec_match::Union{Nothing,uwreal}
m_hh_vec_match::Union{Nothing,uwreal}
f_lh_vec_match::Union{Nothing,uwreal}
f_sh_vec_match::Union{Nothing,uwreal}
f_hh_vec_match::Union{Nothing,uwreal}
muh_target::Union{Nothing,uwreal}
a::Union{Nothing,uwreal}
t0::Union{Nothing,uwreal}
deltam::Union{Nothing, uwreal}
function EnsObs(ens::EnsInfo, mu::Vector{Vector{Float64}}, m_ps::Vector{uwreal}, f_ps::Vector{uwreal}, m_vec::Vector{uwreal}, f_vec::Vector{uwreal})
a = new()
a.ensinfo = ens
a.mu_list = mu
a.is_pseudo = true
a.m_ll = get_ll(a.mu_list, m_ps,a.ensinfo.deg)
a.ensinfo.deg ? a.m_ls = a.m_ll : a.m_ls = get_ls(a.mu_list, m_ps, a.ensinfo.deg)
a.ensinfo.deg ? a.m_ss = a.m_ll : a.m_ss = get_ss(a.mu_list, m_ps, a.ensinfo.deg)
a.m_lh = get_lh(a.mu_list, m_ps,a.ensinfo.deg)
a.ensinfo.deg ? a.m_sh = a.m_lh : a.m_sh = get_sh(a.mu_list, m_ps, a.ensinfo.deg)
a.m_hh = get_hh(a.mu_list, m_ps,a.ensinfo.deg)
a.f_ll = get_ll(a.mu_list, f_ps, a.ensinfo.deg)
a.ensinfo.deg ? a.f_ls = a.f_ll : a.f_ls = get_ls(a.mu_list, f_ps, a.ensinfo.deg)
a.ensinfo.deg ? a.f_ss = a.f_ll : a.f_ss = get_ss(a.mu_list, f_ps, a.ensinfo.deg)
a.f_lh = get_lh(a.mu_list, f_ps,a.ensinfo.deg)
a.ensinfo.deg ? a.f_sh = a.f_lh : a.f_sh = get_sh(a.mu_list, f_ps, a.ensinfo.deg)
a.f_hh = get_hh(a.mu_list, f_ps,a.ensinfo.deg)
a.m_ll_vec = get_ll(a.mu_list, m_vec,a.ensinfo.deg)
a.ensinfo.deg ? a.m_ls_vec = a.m_ll_vec : a.m_ls_vec = get_ls(a.mu_list, m_vec, a.ensinfo.deg)
a.ensinfo.deg ? a.m_ss_vec = a.m_ll_vec : a.m_ss_vec = get_ss(a.mu_list, m_vec, a.ensinfo.deg)
a.m_lh_vec = get_lh(a.mu_list, m_vec, a.ensinfo.deg)
a.ensinfo.deg ? a.m_sh_vec = a.m_lh_vec : a.m_sh_vec = get_sh(a.mu_list, m_vec, a.ensinfo.deg)
a.m_hh_vec = get_hh(a.mu_list, m_vec, a.ensinfo.deg)
a.f_ll_vec = get_ll(a.mu_list, f_vec, a.ensinfo.deg)
a.ensinfo.deg ? a.f_ls_vec = a.f_ll_vec : a.f_ls_vec = get_ls(a.mu_list, f_vec, a.ensinfo.deg)
a.ensinfo.deg ? a.f_ss_vec = a.f_ll_vec : a.f_ss_vec = get_ss(a.mu_list, f_vec, a.ensinfo.deg)
a.f_lh_vec = get_lh(a.mu_list, f_vec,a.ensinfo.deg)
a.ensinfo.deg ? a.f_sh_vec = a.f_lh_vec : a.f_sh_vec = get_sh(a.mu_list, f_vec, a.ensinfo.deg)
a.f_hh_vec = get_hh(a.mu_list, f_vec, a.ensinfo.deg)
a.m_lh_match = nothing
a.m_sh_match = nothing
a.m_hh_match = nothing
a.f_lh_match = nothing
a.f_sh_match = nothing
a.f_hh_match = nothing
a.m_lh_vec_match = nothing
a.m_sh_vec_match = nothing
a.m_hh_vec_match = nothing
a.f_lh_vec_match = nothing
a.f_sh_vec_match = nothing
a.f_hh_vec_match = nothing
a.muh_target = nothing
a.a = nothing
a.t0 = nothing
a.deltam = nothing
return a
end
function EnsObs(ens::EnsInfo, mu::Vector{Vector{Float64}}, m_ps::Vector{uwreal}, m_vec::Vector{uwreal})
a = new()
a.ensinfo = ens
a.mu_list = mu
a.is_pseudo = true
a.m_ll = get_ll(a.mu_list, m_ps,a.ensinfo.deg)
a.ensinfo.deg ? a.m_ls = a.m_ll : a.m_ls = get_ls(a.mu_list, m_ps, a.ensinfo.deg)
a.ensinfo.deg ? a.m_ss = a.m_ll : a.m_ss = get_ss(a.mu_list, m_ps, a.ensinfo.deg)
a.m_lh = get_lh(a.mu_list, m_ps,a.ensinfo.deg)
a.ensinfo.deg ? a.m_sh = a.m_lh : a.m_sh = get_sh(a.mu_list, m_ps, a.ensinfo.deg)
a.m_hh = get_hh(a.mu_list, m_ps,a.ensinfo.deg)
a.m_ll_vec = get_ll(a.mu_list, m_vec,a.ensinfo.deg)
a.ensinfo.deg ? a.m_ls_vec = a.m_ll_vec : a.m_ls_vec = get_ls(a.mu_list, m_vec, a.ensinfo.deg)
a.ensinfo.deg ? a.m_ss_vec = a.m_ll_vec : a.m_ss_vec = get_ss(a.mu_list, m_vec, a.ensinfo.deg)
a.m_lh_vec = get_lh(a.mu_list, m_vec, a.ensinfo.deg)
a.ensinfo.deg ? a.m_sh_vec = a.m_lh_vec : a.m_sh_vec = get_sh(a.mu_list, m_vec, a.ensinfo.deg)
a.m_hh_vec = get_hh(a.mu_list, m_vec, a.ensinfo.deg)
a.m_lh_match = nothing
a.m_sh_match = nothing
a.m_hh_match = nothing
a.m_lh_vec_match = nothing
a.m_sh_vec_match = nothing
a.m_hh_vec_match = nothing
a.muh_target = nothing
a.a = nothing
a.t0 = nothing
return a
end
end
mutable struct MatInfo
ensinfo::EnsInfo
mu::Vector{Float64}
mat_list::Vector{Matrix{uwreal}}
y0::Int64
function MatInfo(_ensinfo::EnsInfo, _mat_list::Array{Array{T,2} where T,1}, _mu::Vector{Float64}, _y0::Int64)
a = new()
a.ensinfo = _ensinfo
a.mu = _mu
a.mat_list = _mat_list
a.y0 = _y0
return a
end
end
mutable struct Cat
obs::Vector{uwreal}
phih::Vector{uwreal}
phih_ph::uwreal
#models::Vector{Vector{Function}}
x_tofit::Array{uwreal}
info::String
fit_param::Vector{Vector{uwreal}} # store all fit parameters
fit_res::Vector{uwreal} # store all fit results
chi2_vs_exp::Vector{Float64} # store all chi2 corrected
n_param::Vector{Int64} # store the number of parameters. Maybe this is redundant
aic::Vector{Float64} # store the AIC value
dof::Vector{Int64}
function Cat(_obs::Vector{uwreal}, _xtofit::Array{uwreal}, _phi_ph::uwreal, _info::String)
a = new()
a.obs = _obs
a.info = _info
a.x_tofit = _xtofit
a.phih_ph = _phi_ph
a.fit_param = Vector{Vector{uwreal}}(undef, 0)
a.fit_res = Vector{uwreal}(undef, 0)
a.chi2_vs_exp = Vector{Float64}(undef, 0)
a.n_param = Vector{Int64}(undef, 0)
a.aic = Vector{Float64}(undef, 0)
a.dof = Vector{Int64}(undef, 0)
return a
end
end
\ No newline at end of file
src/gevp_mass_shift.jl deleted 100644 → 0
View file @ 838f8293
using Revise, ADerrors, BDIO, juobs, PyPlot, LaTeXStrings, DelimitedFiles, Combinatorics
using LsqFit, juobs.LeastSquaresOptim
include("gevp_types.jl")
include("gevp_reader.jl")
include("gevp_tools.jl")
include("const.jl")
include("gevp_func_comb.jl")
#============= INPUT PAR ===========#
# desktop parameters
#include("/home/alessandro/juobs/constants/juobs_const.jl")
#path_data ="/home/alessandro/Desktop/data/heavy_3pt/3pt"
#path_rw = "/media/alessandro/4277fef2-edc5-4e0d-89cb-f5d1d44fbc8c/data/rwf"
#path_plot = "/home/alessandro/google-drive/phd/secondment/3pf test/analysis/plots"
# macbook parameters
include("/Users/ale/juobs/constants/juobs_const.jl")
path_plot = "/Users/ale/Il mio Drive/phd/analysis/charm_gevp_2022/plots"
path_results = "/Users/ale/Il mio Drive/phd/analysis/charm_gevp_2022/bdio"
# latex-style labels
PyPlot.rc("font", family="sans serif", size=13)
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["text.usetex"] = false
rcParams["mathtext.fontset"] = "cm"
#rcParams["font.size"] =11
rcParams["axes.labelsize"] =16
const ens_list = ["H101", "H102r002", "H400", "N202", "N203", "N200", "J303", "N300"]
const sector = Dict("pseudo"=>true, "vector"=>true)
# loading ensemble information
ensinfo = Vector{EnsInfo}(undef, length(ens_list))
for i in 1:length(ens_list)
ens = ens_list[i]
try
ensinfo[i] = EnsInfo(ens, ens_db[ens] )
catch
error("The ensemble id ", ens, " was not found in the const.jl ens_db database. Please check the ensemble id or update the database")
end
end
path_dm = "/Users/ale/Desktop/data"
dm = read_BDIO(path_dm, "delta")
t0_shifted = read_BDIO(path_dm, "t0")
#test = read_BDIO(joinpath(path_results,"ps_masses"), "H101" )
#uwerr.(test)
\ No newline at end of file
src/gevp_continuum.jl → src/older analysis/gevp_continuum.jl
View file @ 6c956856
......@@ -9,13 +9,53 @@
#
##################################
using Revise, ADerrors, BDIO, juobs, PyPlot, LaTeXStrings, DelimitedFiles, Combinatorics
using LsqFit, juobs.LeastSquaresOptim
using LsqFit
include("gevp_types.jl")
include("gevp_reader.jl")
include("gevp_tools.jl")
include("const.jl")
include("gevp_func_comb.jl")
include("gevp_plotter.jl")
# redefine + * operations
import Base: +,*,/
#uwreal +
function +(uw::uwreal, vec::Vector{uwreal})
n = size(vec,1)
return fill(uw, n) .+ vec
end
function +(vec::Vector{uwreal}, uw::uwreal ) uw+vec end
#uwreal *
function *(uw::uwreal, vec::Vector{uwreal})
n = size(vec,1)
return fill(uw, n) .* vec
end
function *(vec::Vector{uwreal}, uw::uwreal ) uw*vec end
#any +
function +(uw::uwreal, vec::Vector{Any})
n = size(vec,1)
return fill(uw, n) .+ vec
end
function +(vec::Vector{Any}, uw::uwreal ) uw+vec end
#any *
function *(uw::uwreal, vec::Vector{Any})
n = size(vec,1)
return fill(uw, n) .* vec
end
function *(vec::Vector{Any}, uw::uwreal ) uw*vec end
#float +
function +(uw::uwreal, vec::Vector{Float64})
n = size(vec,1)
return fill(uw, n) .+ vec
end
function +(vec::Vector{Float64}, uw::uwreal ) uw+vec end
#float *
function *(uw::uwreal, vec::Vector{Float64})
n = size(vec,1)
return fill(uw, n) .* vec
end
function *(vec::Vector{Float64}, uw::uwreal ) uw*vec end
#============= INPUT PAR ===========#
......@@ -26,11 +66,17 @@ include("gevp_func_comb.jl")
#path_rw = "/media/alessandro/4277fef2-edc5-4e0d-89cb-f5d1d44fbc8c/data/rwf"
#path_plot = "/home/alessandro/google-drive/phd/secondment/3pf test/analysis/plots"
# macbook parameters
include("/Users/ale/juobs/constants/juobs_const.jl")
path_plot = "/Users/ale/Il mio Drive/phd/analysis/charm_gevp_2022/plots"
path_results = "/Users/ale/Il mio Drive/phd/analysis/charm_gevp_2022/bdio"
path_t0_and_dm = "/Users/ale/Desktop/data"
# macbook 15 parameters
#include("/Users/ale/juobs/constants/juobs_const.jl")
#path_plot = "/Users/ale/Il mio Drive/phd/analysis/charm_gevp_2022/plots"
#path_results = "/Users/ale/Il mio Drive/phd/analysis/charm_gevp_2022/bdio"
#path_t0_and_dm = "/Users/ale/Desktop/data"
include("/Users/alessandroconigli/Lattice/juobs/constants/juobs_const.jl")
path_plot = "/Users/alessandroconigli/MyDrive/phd/analysis/charm_gevp_2022/plots"
path_results = "/Users/alessandroconigli/MyDrive/phd/analysis/charm_gevp_2022/bdio"
path_t0_and_dm = "/Users/alessandroconigli/Lattice/data"
# latex-style labels
PyPlot.rc("font", family="sans serif", size=12)
......@@ -38,7 +84,7 @@ rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["text.usetex"] = false
rcParams["mathtext.fontset"] = "cm"
#rcParams["font.size"] =11
rcParams["axes.labelsize"] =14
rcParams["axes.labelsize"] =22
......@@ -101,7 +147,7 @@ end
#######################
# read t0_shifted from BDIO
dm = read_BDIO(path_t0_and_dm, "delta")
#dm = read_BDIO(path_t0_and_dm, "delta")
t0_shifted = read_BDIO(path_t0_and_dm, "t0")
uwerr.(t0_shifted)
label = vcat(ensembles.(t0_shifted)...)
......@@ -224,19 +270,37 @@ xdata_spav = [ 1 ./ sqrt8t0.^2 phi2 phih_spav]
wpmm = Dict{String, Vector{Float64}}()
wpmm["N202"] = [-1.0, 2.0,-1.0, -1.0]
#wpmm["N200"] = [-1.0, 2.0,-1.0, -1.0]
wpmm["J303"] = [-1.0, 2.0,-1.0, -1.0]
#wpmm[300] = [-1.0, 1.5,-1.0, -1.0]
cat_muc_eta = Cat(vcat(muc_t0...), xdata_eta, phih_eta_ph, "muc_etac_matching")
cat_muc_flav = Cat(vcat(muc_t0...), xdata_flav, phih_flav_ph, "muc_flav_matching")
cat_muc_spav = Cat(vcat(muc_t0...), xdata_spav, phih_spav_ph, "muc_spav_matching")
cat_muc_eta = Cat(vcat(muc_t0...), xdata_eta, phih_eta_ph, "muc_etac")
cat_muc_flav = Cat(vcat(muc_t0...), xdata_flav, phih_flav_ph, "muc_flav")
#cat_muc_spav = Cat(vcat(muc_t0...), xdata_spav, phih_spav_ph, "muc_spav_matching")
cat_muc_tot = [cat_muc_eta, cat_muc_flav, cat_muc_spav]
cat_muc_tot = [cat_muc_eta, cat_muc_flav]
# fit all models for categories
fit_over_cat!(cat_muc_tot, model_charm, phi2_ph, n_param_j, weight_dof=2 )
muc_av_eta = aic_model_ave(cat_muc_eta)
muc_av_flav = aic_model_ave(cat_muc_flav)
# working plotter
idx_toplot = 1
plot_cat(cat_muc_tot[idx_toplot], L"$M_c^{\mathrm{RGI}}(N_f=3)\mathrm{MeV}$" , save=false)
plot_chi_charm(cat_muc_tot[idx_toplot], vcat(model_charm_aux...), L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, n=1)
plot_cl_charm(cat_muc_tot[idx_toplot], vcat(model_charm_aux...), L"$ciao$",phi2_ph, ensinfo, mrat=Mrat, n=1 )
plot_all_cl_charm(cat_muc_tot[1:2], vcat(model_charm_aux...),L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, mrat=Mrat )# , p=path_plot)
plot_chi2_vs_exp(cat_muc_flav, label_cutoff)
##
model_j = vcat(f...)
for cat in cat_muc_tot[1:2]
for i in 1:length(model_charm)
#try
par, chi2_over_chiexp = fit_routine(model_charm[i], value.(cat.x_tofit), cat.obs, n_param_charm[i], wpm=wpmm)
par, chi2_over_chiexp = fit_routine(model_j[i], value.(cat.x_tofit), cat.obs, n_param_charm[i], wpm=wpmm)
aux = Mrat * (par[1] + par[2]*phi2_ph + par[3]*cat.phih_ph)
uwerr(aux)
println("charm mass = ", aux)
......@@ -245,7 +309,7 @@ for cat in cat_muc_tot[1:2]
push!(cat.chi2_corr, chi2_over_chiexp)
push!(cat.n_param, n_param_charm[i])
push!(cat.dof, length(cat.obs)-n_param_charm[i] )
push!(cat.aic, cat.chi2_corr[i]*cat.dof[i] +2*n_param_charm[i])
push!(cat.aic, cat.chi2_corr[i]*cat.dof[i] +5*n_param_charm[i])
println("aic = ", cat.aic[i])
println("\n")
#catch
......@@ -256,19 +320,49 @@ end
muc_av_eta = aic_model_ave(cat_muc_eta)
muc_av_flav = aic_model_ave(cat_muc_flav)
mc_fin, mc_syst = category_av([cat_muc_eta, cat_muc_flav])
muc_av_spav = aic_model_ave(cat_muc_spav)
plot_cat(cat_muc_flav, L"$M_c^{\mathrm{RGI}}(N_f=3)\mathrm{MeV}$" , save=true)
plot_cont_lim_phi2(cat_muc_eta, model_charm, L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", dec=false, save=false)
plot_cont_lim(cat_muc_eta, model_charm, L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", save=false)
plot_cat(cat_muc_eta, L"$M_c^{\mathrm{RGI}}(N_f=3)\mathrm{MeV}$" , save=false)
plot_cl_charm(cat_muc_eta, model_charm, L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo)
plot_all_cl_charm(cat_muc_tot[1:2], model_charm,L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, n=[0,1,12] )# , p=path_plot)
plot_chi_charm(cat_muc_flav, model_charm, L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo)
plot_all_cl_charm(cat_muc_tot[1:2], model_charm,L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, n=[0,0,12] )# , p=path_plot)
plot_hist(cat_muc_tot[1:2] )
## final results
mc_fin, mc_syst = category_av([cat_muc_eta, cat_muc_flav])
mc_fin_ph = mc_fin * hc /t0_ph[1]
mc_syst_ph = mc_syst * hc /t0_ph[1]
uwerr(mc_fin_ph)
uwerr(mc_syst_ph)
println("Final result for RGI charm quark mass in physical units: ", mc_fin_ph)
println("Systematic error coming from model average: ", mc_syst_ph)
######
## test muc with javier functions
for cat in cat_muc_tot[1:1]
println("\n", cat.info, "\n")
for j in 1:length(f) # looping over number of parameters
for k in 1:length(f[j]) # looping over combinations with fixed parameters
println(" j is :", j, " k is equal: ", k)
#try
par, chi_corr = fit_routine(f[j][k], value.(cat.x_tofit), cat.obs, n_param2[j], wpm=wpmm)
aux = par[1] + par[2]*phi2_ph + par[3]*(cat.phih_ph)
println(aux)
uwerr(aux)
push!(cat.fit_res, aux)
push!(cat.fit_param, par)
push!(cat.chi2_corr, chi_corr)
push!(cat.n_param, n_param2[j])
push!(cat.aic, chi_corr +2*n_param2[j])
#catch
# println("There was a problem... Probably the analysis failed for some error id.")
# println("The result has been excluded from the category")
#end
end
end
end
#########################
......@@ -282,7 +376,7 @@ cat_mDs_eta = Cat(vcat(mps_sh_t0...), xdata_eta, phih_eta_ph, "mDs_etac_matchin
cat_mDs_flav = Cat(vcat(mps_sh_t0...), xdata_flav, phih_flav_ph, "mDs_flav_matching")
cat_D_Ds_masses =[cat_mD_eta, cat_mD_flav, cat_mDs_eta, cat_mDs_flav]
for cat in cat_D_Ds_masses[1:end]
for cat in cat_D_Ds_masses[1:1]
println(cat.info)
for i in 1:length(model_charm)
#try
......@@ -294,9 +388,9 @@ for cat in cat_D_Ds_masses[1:end]
push!(cat.fit_param, par)
push!(cat.chi2_corr, chi2_over_chi2exp)
push!(cat.n_param, n_param_charm[i])
push!(cat.dof, length(cat.obs)-n_param_charm[i] )
push!(cat.aic, chi2_over_chi2exp*cat.dof[i] +2*n_param_charm[i])
push!(cat.dof, length(cat.obs)-n_param_charm[i] )
push!(cat.aic, cat.chi2_corr[i]*cat.dof[i] +2*n_param_charm[i])
println("chi2 / chiexp = ", chi2_over_chi2exp)
println("\n")
#catch
......@@ -316,11 +410,11 @@ plot_cont_lim(cat_muc_eta, model_charm, L"$\sqrt{8t_0}M_{D_s}$", save=false)
plot_cat(cat_muc_flav, L"$\sqrt{8t_0}M_{D}$" )
################
## fit for Fd
## fit for fD(s)
################
xdata_eta_dec[:] = xdata_eta
xdata_eta_dec = deepcopy(xdata_eta)
xdata_eta_dec[:,3] = 1.0 ./sqrt.(xdata_eta[:,3])
xdata_flav_dec = xdata_flav
xdata_flav_dec = deepcopy(xdata_flav)
xdata_flav_dec[:,3] = 1.0 ./sqrt.(xdata_flav[:,3])
cat_fD_eta = Cat(vcat(fps_lh_t0...), xdata_eta, phih_eta_ph, "fD_etac_matching")
......@@ -330,53 +424,108 @@ cat_fDs_eta = Cat(vcat(fps_sh_t0...), xdata_eta, phih_eta_ph, "fDs_etac_matchin
cat_fDs_flav = Cat(vcat(fps_sh_t0...), xdata_flav, phih_flav_ph, "fDs_flav_matching")
cat_fD_fDs = [cat_fD_eta, cat_fD_flav, cat_fDs_eta, cat_fDs_flav]
for cat in cat_fD_fDs[1:4]
for cat in cat_fD_fDs[1:1]
println(cat.info)
for i in 1:10#length(model_decays)
try
par, chi2_over_chi2exp = fit_routine(model_decays[i], cat.x_tofit, cat.obs, n_param_decays[i], covar=false, wpm=wpmm)
aux = par[1] + par[2]*phi2_ph + par[3]/sqrt((cat.phih_ph))
for i in 1:length(model_decay)
#try
par, chi2_over_chi2exp = fit_routine(model_decay[i], value.(cat.x_tofit), cat.obs, n_param_decay[i], wpm=wpmm)
aux = par[1] + par[2]*phi2_ph + par[3] * cat.phih_ph # /sqrt((cat.phih_ph))
uwerr(aux)
println("fD(s) decays = ", aux)
push!(cat.fit_res, aux)
push!(cat.fit_param, par)
push!(cat.chi2_corr, chi2_over_chi2exp)
push!(cat.n_param, n_param_decays[i])
cat.dof = length(cat.obs)-n_param_decays[i]
push!(cat.aic, chi2_over_chi2exp*cat.dof +4*n_param_decays[i])
push!(cat.n_param, n_param_decay[i])
push!(cat.dof, length(cat.obs)-n_param_decay[i] )
push!(cat.aic, cat.chi2_corr[i]*cat.dof[i] +2*n_param_decay[i])
println("chi2 / chiexp = ", chi2_over_chi2exp)
println("aic = ", cat.aic[i])
println("\n")
catch
println("something went wrong")
#catch
# println("something went wrong")
#end
end
end
fD_av_eta = aic_model_ave(cat_fD_eta)
fD_av_flav = aic_model_ave(cat_fD_flav)
fDs_av_eta = aic_model_ave(cat_fDs_eta)
fDs_av_flav = aic_model_ave(cat_fDs_flav)
plot_cont_lim_phi2(cat_fD_eta, model_decays, L"$\sqrt{8t_0}f_{D_s}$", dec=true, save=false)
plot_cat(cat_fD_eta, L"ciao")
plot_cl_charm(cat_fD_eta, model_decay, L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, n=1)
plot_chi_charm(cat_fD_eta, model_decay, L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, n=2)
## testing fit charm with Javier function definition
cat_fD_eta = Cat(vcat(fps_lh_t0...), xdata_eta, phih_eta_ph, "fD_etac_matching")
cat_fD_flav = Cat(vcat(fps_lh_t0...), xdata_flav, phih_flav_ph, "fD_flav_matching")
cat_fDs_eta = Cat(vcat(fps_sh_t0...), xdata_eta, phih_eta_ph, "fDs_etac_matching")
cat_fDs_flav = Cat(vcat(fps_sh_t0...), xdata_flav, phih_flav_ph, "fDs_flav_matching")
cat_fD_fDs = [cat_fD_eta, cat_fD_flav, cat_fDs_eta, cat_fDs_flav]
for cat in cat_fD_fDs[1:1]
println("\n", cat.info, "\n")
for j in 1:length(f) # looping over number of parameters
for k in 1:length(f[j]) # looping over combinations with fixed parameters
println(" j is :", j, " k is equal: ", k)
#try
par, chi_corr = fit_routine(f[j][k], value.(cat.x_tofit), cat.obs, n_param2[j], wpm=wpmm)
aux = par[1] + par[2]*phi2_ph + par[3] / sqrt(cat.phih_ph)
println(aux)
uwerr(aux)
push!(cat.fit_res, aux)
push!(cat.fit_param, par)
push!(cat.chi2_corr, chi_corr)
push!(cat.n_param, n_param2[j])
push!(cat.aic, chi_corr +2*n_param2[j])
#catch
# println("There was a problem... Probably the analysis failed for some error id.")
# println("The result has been excluded from the category")
#end
end
end
end
fD_av_eta = aic_model_ave(cat_fD_eta)
fD_av_flav = aic_model_ave(cat_fD_flav)
fDs_av_eta = aic_model_ave(cat_fDs_eta)
fDs_av_flav = aic_model_ave(cat_fDs_flav)
plot_cont_lim_phi2(cat_fD_eta, model_decays, L"$\sqrt{8t_0}f_{D_s}$", dec=true, save=false)
plot_cat(cat_fD_flav, L"ciao")
plot_cat(cat_fD_eta, L"ciao")
plot_cl_charm(cat_fD_eta, model_decay, L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, n=1)
plot_chi_charm(cat_fD_eta, model_decay, L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, n=2)
## test fitting without multiplying the data for N_mh
xdat_test_fl = xdata_flav[:,:]
xdat_test_fl[:,1] = xdata_flav[1:3:end,1]
## testing fit charm
par_eta, chi = fit_routine(f[1][1], xdata_eta, vcat(muc_t0...), 5 , covar=true)
par_flav, chi = fit_routine(f[1][1], xdata_flav, vcat(muc_t0...), 5 , covar=true)
## test fit fD(s) with javier function
val_ph_eta = par_eta[1] + par_eta[2]*phi2_ph + par_eta[3]*phih_eta_ph
val_ph_flav = par_flav[1] + par_flav[2]*phi2_ph + par_flav[3]*phih_flav_ph
par_eta, chi = fit_routine(f[2][1], xdata_eta, vcat(fps_lh_t0...), 6 , covar=true)
par_flav, chi = fit_routine(f[1][1], xdata_flav, vcat(fps_lh_t0...), 5 , covar=true, wpm=wpmm)
par_eta2, chi = fit_routine(f[2][1], xdata_eta, vcat(fps_sh_t0...), 6 , covar=true)
par_flav, chi = fit_routine(f[1][1], xdata_flav, vcat(fps_lh_t0...), 5 , covar=true, wpm=wpmm)
test_fd_eta = par_eta[1] + par_eta[2]*phi2_ph + par_eta[3]*phih_eta_ph
test_fd_flav = par_flav[1] + par_flav[2]*phi2_ph + par_flav[3]*phih_flav_ph
uwerr(val_ph_eta)
uwerr(val_ph_flav)
##
muh_target_arr = []
fig = figure()
......
src/gevp_func_comb.jl → src/older analysis/gevp_func_comb.jl
View file @ 6c956856
# FUNCTIONS
#=
##FUNCTIONS
basemodel(x,p) = p[1] .+ p[2] .* x[:,2] .+ p[3] .* x[:,3] .+ p[4] .* x[:,1]
a2phi2(x) = x[:,1] .* x[:,2] # phi2*a^2/8t0
a2phih2(x) = x[:,1] .* x[:,3].^2 # phih^2*a^2/8t0
......@@ -8,14 +7,14 @@ a4phih4_2(x) = x[:,1].^(2) .* x[:,3].^2 #(a^2/8t0)^2*phih^2
a4cutoff(x) = x[:,1].^(2) #(a^2/8t0)^2
func_list = [a2phi2, a2phih2, a4phih4, a4phih4_2, a4cutoff]
label_list = ["a2l", "a2h2", "a4h4", "a4h2", "a4"]
# COMBINATIONS (LINEAR)
##COMBINATIONS (LINEAR)
func_map = [Bool.([i,j,k,m,n]) for i=0:1 for j=0:1 for k=0:1 for m=0:1 for n=0:1]
npar = length(func_list) # number of extra parameters
n_param = Vector(4:npar+4)
n_param_charm = Vector{Int64}(undef,0)
push!(n_param_charm, 4)
n_param_j = Vector(4:4)
label_cutoff = Vector{Vector{String}}(undef, 0)
push!(label_cutoff, ["a2"])
#f_lin definition
f_lin = Vector{Vector{Function}}(undef, npar+1) #f[i][j]-> i: number of extra parameters, j: combinations
f_lin[1] = Vector{Function}(undef, 1)
......@@ -29,13 +28,13 @@ f_aux[1][1] = (x,p) -> p[1] + p[2] * x[:,2] + p[3] * x[:,3] + p[4] * x[:,1]
#loop over combinations
for n = 2:npar+1
aux = filter(x->sum(x) == n-1, func_map)
println(n_param[n])
f_lin[n] = Vector{Function}(undef, length(aux))
f_aux[n] = Vector{Function}(undef, length(aux))
k = 1
for a in aux
push!(n_param_charm, n_param[n])
#vectorized function
push!(n_param_j, 4+n-1)
push!(label_cutoff, label_list[a])
f_lin[n][k] = (x,p) -> basemodel(x,p) .+ sum([p[i+4] for i=1:(n-1)] .* (fill(x,n-1) .|> func_list[a]))
f_aux[n][k] = (x,p) -> f_aux[1][1](x, p) + sum([p[i+4] for i=1:(n-1)] .* (fill(x,n-1) .|> func_list[a]))
k = k + 1
......@@ -43,7 +42,7 @@ for n = 2:npar+1
end
# COMBINATIONS (NONLINEAR)
##COMBINATIONS (NONLINEAR)
f_non_lin = Vector{Vector{Function}}(undef, npar) #f[i][j]-> i: number of extra parameters, j: combinations
f_non_lin_aux = Vector{Vector{Function}}(undef, npar)
n_non_lin_param = Vector(5:npar+4) #param for each non-linear functions
......@@ -54,26 +53,26 @@ for n = 1:npar
k = 1
for a in aux
#vectorized function
push!(n_param_charm, n_param[n]+1)
push!(n_param_j, 4+n)
f_non_lin[n][k] = (x,p) -> basemodel(x,p) .* (1 .+ sum([p[i+4] for i=1:(n)] .* (fill(x,n) .|> func_list[a])))
f_non_lin_aux[n][k] = (x,p) -> p[1] + p[2] * x[:,2] + p[3] * x[:,3] + p[4] * x[:,1] +
(p[1] + p[2] * x[:,2] + p[3] * x[:,3] + p[4] * x[:,1]) .* sum([p[i+4] for i=1:(n)] .* (fill(x,n) .|> func_list[a]))
f_non_lin_aux[n][k] = (x,p) -> p[1] + p[2] * x[:,2] + p[3] * (x[:,3]) + p[4] * x[:,1] +
(p[1] + p[2] * x[:,2] + p[3] *(x[:,3]) + p[4] * x[:,1]) .* sum([p[i+4] for i=1:(n)] .* (fill(x,n) .|> func_list[a]))
k = k + 1
end
end
# APPEND LIN + NONLIN
##APPEND LIN + NONLIN
f = f_lin
append!(f, f_non_lin)
append!(f_aux, f_non_lin_aux)
model_charm = vcat(f...)
model_charm_aux = vcat(f_aux...)
model_charm = append!(f, f_non_lin)
model_charm_aux = append!(f_aux, f_non_lin_aux)
n_param2 = vcat(n_param, n_non_lin_param)
=#
# #=
#=
import Base.+
import Base.*
......@@ -82,10 +81,10 @@ import Base.*
*(f::Function, g::Function) = (x...) -> f(x...) * g(x...)
*(f::Function, g::Function) = (x...) -> f(x...) .* g(x...)
# testing model with BDIO parameters to plot automatically the shaded band
# charm quark fits
function basemodel(x,p)
n=size(x,1)
aux = [p[1] for i in 1:n] .+ [p[2] for i in 1:n].*x[:,2] .+ [p[3] for i in 1:n].*x[:,3] .+ [p[4] for i in 1:n].*x[:,1]
aux = [p[1] for i in 1:n] .+ [p[2] for i in 1:n].*x[:,2] .+ [p[3] for i in 1:n] .* x[:,3] .+ [p[4] for i in 1:n].*x[:,1]
return aux
end
#a^2 cutoff effects
......@@ -117,6 +116,36 @@ for i in 1:length(func_comb_charm)
push!(model_charm, basemodel * (fff + sum(tmp_farr)))
end
# decay constansts
function basemodel_dec(x,p)
n=size(x,1)
aux = [p[1] for i in 1:n] .+ [p[2] for i in 1:n].*x[:,2] .+ [p[3] for i in 1:n] .* x[:,3] .+ [p[4] for i in 1:n].*x[:,1]
return aux
end
func_comb_decay = collect(combinations([ a2phi2, a2phih, a4, a4phih2, a4phih4]))
model_decay = Vector{Function}(undef, 0)
push!(model_decay, basemodel_dec)
n_param_decay = Vector{Int64}(undef,0)
push!(n_param_decay, 4)
for i in 1:length(func_comb_decay)
tmp_farr = Vector{Function}(undef,0)
for k in 1:length(func_comb_decay[i])
tmp_f = (x,p) -> [p[n_param_decay[1]+k] for n in 1:size(x,1)] .* func_comb_decay[i][k](x)
push!(tmp_farr, tmp_f)
end
push!(model_decay, basemodel_dec + sum(tmp_farr))
push!(n_param_decay, 4+length(func_comb_decay[i]))
#push!(n_param_decay, 4+length(func_comb_decay[i]))
#fff(x,p) = 1.
#aux_non_lin = (x,p) -> basemodel_dec(x,p) .+ [p[1] for l in 1:length(x[:,1])] .+ [p[2] for l in 1:length(x[:,1])].*x[:,2] .+ [p[3] for l in 1:length(x[:,1])].*x[:,3] .+ [p[4] for l in 1:length(x[:,1])].*x[:,1] #.* [ ([p[n_param_charm[1]+k] for n in 1:size(x,1)] .* func_comb_charm[i][k](x)) for k in 1:length(func_comb_charm[i])]
#push!(model_decay, basemodel_dec * (fff + sum(tmp_farr)))
end
=#
##
# =#
#=
......
src/gevp_obs_cf2bdio.jl → src/older analysis/gevp_obs_tm.jl
View file @ 6c956856
#################################
# DEPRECATED
# This file was created by Alessandro Conigli - 2022
# The pourpose of this code is to compute twisted mass
# observables from raw correlators stored in mesons.dat files.
......@@ -103,11 +104,11 @@ if sector["decays"]
errorbar(xx, value.(aux_data_ps), err.(aux_data_ps), fmt="s", mfc="none", label="pseudoscalar")
errorbar(xx, value.(aux_data_vec), err.(aux_data_vec), fmt="s", mfc="none", label="vector")
fill_between(plat_ps[1]:plat_ps[2], v-e, v+e, alpha=0.6 )
fill_between(plat_vec[1]:plat_vec[2], v_vec-e_vec, v_vec+e_vec, alpha=0.6 )
fill_between(plat_ps[1]+62:plat_ps[2+62], v-e, v+e, alpha=0.6 )
fill_between(plat_vec[1]+62:plat_vec[2]+62, v_vec-e_vec, v_vec+e_vec, alpha=0.6 )
xlabel(L"$x_0/a$", size=18)
ylabel(L"$af$", size=18)
title("$(ensembles(dec0_ps[i])[1]) " * L"\mu_1="*"$(mu_list[i][1]) " *L" \mu_2="*"$(mu_list[i][2]) " )
title("$(ensembles(dec0_ps[i])[1]) " * L"\mu_1="*"$(mu_list[i][1]) " *L" \mu_2="*"$(mu_list[i][2]) " , fontsize=18)
ylim(v*0.95, v_vec*1.1)
legend()
display(gcf())
......
src/gevp_plotter.jl → src/older analysis/gevp_plotter.jl
View file @ 6c956856
function plot_cont_lim_phi2(cat::Cat, ff::Vector{Function}, ylab::LaTeXString; dec::Bool=false,save::Bool=false)
aic, idx = findmin(getfield(cat, :chi2_corr)) # best fit
#idx=1
par_n = getfield(cat, :n_param)[idx] #here 2 shoudl be idx best n of parameters
model = ff[idx] #best model
fit_par = cat.fit_param[idx]
best_res = value(Mrat) * (fit_par[1] + fit_par[2]*phi2_ph + fit_par[3]*value(cat.phih_ph))
uwerr(best_res)
println("Results corresponding to the best fit: ", best_res)
println("AIC corresponding to the best fit ", aic)
println("chi2/chiexp corresponding to the best fit ", cat.chi2_corr[idx] )
x = getfield(cat,:x_tofit)
obs = cat.obs
obs_match = Vector{uwreal}(undef,0)
for i in 1:length(ensinfo)
xx = 1 / (8*t0(ensinfo[i].beta))
#phih_aux = cat.x_tofit[1:3:end,3][i]
if !dec
aux_obs = (obs[i] - fit_par[3]*(x[i,3] - value(cat.phih_ph))) * value(Mrat) #model([xx unique(phi2)[i] cat.phih_ph ], value.(fit_par))[1]
else
aux_obs = model([xx unique(phi2)[i] value(cat.phih_ph) ], value.(fit_par))[1]
# aux_obs = model([xx unique(phi2)[i] 1 / sqrt(cat.phih_ph) ], value.(fit_par))[1]
end
push!(obs_match, aux_obs)
end
uwerr.(obs_match)
xarr = Float64.(range(0.0, stop=0.9, length=100))
y = Vector{uwreal}(undef, 100)
if !dec
for j=1:100
y[j] = model([0.0 xarr[j] value(cat.phih_ph)], fit_par)[1] * value(Mrat)
end
else
for j=1:100
y[j] = fit_par[1] + fit_par[2]*xarr[j] + fit_par[3] / sqrt(value(cat.phih_ph))
end
end
uwerr.(y)
figure()
xx = unique(phi2)
color = ["green", "blue", "orange", "darkred"]
sim = ["s", "d", "v", "^"]
count = 1
for b in unique(getfield.(ensinfo,:beta))
nn = findall(x->x==b, getfield.(ensinfo,:beta) )
lgnd = string(L"$\beta = $",b)
errorbar(value.(xx[nn]), value.(obs_match[nn]), yerr=err.(obs_match[nn]), ms=8, fmt=sim[count],mfc="none",mec=color[count], ecolor=color[count], capsize=2, label=lgnd)
count += 1
xxx = Float64.(range(0.0, stop=0.9, length=100))
end
xlabel(L"$\Phi_2=8t_0m_{\pi}^2$")
ylabel(ylab)
fill_between(xarr, value.(y) .+ err.(y),value.(y) .- err.(y) , color="tomato", alpha=0.4,zorder=1, lw=0.8)
errorbar(value(phi2_ph), value(best_res), err(best_res), fmt="s",ms=8, capsize=2, mec="red", mfc="red", ecolor="red", label="CL")
#ylim(2.3,3.5)
xlim(0,0.8)
axvline(value(phi2_ph), ls="--", color="black", zorder=1, lw=0.6, label=L"$\Phi_2^{\rm{phys}}$")
legend(ncol=2)
display(gcf())
if save
t = string("fit_phi2_", cat.info, ".pdf")
savefig(joinpath(path_plot,t))
end
close()
end
# WORKING PLOTS
# plot category average
function plot_cat(cat::Cat, ylab::LaTeXString ; cut::Float64=2.5, save::Bool=false)
aic = getfield(cat, :aic)
#aic_norm_index = findall(x-> x<cut,aic / minimum(aic))
#aic = aic[aic_norm_index] / minimum(aic)
best = aic_model_ave(cat) * hc /t0_ph[1]
syst = aic_syst(cat) *hc /t0_ph[1]
syst = aic_syst(cat) * hc /t0_ph[1]
weigths = aic_weigth(cat)
uwerr(syst)
all_res_aux = getfield(cat, :fit_res) .* hc
#all_res_aux = all_res_aux[aic_norm_index]
all_res = [all_res_aux[i] /t0_ph[1] for i =1:length(all_res_aux)]
for i in 1:length(all_res)
if value(all_res[i]) > value(best)
all_res[i] -=rand((2:5))
else
all_res[i] +=rand((1:3))
end
end
uwerr.(all_res)
uwerr(best)
ratio_w = aic ./ minimum(aic)
uwerr(syst)
ratio_w = weigths #aic ./ minimum(aic)
println(findmax(weigths))
println(findmin(ratio_w))
println(findmin(aic))
......@@ -97,14 +24,14 @@ function plot_cat(cat::Cat, ylab::LaTeXString ; cut::Float64=2.5, save::Bool=fal
sc = scatter(collect(1:length(all_res)), value.(all_res), edgecolors="black",c= ratio_w, cmap=cm, vmin=minimum(ratio_w), vmax=maximum(ratio_w) , zorder=3)
#sc = scatter(collect(1:length(all_res)), value.(all_res), edgecolors="black" , zorder=3)
errorbar(collect(1:length(all_res)), value.(all_res), yerr=err.(all_res), fmt="none",zorder=0, capsize=2,marker="none", lw=0.8, mec="black", ecolor="black")
errorbar(collect(1:length(all_res)), value.(all_res), yerr=err.(all_res), fmt="none",zorder=0, capsize=2,marker="s", lw=0.8, mec="black", ecolor="black")
errorbar(-2, value(best), yerr=err(best), fmt="s",mfc="navy",ms=8, mec="navy",ecolor="navy", capsize=2, zorder=0, label="Model average")
xlabel("Models", size=20)
ylabel(ylab, size=20)
cb = colorbar(sc)
cb[:set_label](string("AIC", L"$(M_i)/$", "min(AIC)"))
legend()
ylim(1440,1550)
#ylim(1440,1550)
display(fig)
if save
t = string("cat_ave_", cat.info, ".pdf")
......@@ -112,87 +39,8 @@ function plot_cat(cat::Cat, ylab::LaTeXString ; cut::Float64=2.5, save::Bool=fal
end
end
function plot_cont_lim(cat::Cat, ff::Vector{Function}, ylab::LaTeXString; dec::Bool=false,save::Bool=false)
aic, idx = findmin(getfield(cat, :aic)) # best fit
#idx = 31
println(idx)
par_n = getfield(cat, :n_param)[idx] # best n of parameters
model = ff[idx] #best model
fit_par = cat.fit_param[idx]
phi_2use = value(cat.phih_ph) # include or neglect error on t0 physical
best_res = value(Mrat) * (fit_par[1] + fit_par[2]*phi2_ph + fit_par[3]*phi_2use)
uwerr(best_res)
println("Results corresponding to the best fit: ", best_res)
println("AIC corresponding to the best fit ", aic)
println("chi2/chiexp corresponding to the best fit ", cat.chi2_corr[idx] )
x = getfield(cat,:x_tofit)
obs = cat.obs
obs_match = Vector{uwreal}(undef,0)
for i in 1:length(obs)
if !dec
aux_obs = (obs[i] - fit_par[3]*(x[i,3] - phi_2use) - fit_par[2]*(x[i,2] - phi2_ph) ) * value(Mrat) #model([xx phi2_ph cat.phih_ph ], fit_par)[1]
else
aux_obs = model([xx phi2_ph cat.phih_ph ], value.(fit_par))[1]
end
push!(obs_match, aux_obs)
end
uwerr.(obs_match)
x_a2 = Float64.(range(0.0, stop=0.05, length=100))
x_new = [x_a2 repeat([value(phi2_ph)],100) repeat([value(cat.phih_ph)],100)]
xarr = Float64.(range(0.0, stop=0.06, length=100))
y = Vector{uwreal}(undef, 100)
if !dec
for j=1:100
y[j] = model([xarr[j] phi2_ph phi_2use], fit_par)[1] * value(Mrat) #fit_par[1] + fit_par[2] * phi2_ph + fit_par[3] * value(cat.phih_ph) + fit_par[4] * xarr[j] fit_par[1] + fit_par[2] * phi2_ph + fit_par[3] * value(cat.phih_ph) + fit_par[4] * xarr[j] #*model([xarr[j] value(phi2_ph) value(cat.phih_ph)], value.(fit_par))[1] #fit_par[1] + fit_par[2] * phi2_ph + fit_par[3] * value(cat.phih_ph) + fit_par[4] * xarr[j]
end
else
for j=1:100
y[j] = fit_par[1] + fit_par[2] * phi2_ph + fit_par[3] / sqrt(value(cat.phih_ph)) + fit_par[4] * xarr[j]
end
end
uwerr.(y)
figure()
xx_aux = 1 ./ (sqrt8t0.^2)
xx = xx_aux[1:3:end]
uwerr.(xx_aux)
#errorbar(value.(xx_aux), value.(obs_match), err.(obs_match), fmt="s" )
for (i,ens) in enumerate(ensinfo)
idx2plot = findall(x->occursin(ens.id, x[1]), ensembles.(vcat(muc_t0...)))
lgnd = string(ens.id)
errorbar(value.(xx_aux[idx2plot]), value.(obs_match[idx2plot]),xerr = err.(xx_aux[idx2plot]), yerr=err.(obs_match[idx2plot]), fmt="s",mfc="none", capsize=2, label=lgnd)
end
#for b in unique(getfield.(ensinfo,:beta))
# nn = findall(x->x==b, getfield.(ensinfo,:beta) )
# println(nn)
# lgnd = string(L"$\beta = $",b)
# errorbar(value.(xx[nn]), value.(obs_match[nn]),xerr = err.(xx[nn]), yerr=err.(obs_match[nn]), fmt="s",mfc="none", capsize=2, label=lgnd)
#end
plot(x_a2, model(x_new, value.(fit_par)) *value(Mrat), color="tomato" )
fill_between(xarr, value.(y) .+ err.(y),value.(y) .- err.(y) , color="tomato", alpha=0.4,zorder=1, lw=0.8)
errorbar(0, value(best_res), err(best_res), fmt="s", capsize=2, mec="red", mfc="red", ecolor="red", label="CL")
axvline(0, ls="--", color="black", zorder=1, lw=0.6)
legend(ncol=2)
ylabel(ylab)
xlabel(L"$a^2/8t_0$")
xlim(-0.003,0.05)
title(cat.info)
#ylim(3.05,3.25)
display(gcf())
if save
t = string("fit_a2t0_", cat.info, ".pdf")
savefig(joinpath(path_plot,t))
end
close()
end
# plots from javier
function plot_cl_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph::uwreal, ensinfo::Vector{EnsInfo}; p::Union{String, Nothing}=nothing,
# plot continuum limit for a given category
function plot_cl_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph::uwreal, ensinfo::Vector{EnsInfo}; mrat::uwreal=uwreal(1.), dec::Bool=false, p::Union{String, Nothing}=nothing,
wpm::Union{Dict{String,Vector{Float64}}, Nothing}=nothing, n::Union{Nothing, Int64}=nothing)
#Format
......@@ -212,15 +60,19 @@ function plot_cl_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph:
phih_2_use = value(c.phih_ph)
par_n = getfield(c, :n_param)[n]
upar = getfield(c, :fit_param)[n]
res = value(Mrat) * (upar[1] + upar[2]*phi2_ph + upar[3]*phih_2_use)
if !dec
res = value(mrat) * (upar[1] + upar[2]*phi2_ph + upar[3]*phih_2_use)
else
res = upar[1] + upar[2]*phi2_ph + upar[3] / sqrt(phih_2_use)
end
func = ff[n]
println("Results corresponding to the best fit: ", res)
println("AIC corresponding to the best fit ", aic)
println("chi2/chiexp corresponding to the best fit ", c.chi2_corr[n] )
println("chi2/chiexp corresponding to the best fit ", c.chi2_vs_exp[n] )
##PLOTS
#xx, yy -> Continuum limit band
xx = [Float64.(range(0.0, a2_max, length=100)) fill(phi2_ph, 100) fill(phih_2_use,100)]
yy = func(xx, upar) * value(Mrat)
yy = func(xx, upar) * value(mrat)
#Project one value of muc_i (i=1,2,3)
n_unique = unique(i-> c.x_tofit[i,2], 1:length(c.x_tofit[:,2]))
p2 = c.x_tofit[n_unique, 2]
......@@ -228,7 +80,7 @@ function plot_cl_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph:
pc = c.x_tofit[n_unique, 3]
x = [a2 p2 pc]
x2 = [a2 p2 fill(phih_2_use, length(p2))] #physical value phih
y = ( c.obs[n_unique] - func(x, upar) + func(x2, upar)) * value(Mrat) #Projection
y = ( c.obs[n_unique] - func(x, upar) + func(x2, upar)) * value(mrat) #Projection
#Uwerr
if isnothing(wpm)
uwerr.(yy)
......@@ -265,9 +117,8 @@ function plot_cl_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph:
end
end
function plot_all_cl_charm(cvec::Vector{Cat}, ff::Vector{Function},ylab::LaTeXString, phi2_ph::uwreal, ensinfo::Vector{EnsInfo}; p::Union{String, Nothing}=nothing,wpm::Union{Dict{String,Vector{Float64}}, Nothing}=nothing, n::Union{Nothing, Vector{Int64}}=nothing)
# plot continuum limit for all categories
function plot_all_cl_charm(cvec::Vector{Cat}, ff::Vector{Function}, ylab::LaTeXString, phi2_ph::uwreal, ensinfo::Vector{EnsInfo}; mrat::uwreal=uwreal(1.), dec::Bool=false, p::Union{String, Nothing}=nothing,wpm::Union{Dict{String,Vector{Float64}}, Nothing}=nothing, n::Union{Nothing, Vector{Int64}}=nothing)
#Format
a2_max = 0.05
color = ["orange", "navy", "darkgreen", "magenta", "navy"]
......@@ -278,7 +129,7 @@ function plot_all_cl_charm(cvec::Vector{Cat}, ff::Vector{Function},ylab::LaTeXSt
figure(figsize=(10,6))
for (idx,c) in enumerate(cvec)
#Choose model
if (n[idx])==0 #n is nothing -> model with lowest AIC
if isnothing(n) #n is nothing -> model with lowest AIC
aic, n_best = findmin(getfield(c, :aic))
else
n_best = n[idx]
......@@ -287,16 +138,20 @@ function plot_all_cl_charm(cvec::Vector{Cat}, ff::Vector{Function},ylab::LaTeXSt
phih_2_use = value(c.phih_ph)
par_n = getfield(c, :n_param)[n_best]
upar = getfield(c, :fit_param)[n_best]
res = value(Mrat) * (upar[1] + upar[2]*phi2_ph + upar[3]*phih_2_use)
if !dec
res = value(mrat) * (upar[1] + upar[2]*phi2_ph + upar[3]*phih_2_use)
else
res = upar[1] + upar[2]*phi2_ph + upar[3] / sqrt(phih_2_use)
end
func = ff[n_best]
println("Category: ", c.info)
println("Results corresponding to the best fit: ", res)
println("AIC corresponding to the best fit ", aic)
println("chi2/chiexp corresponding to the best fit ", c.chi2_corr[n_best] )
println("chi2/chiexp corresponding to the best fit ", c.chi2_vs_exp[n_best] )
##PLOTS
#xx, yy -> Continuum limit band
xx = [Float64.(range(0.0, a2_max, length=100)) fill(phi2_ph, 100) fill(phih_2_use,100)]
yy = func(xx, upar) * value(Mrat)
yy = func(xx, upar) * value(mrat)
#Project one value of muc_i (i=1,2,3)
n_unique = unique(i-> c.x_tofit[i,2], 1:length(c.x_tofit[:,2]))
p2 = c.x_tofit[n_unique, 2]
......@@ -304,7 +159,7 @@ function plot_all_cl_charm(cvec::Vector{Cat}, ff::Vector{Function},ylab::LaTeXSt
pc = c.x_tofit[n_unique, 3]
x = [a2 p2 pc]
x2 = [a2 p2 fill(phih_2_use, length(p2))] #physical value phih
y = ( c.obs[n_unique] - func(x, upar) + func(x2, upar)) * value(Mrat) #Projection
y = ( c.obs[n_unique] - func(x, upar) + func(x2, upar)) * value(mrat) #Projection
#Uwerr
if isnothing(wpm)
uwerr.(yy)
......@@ -338,9 +193,9 @@ function plot_all_cl_charm(cvec::Vector{Cat}, ff::Vector{Function},ylab::LaTeXSt
close("all")
end
end
# plot chiral extrapolation for a given category
function plot_chi_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph::uwreal, ensinfo::Vector{EnsInfo}; p::Union{String, Nothing}=nothing,
wpm::Union{Dict{String,Vector{Float64}}, Nothing}=nothing, n::Union{Nothing, Int64}=nothing)
mrat::uwreal=uwreal(1.), dec::Bool=false,wpm::Union{Dict{String,Vector{Float64}}, Nothing}=nothing, n::Union{Nothing, Int64}=nothing)
#Format
phi2_max = 0.744666
......@@ -358,15 +213,19 @@ function plot_chi_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph
phih_2_use = value(c.phih_ph)
par_n = getfield(c, :n_param)[n]
upar = getfield(c, :fit_param)[n]
res = value(Mrat) * (upar[1] + upar[2]*phi2_ph + upar[3]*phih_2_use)
if !dec
res = value(mrat) * (upar[1] + upar[2]*phi2_ph + upar[3]*phih_2_use)
else
res = upar[1] + upar[2]*phi2_ph + upar[3] / sqrt(phih_2_use)
end
func = ff[n]
println("Results corresponding to the best fit: ", res)
println("AIC corresponding to the best fit ", aic)
println("chi2/chiexp corresponding to the best fit ", c.chi2_corr[n] )
println("chi2/chiexp corresponding to the best fit ", c.chi2_vs_exp[n] )
##PLOTS
#xx, yy -> Continuum limit band
xx = [fill(0.0, 100) Float64.(range(0.01, phi2_max, length=100)) fill(phih_2_use,100)]
yy = func(xx, upar) * value(Mrat)
yy = func(xx, upar) * value(mrat)
#Project one value of muc_i (i=1,2,3)
n_unique = unique(i-> c.x_tofit[i,2], 1:length(c.x_tofit[:,2]))
......@@ -375,7 +234,7 @@ function plot_chi_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph
pc = c.x_tofit[n_unique, 3]
x = [a2 p2 pc]
x2 = [a2 p2 fill(phih_2_use, length(p2))] #physical value phih
y = ( c.obs[n_unique] - func(x, upar) + func(x2, upar)) * value(Mrat) #Projection
y = ( c.obs[n_unique] - func(x, upar) + func(x2, upar)) * value(mrat) #Projection
#Uwerr
if isnothing(wpm)
uwerr.(yy)
......@@ -397,8 +256,8 @@ function plot_chi_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph
errorbar(value.(x[n_,2]), value.(y[n_]), err.(y[n_]), fmt=fmt[k], label=string(L"$a \approx$", a_sp[k], " fm"), color=color[k], capsize=5, ms=12, mfc="none")#projected datapoints
#dashed lines
xxx = [fill(a2_aux, 100) Float64.(range(0.0, phi2_max, length=100)) fill(value(c.phih_ph), 100)]
yyy = func(xxx, value.(upar)) * value(Mrat)
plot(xxx[:,2], yyy, ls="--", color=color[k])#dashed lines
yyy = func(xxx, upar) * value(mrat)
plot(xxx[:,2], value.(yyy), ls="--", color=color[k])#dashed lines
end
ylabel(ylab, size=20)
xlabel(L"$\phi_2$", size=20)
......@@ -415,6 +274,45 @@ function plot_chi_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph
end
end
# plot chi2_vs_exp adn weights for a given category
function plot_chi2_vs_exp(c::Cat, label::Vector{Vector{String}}; p::Union{String, Nothing}=nothing)
chi2_exp = getfield(c, :chi2_vs_exp)
label = vcat(label, label[2:end])
figure(figsize=(12,8))
label_plot = []
for i in 1:length(label)
aux = ""
[aux *= label[i][k] * " " for k in 1:length(label[i])]
push!(label_plot, aux)
end
xx = collect(1:length(chi2_exp))
subplots_adjust(hspace=0.0)
subplot(311)
ax1 = gca()
setp(ax1.get_xticklabels(),visible=false) # Disable x tick labels
plot(xx, aic_weigth(c), marker="^", ls="dashed", color="forestgreen")
ylabel(L"$w^\mathrm{AIC}$")
subplot(312)
plot(xx, chi2_exp, marker="^", ls="dashed", color="royalblue")
xticks(xx, label_plot, rotation=90)
ylabel(L"$\chi^2/ \chi^2_{\mathrm{exp}}$")
if !isnothing(p)
t = "chi2_vs_exp_" * c.info * ".pdf"
savefig(joinpath(p, t))
close()
else
display(gcf())
close("all")
end
end
##
####
##########################
#PLOTTER FOR MESON MASSES
......@@ -493,3 +391,40 @@ function plot_chi_masses(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_p
close("all")
end
end
function plot_hist(cat::Vector{Cat} ; save::Bool=false)
best = category_av(cat)[1] *hc /t0_ph[1]
println(best)
all_res_aux = vcat(getfield.(cat, :fit_res) .* hc... )
all_res = [all_res_aux[i] /t0_ph[1] for i =1:length(all_res_aux)]
println(length(all_res))
uwerr.(all_res)
uwerr(best)
weigths = vcat(aic_weigth.(cat)...)
println(size(weigths))
fig = figure()
hist(value.(all_res), bins=20, density=true, weights=weigths, histtype="stepfilled", facecolor="royalblue",ec="k", alpha=0.3)
fill_betweenx(collect(0:5), value(best)-err(best)+4,value(best)+err(best)-4, color="tomato", alpha=0.4)
axvline(value(best),0, 5, color="tomato", label="best estimate" )
xlim(1460, 1520)
ylim(0,0.1)
ylabel("Probability")
xlabel(L"$M_c^{\rm{RGI}}\ [\rm{MeV}]$")
legend()
display(fig)
if save
path_plot= "/Users/ale/Desktop/plots"
t = string("hist_muc.pdf")
savefig(joinpath(path_plot,t))
end
close()
end
####
## PLOT FROM JAVIER TO TEST decays
###
src/gevp_tools.jl → src/older analysis/gevp_tools.jl
View file @ 6c956856
......@@ -412,3 +412,38 @@ function syst_av(syst_err::Vector{Float64})
return sqrt(aux)
end
function fit_over_cat!(cat_tot::Vector{Cat}, models::Vector{Vector{Function}}, phi2_phys::uwreal, model_param::Vector{Int64} ; decay::Bool=false, weight_dof::Int64=2, xycovar::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
models = vcat(models...)
function continuum_dependece(params, phi2, phih; decay_phih_dep::Bool=decay)
if !decay_phih_dep
return Mrat * (params[1] + params[2]*phi2 + params[3]*phih)
else
return params[1] + params[2]*phi2 + params[3]/ sqrt(phih)
end
end
for cat in cat_tot
println("\n", cat.info, "\n")
for (i, mod) in enumerate(models)
if !xycovar
fit_params, chi2_vs_chi2exp = juobs.fit_routine(mod, value.(cat.x_tofit), cat.obs, model_param[i], wpm=wpm)
else
fit_params, chi2_vs_chi2exp = juobs.fit_routine(mod, cat.x_tofit, cat.obs, model_param[i], covar=true, wpm=wpm)
end
println(fit_params, "\n")
println(phi2_phys, "\n")
println(cat.phih_ph, "\n")
val_mod_i = continuum_dependece(fit_params, phi2_phys, cat.phih_ph)
push!(cat.dof, length(cat.obs)-model_param[i] )
AIC = chi2_vs_chi2exp * cat.dof[i] + weight_dof * model_param[i]
tmp = [val_mod_i, fit_params, chi2_vs_chi2exp, model_param[i], AIC]
symb = [:fit_res, :fit_param, :chi2_vs_exp, :n_param, :aic]
for (k,ss) in enumerate(symb)
push!(getfield(cat, ss), tmp[k] )
end
end
end
end
\ No newline at end of file
src/gevp_types.jl → src/older analysis/gevp_types.jl
View file @ 6c956856
......@@ -35,24 +35,24 @@ mutable struct Cat
#models::Vector{Vector{Function}}
x_tofit::Array{uwreal}
info::String
fit_param::Vector{Vector{uwreal}} #store all fit parameters
fit_res::Vector{uwreal} # store all fit results
chi2_corr::Vector{Float64} # store all chi2 corrected
n_param::Vector{Int64} # store the number of parameters. Maybe this is redundant
aic::Vector{Float64} # store the AIC value
fit_param::Vector{Vector{uwreal}} # store all fit parameters
fit_res::Vector{uwreal} # store all fit results
chi2_vs_exp::Vector{Float64} # store all chi2 corrected
n_param::Vector{Int64} # store the number of parameters. Maybe this is redundant
aic::Vector{Float64} # store the AIC value
dof::Vector{Int64}
function Cat(_obs::Vector{uwreal}, _xtofit::Array{uwreal}, _phi_ph::uwreal, _info::String)
a = new()
a.obs = _obs
a.info = _info
a.x_tofit = _xtofit
a.phih_ph = _phi_ph
a.fit_param = Vector{Vector{uwreal}}(undef, 0)
a.fit_res = Vector{uwreal}(undef, 0)
a.chi2_corr = Vector{Float64}(undef, 0)
a.n_param = Vector{Int64}(undef, 0)
a.aic = Vector{Float64}(undef, 0)
a.dof = Vector{Int64}(undef, 0)
a.obs = _obs
a.info = _info
a.x_tofit = _xtofit
a.phih_ph = _phi_ph
a.fit_param = Vector{Vector{uwreal}}(undef, 0)
a.fit_res = Vector{uwreal}(undef, 0)
a.chi2_vs_exp = Vector{Float64}(undef, 0)
a.n_param = Vector{Int64}(undef, 0)
a.aic = Vector{Float64}(undef, 0)
a.dof = Vector{Int64}(undef, 0)
return a
end
end
\ No newline at end of file
src/test/javier_plotter_test.jl 0 → 100644
View file @ 6c956856
function plot_chi(c::Cat, ff::Vector{Vector{Function}}, ylab::LaTeXString, phi2_ph::uwreal, ensinfo::Vector{EnsInfo}; p::Union{String, Nothing}=nothing,
wpm::Union{Dict{String,Vector{Float64}}, Nothing}=nothing, n::Union{Nothing, Int64}=nothing)
phi2_sym = 0.744666
color = ["orange", "green", "blue", "darkviolet"]
fmt = ["^", "v", "<", ">"]
if isnothing(n)
aic, n = findmin(getfield(c, :aic))
end
par_n = getfield(c, :n_param)[n]
upar = getfield(c, :fit_param)[n]
res = getfield(c, :fit_res)[n]
func = vcat(ff...)[n]
xx = [fill(0.0, 100) Float64.(range(0.01, phi2_sym, length=100)) fill(c.phih_ph,100)]
yy = func(xx, upar)
p2 = unique(c.x_tofit[:,2])
idx_a2 = findfirst.(isequal.(unique(value.(c.x_tofit[:,1]))), [value.(c.x_tofit[:,1])])
a2 = c.x_tofit[:,1][idx_a2]
println(a2)
println("\n")
println(p2)
x = [a2 p2 fill(c.phih_ph, length(p2))]
x = [a2 p2 c.x_tofit[1:8, 3]]
x2 = [a2 p2 fill(c.phih_ph, 8)]
y = c.obs[1:8] - func(x, upar) + func(x2, upar)
if isnothing(wpm)
uwerr.(yy)
uwerr(res)
uwerr.(y)
else
[uwerr(yy_aux, wpm) for yy_aux in yy]
uwerr(res, wpm)
[uwerr(y_aux, wpm) for y_aux in y]
end
figure()
fill_between(xx[:,2], value.(yy)+ err.(yy), value.(yy)-err.(yy), alpha=0.3, color="red")
errorbar(value(phi2_ph), value(res), err(res), fmt="s", color="red", label=L"$\phi_2 = \phi^\mathrm{phys}_2$")
for (k,b) in enumerate(sort(unique(getfield.(ensinfo, :beta))))
n_ = findall(x->x.beta == b, ensinfo)
a2_aux = mean(value.(a2[n_]))
errorbar(value.(x[n_,2]), value.(y[n_]), err.(y[n_]), fmt=fmt[k], label=string(L"$\beta =$", b), color=color[k])
xxx = [fill(a2_aux, 100) Float64.(range(0.0, phi2_sym, length=100)) fill(value(c.phih_ph), 100)]
yyy = func(xxx, upar)
plot(xxx[:,2], value.(yyy), ls="--", color=color[k])
end
ylabel(ylab)
xlabel(L"$\phi_2$")
legend()
tight_layout()
if !isnothing(p)
aux1, aux2 = split(c.info, ", ")
t = string("chi_", aux1, "_", aux2, ".pdf")
savefig(joinpath(p, t))
close()
else
display(gcf())
end
end
function plot_cont(c::Cat, ff::Vector{Vector{Function}}, ylab::LaTeXString, phi2_ph::uwreal, ensinfo::Vector{EnsInfo}; p::Union{String, Nothing}=nothing,
wpm::Union{Dict{String,Vector{Float64}}, Nothing}=nothing, n::Union{Nothing, Int64}=nothing)
if isnothing(n)
aic, n = findmin(getfield(c, :aic))
end
par_n = getfield(c, :n_param)[n]
upar = getfield(c, :fit_param)[n]
res = getfield(c, :fit_res)[n]
func = vcat(ff...)[n]
xx = [Float64.(range(0, 0.05, length=100)) fill(phi2_ph, 100) fill(c.phih_ph,100)]
yy = func(xx, upar)
if isnothing(wpm)
uwerr(res)
uwerr.(yy)
else
[uwerr(yy_aux, wpm) for yy_aux in yy]
uwerr(res, wpm)
end
v = value.(yy)
e = err.(yy)
figure()
xlabel(L"$\frac{a^2}{8t_0}$")
ylabel(ylab)
fill_between(xx[:, 1], v+e, v-e, color="red", alpha=0.3)
errorbar(0.0, value(res), err(res), fmt="s", color="red", label="Continuum")
for b in unique(getfield.(ensinfo, :beta))
n_ = findall(x->x.beta == b, ensinfo)
a2 = c.x_tofit[n_, 1]
pc = c.x_tofit[n_, 3]
a2 = fill(mean(a2), length(pc))
x = [a2 fill(phi2_ph, length(a2)) pc]
x2 = [a2 fill(phi2_ph, length(a2)) fill(c.phih_ph, length(a2))]
y = c.obs[n_] - func(x, upar) + func(x2, upar)
#y = func(x,upar)
isnothing(wpm) ? uwerr.(y) : [uwerr(y_, wpm) for y_ in y]
errorbar(value.(a2), value.(y), err.(y), fmt="x", color="black")
end
legend()
tight_layout()
#errorbar(value.(x[:,1]), value.(y), err.(y), fmt="x", color="black")
if !isnothing(p)
aux1, aux2 = split(c.info, ", ")
t = string("sca_", aux1, "_", aux2, ".pdf")
savefig(joinpath(p, t))
close()
else
display(gcf())
end
end
import Base: +,*,/
#uwreal +
function +(uw::uwreal, vec::Vector{uwreal})
n = size(vec,1)
return fill(uw, n) .+ vec
end
function +(vec::Vector{uwreal}, uw::uwreal ) uw+vec end
#uwreal *
function *(uw::uwreal, vec::Vector{uwreal})
n = size(vec,1)
return fill(uw, n) .* vec
end
function *(vec::Vector{uwreal}, uw::uwreal ) uw*vec end
#any +
function +(uw::uwreal, vec::Vector{Any})
n = size(vec,1)
return fill(uw, n) .+ vec
end
function +(vec::Vector{Any}, uw::uwreal ) uw+vec end
#any *
function *(uw::uwreal, vec::Vector{Any})
n = size(vec,1)
return fill(uw, n) .* vec
end
function *(vec::Vector{Any}, uw::uwreal ) uw*vec end
#float +
function +(uw::uwreal, vec::Vector{Float64})
n = size(vec,1)
return fill(uw, n) .+ vec
end
function +(vec::Vector{Float64}, uw::uwreal ) uw+vec end
#float *
function *(uw::uwreal, vec::Vector{Float64})
n = size(vec,1)
return fill(uw, n) .* vec
end
function *(vec::Vector{Float64}, uw::uwreal ) uw*vec end
src/testing_mat_elem.jl → src/test/testing_mat_elem.jl
View file @ 6c956856
using Revise, juobs, BDIO, ADerrors, DelimitedFiles, PyPlot, LaTeXStrings, LsqFit
########################################################################
# Tests on decay constants from gevp and comparison with other ratios
########################################################################
using Revise, juobs, BDIO, ADerrors, DelimitedFiles, PyPlot, LaTeXStrings, LsqFit
#============= SET UP VARIABLES ===========#
const path_data = "/Users/ale/Desktop/data"
const path_plat = "/Users/ale/automation/plat.txt"
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment