Implemented lagrange 2pt and IO. Tests needed.

b2811423 · Fernando Pérez Panadero · 6dcebc0d · b2811423 · b2811423 · b2811423
Commit b2811423 authored Oct 31, 2024 by Fernando Pérez Panadero
Hide whitespace changes
Inline Side-by-side

Showing with 133 additions and 81 deletions

main.jl main.jl +2 -2

src/io.jl src/io.jl +5 -3

src/meas.jl src/meas.jl +20 -20

src/utils.jl src/utils.jl +106 -56

No files found.
--- a/main.jl
+++ b/main.jl
@@ -29,9 +29,9 @@ while RUN_ON
    @timeit "HMC" Gauge_update()
-    @timeit "Two point function" begin
+    @timeit "Two point functions" begin
        Frontflow_pt()
-        two_pt_lagrange()
+        Two_pt_lagrange()
    end
    @timeit "One point function" Backflow_pt()

--- a/src/io.jl
+++ b/src/io.jl
@@ -135,6 +135,10 @@ function save_data()
        BDIO_write!(fb, [convert(Int32,params["Frontflow"]["nsteps"])])
        BDIO_write!(fb, [convert(Int32,params["Backflow"]["N_noise"])])
        BDIO_write!(fb, [convert(Int32,length(flow_times))])
+        BDIO_write!(fb, [convert(Int32,params["Frontflow"]["t_zero"])])
+        BDIO_write!(fb, [convert(Int32,params["Frontflow"]["epsilon"])])
+        BDIO_write!(fb, [convert(Int32,params["Frontflow"]["tsource"])])
+        BDIO_write!(fb, [convert(Int32,params["Backflow"]["tsource"])])
        length(flow_times) > 0 ? BDIO_write!(fb, flow_times) : nothing
        BDIO_write_hash!(fb)
@@ -151,7 +155,7 @@ function save_data()
            BDIO_write!(fb,pp_corr_t[:,noi,fl])
            BDIO_write!(fb,ap_corr_t[:,noi,fl])
-            BDIO_write!(fb,pp_hat_t[:,noi,fl])
+            BDIO_write!(fb,pphat_t[:,noi,fl])
            BDIO_write!(fb,pptilde_t[:,noi,fl])
        end
    end
@@ -172,8 +176,6 @@ function save_data()
            BDIO_write!(fb,ChiDchi[:,noi,fl])
            BDIO_write!(fb,ChiDchi_cfl[:,noi,fl])
-            BDIO_write!(fb,[phat_t[noi,fl]])
        end
    end

--- a/src/meas.jl
+++ b/src/meas.jl
@@ -69,7 +69,7 @@ function Frontflow_pt() # Will be Frontflow
        pp_corr_t0[:,noi] .= zero(Float64)
        ap_corr_t0[:,noi] .= zero(ComplexF64)
-        for t in 1:lp.iL[4]
+        @timeit "Volume sum" for t in 1:lp.iL[4]
            for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
                b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
                pp_corr_t0[t,noi] += pp_density[b,r]
@@ -94,7 +94,7 @@ function Frontflow_pt() # Will be Frontflow
            ap_density .= Array(dot.(psi,dmul.(Gamma{4},psi)))
            pp_corr_t[:,noi,tstep] .= zero(Float64)
            ap_corr_t[:,noi,tstep] .= zero(ComplexF64)
-            for t in 1:lp.iL[4]
+            @timeit "Volume sum" for t in 1:lp.iL[4]
                for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
                    b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
                    pp_corr_t[t,noi,tstep] += pp_density[b,r]
@@ -109,7 +109,7 @@ function Frontflow_pt() # Will be Frontflow
        ap_density .= Array(dot.(psi,dmul.(Gamma{4},psi)))
        pp_corr_t[:,noi,end] .= zero(Float64)
        ap_corr_t[:,noi,end] .= zero(ComplexF64)
-        for t in 1:lp.iL[4]
+        @timeit "Volume sum" for t in 1:lp.iL[4]
            for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
                b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
                pp_corr_t[t,noi,end] += pp_density[b,r]
@@ -172,7 +172,7 @@ function Backflow_pt() # Will be backflow
            Quark_cond[:,noi,fl] .= zero(ComplexF64)
            ap_density .= Array(dot(dws.st[b,r],psi[b,r]))
-            for t in 1:lp.iL[4]
+            @timeit "Volume sum" for t in 1:lp.iL[4]
                for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
                    b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
                    Quark_cond[t,noi,fl] += ap_density[b,r]
@@ -181,7 +181,7 @@ function Backflow_pt() # Will be backflow
            Quark_cond_cfl[:,noi,fl] .= zero(ComplexF64)
            pp_density .= Array(norm2.(dws.st))
-            for t in 1:lp.iL[4]
+            @timeit "Volume sum" for t in 1:lp.iL[4]
                for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
                    b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
                    Quark_cond_cfl[t,noi,fl] += -pp_density[b,r]
@@ -192,7 +192,7 @@ function Backflow_pt() # Will be backflow
            ap_density .= Array(dot.(psi,dmul.(Gamma{4},psi)))
            pp_corr_tfl[:,noi,fl] .= zero(Float64)
            ap_corr_tfl[:,noi,fl] .= zero(ComplexF64)
-            for t in 1:lp.iL[4]
+            @timeit "Volume sum" for t in 1:lp.iL[4]
                for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
                    b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
                    pp_corr_tfl[t,noi,fl] += pp_density[b,r]
@@ -211,7 +211,7 @@ function Backflow_pt() # Will be backflow
                ChiDchi[:,noi,fl] .= zero(Float64)
                ap_density .= Array(dot(dws.st[b,r],dws.sp[b,r]))
-                for t in 1:lp.iL[4]
+                @timeit "Volume sum" for t in 1:lp.iL[4]
                    for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
                        b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
                        ChiDchi[t,noi,fl] += ap_density[b,r]
@@ -219,7 +219,7 @@ function Backflow_pt() # Will be backflow
                end
                Dw!(dws.sp,U,dws.st,dpar,dws,lp)
                ap_density .= Array(dot(dws.sp[b,r],psi[b,r]))
-                for t in 1:lp.iL[4]
+                @timeit "Volume sum" for t in 1:lp.iL[4]
                    for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
                        b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
                        ChiDchi[t,noi,fl] += -ap_density[b,r]
@@ -228,7 +228,7 @@ function Backflow_pt() # Will be backflow
                Quark_cond2[:,noi,fl] .= zero(ComplexF64)
                ap_density .= Array(dot(dws.st[b,r],psi[b,r]))
-                for t in 1:lp.iL[4]
+                @timeit "Volume sum" for t in 1:lp.iL[4]
                    for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
                        b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
                        Quark_cond2[t,noi,fl] += ap_density[b,r]
@@ -243,7 +243,7 @@ function Backflow_pt() # Will be backflow
                ChiDchi_cfl[:,noi,fl] .= zero(ComplexF64)
                ap_density .= Array(dot(dws.st[b,r],dws.sp[b,r]))
-                for t in 1:lp.iL[4]
+                @timeit "Volume sum" for t in 1:lp.iL[4]
                    for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
                        b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
                        ChiDchi_cfl[t,noi,fl] += ap_density[b,r]
@@ -251,7 +251,7 @@ function Backflow_pt() # Will be backflow
                end
                Dw!(dws.sp,U,dws.st,dpar,dws,lp)
                ap_density .= Array(dot(dws.sp[b,r],psi[b,r]))
-                for t in 1:lp.iL[4]
+                @timeit "Volume sum" for t in 1:lp.iL[4]
                    for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
                        b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
                        ChiDchi_cfl[t,noi,fl] += -ap_density[b,r]
@@ -260,7 +260,7 @@ function Backflow_pt() # Will be backflow
                Quark_cond2_cfl[:,noi,fl] .= zero(ComplexF64)
                ap_density .= Array(dot(dws.st[b,r],psi[b,r]))
-                for t in 1:lp.iL[4]
+                @timeit "Volume sum" for t in 1:lp.iL[4]
                    for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
                        b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
                        Quark_cond2_cfl[t,noi,fl] += ap_density[b,r]
@@ -295,11 +295,11 @@ end
 """
 function two_pt_lagrange()
-Computes the two point functions containing Lagrange multiplier fields using front flow.
+Computes the two point functions containing Lagrange multiplier fields using front flow. The convention is such that the quark condentate is 'pptilde_t' + c_fl*'pphat_t'
 Stores the two point functions in the variables 'pphat_t' and 'pptilde_t'.
 """
-function two_pt_lagrange() # Will be 2pt lagrange mult
+function Two_pt_lagrange() # Will be 2pt lagrange mult
    @timeit "CPU to GPU" copyto!(U,U_CPU)
@@ -337,15 +337,15 @@ function two_pt_lagrange() # Will be 2pt lagrange mult
            pphat_t[:,noi,fl] .= zero(Float64)
            pp_density .= Array(sum(norm2.(psi)))
-            for t in 1:lp.iL[4]
+            @timeit "Volume sum" for t in 1:lp.iL[4]
                for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
                    b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
-                    pphat_t[t,noi,fl] += pp_density[b,r]
+                    pphat_t[t,noi,fl] += -pp_density[b,r]
                end end end
            end
            pptilde_t[:,noi,fl] .= zero(Float64)
            ap_density .= Array(dot.(dws.st,psi))
-            for t in 1:lp.iL[4]
+            @timeit "Volume sum" for t in 1:lp.iL[4]
                for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
                    b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
                    pptilde_t[t,noi,fl] += ap_density[b,r]
@@ -360,15 +360,15 @@ function two_pt_lagrange() # Will be 2pt lagrange mult
        end
        pphat_t[:,noi,end] .= zero(Float64)
        pp_density .= Array(sum(norm2.(psi)))
-        for t in 1:lp.iL[4]
+        @timeit "Volume sum" for t in 1:lp.iL[4]
            for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
                b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
-                pphat_t[t,noi,end] += pp_density[b,r]
+                pphat_t[t,noi,end] += -pp_density[b,r]
            end end end
        end
        pptilde_t[:,noi,end] .= zero(Complex{Float64})
        ap_density .= Array(dot.(dws.st,psi))
-        for t in 1:lp.iL[4]
+        @timeit "Volume sum" for t in 1:lp.iL[4]
            for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
                b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
                pptilde_t[t,noi,end] += ap_density[b,r]

--- a/src/utils.jl
+++ b/src/utils.jl
@@ -8,62 +8,98 @@ function read_ff(name::String)
    file = BDIO_open(name,"r","FerFlow correlators")
+    while BDIO_get_uinfo(file) != 14
+        BDIO_seek!(file)
+    end
+    ihdr = Vector{Int32}(undef,1)
+    BDIO_read(file, ihdr)
+    if ihdr[1] == 1730280201
+        println("Reading file ",name)
+    else
+        error("Wrong IHDR in file ", name)
+    end
    while BDIO_get_uinfo(file) != 1
        BDIO_seek!(file)
    end
    dims=Vector{Int32}(undef,1)
    iL=Vector{Int32}(undef,4)
    dpars=Vector{Float64}(undef,2) # m0,csw
    th=Vector{ComplexF64}(undef,4)
+    N_noiseff=Vector{Int32}(undef,1)
+    nsteps=Vector{Int32}(undef,1)
+    N_noisebf=Vector{Int32}(undef,1)
    nflow=Vector{Int32}(undef,1)
-    N_noise=Vector{Int32}(undef,1)
+    tzeroff=Vector{Int32}(undef,1)
+    espilon=Vector{Int32}(undef,1)
+    tsource_ff=Vector{Int32}(undef,1)
+    tsource_bf=Vector{Int32}(undef,1)
    BDIO_read(file, dims)
    BDIO_read(file, iL)
    BDIO_read(file, dpars)
    BDIO_read(file, th)
+    BDIO_read(file, N_noiseff)
+    BDIO_read(file, nsteps)
+    BDIO_read(file, N_noisebf)
    BDIO_read(file, nflow)
+    BDIO_read(file, tzeroff)
+    BDIO_read(file, epsilon)
+    BDIO_read(file, tsourceff)
+    BDIO_read(file, tsourcebf)
    flow_times=Vector{Float64}(undef,nflow[1])
    nflow[1] > 0 ? BDIO_read(file, flow_times) : nothing
-    BDIO_read(file, N_noise)
-    N_noise = N_noise[1]
+    println("Reading data from a ",iL, " lattice.")
+    println("Frontflow: ",N_noiseff[1], " sources at euclidean time ", tsourceff[1] ,". Flowtime starting at ",tzeroff[1], " with ",nsteps[1]," steps of size ", epsilon[1])
+    println("Backflow: ",N_noisebf[1]," sources at euclidean time ",tsourcebf[1],". Flow times are ", flow_times)
+    N_noiseff = N_noiseff[1]
+    N_noisebf = N_noisebf[1]
    nflow = nflow[1]
+    nsteps = nsteps[1]
    TvecR = Vector{Float64}(undef,iL[4]);
    TvecC = Vector{ComplexF64}(undef,iL[4]);
    CN = Vector{ComplexF64}(undef,1);
    RN = Vector{Float64}(undef,1);
    ET = Vector{ComplexF64}(undef,1+nflow);
-    pp_corr = Array{Float64}(undef,1,iL[4],N_noise);
+    pp_corr = Array{Float64}(undef,1,iL[4],N_noiseff);
-    ap_corr = Array{ComplexF64}(undef,1,iL[4],N_noise);
+    ap_corr = Array{ComplexF64}(undef,1,iL[4],N_noiseff);
-    pp_corr_t = Array{Float64}(undef,1,iL[4],N_noise,nflow);
+    pp_corr_t = Array{Float64}(undef,1,iL[4],N_noiseff,nsteps+1);
-    ap_corr_t = Array{ComplexF64}(undef,1,iL[4],N_noise,nflow);
+    ap_corr_t = Array{ComplexF64}(undef,1,iL[4],N_noiseff,nsteps+1);
-    pp_corr_tfl = Array{Float64}(undef,1,iL[4],N_noise,nflow);
+    pphat_t = Array{Float64}(undef,1,N_noiseff,nsteps+1);
-    ap_corr_tfl = Array{ComplexF64}(undef,1,iL[4],N_noise,nflow);
+    pptilde_t = Array{ComplexF64}(undef,1,N_noise,nsteps+1);
+    Eoft = Array{Complex{Float64}}(undef,1,2+nsteps,iL[4]);
+    pp_corr_tfl = Array{Float64}(undef,1,iL[4],N_noisebf,nflow);
+    ap_corr_tfl = Array{ComplexF64}(undef,1,iL[4],N_noisebf,nflow);
    Sigma = Array{ComplexF64}(undef,1,iL[4],N_noise,nflow);
    Sigma_cfl = Array{ComplexF64}(undef,1,iL[4],N_noise,nflow);
    Sigma2 = Array{ComplexF64}(undef,1,iL[4],N_noise,nflow);
    Sigma2_cfl = Array{ComplexF64}(undef,1,iL[4],N_noise,nflow);
    ChiDchi = Array{ComplexF64}(undef,1,iL[4],N_noise,nflow);
    ChiDchi_cfl = Array{ComplexF64}(undef,1,iL[4],N_noise,nflow);
-    Phat_t = Array{Float64}(undef,1,N_noise,nflow);
-    Eoft = Array{Complex{Float64}}(undef,1,1+nflow,iL[4]);
+    pp_corr_mc = Array{Float64}(undef,0,iL[4],N_noiseff);
+    ap_corr_mc = Array{ComplexF64}(undef,0,iL[4],N_noiseff);
-    pp_corr_mc = Array{Float64}(undef,0,iL[4],N_noise);
+    pphat_t_mc = Array{Float64}(undef,0,N_noiseff,nsteps+1);
-    ap_corr_mc = Array{ComplexF64}(undef,0,iL[4],N_noise);
+    pptilde_t_mc = Array{ComplexF64}(undef,0,N_noiseff,nsteps+1);
-    pp_corr_t_mc = Array{Float64}(undef,0,iL[4],N_noise,nflow);
+    Eoft_mc = Array{Complex{Float64}}(undef,0,2+nsteps,iL[4]);
-    ap_corr_t_mc = Array{ComplexF64}(undef,0,iL[4],N_noise,nflow);
+    pp_corr_t_mc = Array{Float64}(undef,0,iL[4],N_noiseff,nsteps+1);
-    pp_corr_tfl_mc = Array{Float64}(undef,0,iL[4],N_noise,nflow);
+    ap_corr_t_mc = Array{ComplexF64}(undef,0,iL[4],N_noiseff,nsteps+1);
-    ap_corr_tfl_mc = Array{ComplexF64}(undef,0,iL[4],N_noise,nflow);
-    Sigma_mc = Array{ComplexF64}(undef,0,iL[4],N_noise,nflow);
+    pp_corr_tfl_mc = Array{Float64}(undef,0,iL[4],N_noisebf,nflow);
-    Sigma_cfl_mc = Array{ComplexF64}(undef,0,iL[4],N_noise,nflow);
+    ap_corr_tfl_mc = Array{ComplexF64}(undef,0,iL[4],N_noisebf,nflow);
-    Sigma2_mc = Array{ComplexF64}(undef,0,iL[4],N_noise,nflow);
+    Sigma_mc = Array{ComplexF64}(undef,0,iL[4],N_noisebf,nflow);
-    Sigma2_cfl_mc = Array{ComplexF64}(undef,0,iL[4],N_noise,nflow);
+    Sigma_cfl_mc = Array{ComplexF64}(undef,0,iL[4],N_noisebf,nflow);
-    ChiDchi_mc = Array{ComplexF64}(undef,0,iL[4],N_noise,nflow);
+    Sigma2_mc = Array{ComplexF64}(undef,0,iL[4],N_noisebf,nflow);
-    ChiDchi_cfl_mc = Array{ComplexF64}(undef,0,iL[4],N_noise,nflow);
+    Sigma2_cfl_mc = Array{ComplexF64}(undef,0,iL[4],N_noisebf,nflow);
-    Phat_t_mc = Array{Float64}(undef,0,N_noise,nflow);
+    ChiDchi_mc = Array{ComplexF64}(undef,0,iL[4],N_noisebf,nflow);
-    Eoft_mc = Array{Complex{Float64}}(undef,0,1+nflow,iL[4]);
+    ChiDchi_cfl_mc = Array{ComplexF64}(undef,0,iL[4],N_noisebf,nflow);
    while BDIO_get_uinfo(file) != 8
        BDIO_seek!(file)
@@ -71,20 +107,33 @@ function read_ff(name::String)
    while BDIO_get_uinfo(file) == 8
-        for noi in 1:N_noise
+        for noi in 1:N_noiseff
            BDIO_read(file,TvecR)
            pp_corr[1,:,noi] .= TvecR
            BDIO_read(file,TvecC)
            ap_corr[1,:,noi] .= TvecC
-            for fl in 1:nflow
+            for fl in 1:nsteps+1
                BDIO_read(file,TvecR)
                pp_corr_t[1,:,noi,fl].= TvecR
                BDIO_read(file,TvecC)
                ap_corr_t[1,:,noi,fl].= TvecC
+                BDIO_read(file,TvecR)
+                pphat_t[1,:,noi,fl].= TvecR
+                BDIO_read(file,TvecC)
+                pptilde_t[1,:,noi,fl].= TvecC
+            end
+        end
+        for time in 1:iL[4]
+            BDIO_read(file,ET)
+            Eoft[1,:,time] .= ET
+        end
+        for noi in 1:N_noisebf
+            for fl in 1:nflow
                BDIO_read(file,TvecR)
                pp_corr_tfl[1,:,noi,fl].= TvecR
@@ -98,28 +147,27 @@ function read_ff(name::String)
                BDIO_read(file,TvecC)
                Sigma2[1,:,noi,fl] .= TvecC
                BDIO_read(file,TvecC)
                Sigma2_cfl[1,:,noi,fl] .= TvecC
                BDIO_read(file,TvecC)
                ChiDchi[1,:,noi,fl] .= TvecC
                BDIO_read(file,TvecC)
                ChiDchi_cfl[1,:,noi,fl] .= TvecC
-                BDIO_read(file,RN)
-                Phat_t[1,noi,fl]= RN[1]
            end
        end
-        for time in 1:iL[4]
-            BDIO_read(file,ET)
-            Eoft[1,:,time] .= ET
-        end
        Eoft_mc=[Eoft_mc;Eoft]
        pp_corr_mc=[pp_corr_mc;pp_corr]
        ap_corr_mc=[ap_corr_mc;ap_corr]
        pp_corr_t_mc=[pp_corr_t_mc;pp_corr_t]
        ap_corr_t_mc=[ap_corr_t_mc;ap_corr_t]
+        pphat_t_mc=[pphat_t_mc;pphat_t]
+        pptilde_t_mc=[pptilde_t_mc;pptilde_t]
        pp_corr_tfl_mc=[pp_corr_tfl_mc;pp_corr_tfl]
        ap_corr_tfl_mc=[ap_corr_tfl_mc;ap_corr_tfl]
        Sigma_mc=[Sigma_mc;Sigma]
@@ -128,7 +176,6 @@ function read_ff(name::String)
        Sigma2_cfl_mc=[Sigma2_cfl_mc;Sigma2_cfl]
        ChiDchi_mc=[ChiDchi_mc;ChiDchi]
        ChiDchi_cfl_mc=[ChiDchi_cfl_mc;ChiDchi_cfl]
-        Phat_t_mc=[Phat_t_mc;Phat_t]
        BDIO_seek!(file)
        BDIO_seek!(file)
@@ -136,44 +183,47 @@ function read_ff(name::String)
    BDIO_close!(file)
-    return Eoft_mc,pp_corr_mc,ap_corr_mc,pp_corr_t_mc,ap_corr_t_mc,pp_corr_tfl_mc,ap_corr_tfl_mc,Sigma_mc,Sigma_cfl_mc,Sigma2_mc,Sigma2_cfl_mc,ChiDchi_mc,ChiDchi_cfl_mc,Phat_t_mc
+    return Eoft_mc,pp_corr_mc,ap_corr_mc,pp_corr_t_mc,ap_corr_t_mc,pp_corr_tfl_mc,ap_corr_tfl_mc,Sigma_mc,Sigma_cfl_mc,Sigma2_mc,Sigma2_cfl_mc,ChiDchi_mc,ChiDchi_cfl_mc,pphat_t_mc,pptilde_t_mc
 end
 """
 function uwff(file::String)
-Loads in the variables 'Eoft', 'pp_corr', 'ap_corr', 'pp_corr_t','ap_corr_t', 'Sigma' and 'Sigma_cfl' the corresponding
+Loads in the variables 'Eoft', 'pp_corr', 'ap_corr', 'pp_corr_t','ap_corr_t', 'pphat_t', 'pptilde_t','pp_corr_tfl','ap_corr_tfl' 'Sigma'(v1 and v2), 'Sigma_cfl'(v1 and v2), 'ChiDchi' and 'ChiDchi_cfl' the corresponding
-'uwreal' quantities. The indices are, when corresponding, eucliden time and flow times
+'uwreal' quantities. Noise average has been performed. The indices are, when corresponding, eucliden time and flow times
 """
 function uwff(file::String)
-    Eoft_mc,pp_corr_mc,ap_corr_mc,pp_corr_t_mc,ap_corr_t_mc,pp_corr_tfl_mc,ap_corr_tfl_mc,Sigma_mc,Sigma_cfl_mc,Sigma2_mc,Sigma2_cfl_mc,ChiDchi_mc,ChiDchi_cfl_mc,Phat_t_mc = read_ff(file)
+    Eoft_mc,pp_corr_mc,ap_corr_mc,pp_corr_t_mc,ap_corr_t_mc,pp_corr_tfl_mc,ap_corr_tfl_mc,Sigma_mc,Sigma_cfl_mc,Sigma2_mc,Sigma2_cfl_mc,ChiDchi_mc,ChiDchi_cfl_mc,pphat_t_mc,pptilde_t_mc = read_ff(file)
    runame = String(split(split(file,"/")[end],".")[1])
    T = size(pp_corr_mc)[2]
-    Ns = size(pp_corr_mc)[3]
+    Nsff = size(pp_corr_mc)[3]
-    Nfl = size(pp_corr_t_mc)[4]
+    nsteps = size(pp_corr_t_mc)[4]
+    Nsbf = size(pp_corr_tfl_mc)[3]
+    Nfl = size(pp_corr_tfl_mc)[4]
-    global Eoft = [uwreal(real.(Eoft_mc[:,i,j]),runame) for j in 1:T, i in 1:Nfl+1]
+    global Eoft = [uwreal(real.(Eoft_mc[:,i,j]),runame) for j in 1:T, i in 1:nsteps+2]
-    global pp_corr = [uwreal(sum(pp_corr_mc,dims = 3)[:,i,1]./Ns,runame) for i in 1:T]
+    global pp_corr = [uwreal(sum(pp_corr_mc,dims = 3)[:,i,1]./Nsff,runame) for i in 1:T]
-    global ap_corr = [uwreal(real.(sum(ap_corr_mc,dims = 3)[:,i,1]./Ns),runame) for i in 1:T]
+    global ap_corr = [uwreal(real.(sum(ap_corr_mc,dims = 3)[:,i,1]./Nsff),runame) for i in 1:T]
-    global pp_corr_t = [uwreal(sum(pp_corr_t_mc,dims = 3)[:,i,1,k]./Ns,runame) for i in 1:T, k in 1:Nfl]
+    global pp_corr_t = [uwreal(sum(pp_corr_t_mc,dims = 3)[:,i,1,k]./Nsff,runame) for i in 1:T, k in 1:nsteps+1]
-    global ap_corr_t = [uwreal(real.(sum(ap_corr_t_mc,dims = 3)[:,i,1,k]./Ns),runame) for i in 1:T,k in 1:Nfl]
+    global ap_corr_t = [uwreal(real.(sum(ap_corr_t_mc,dims = 3)[:,i,1,k]./Nsff),runame) for i in 1:T,k in 1:nsteps+1]
-    global pp_corr_tfl = [uwreal(sum(pp_corr_tfl_mc,dims = 3)[:,i,1,k]./Ns,runame) for i in 1:T, k in 1:Nfl]
+    global pphat_t = [uwreal(sum(pphat_t_mc,dims = 3)[:,i,1,k]./Nsff,runame) for i in 1:T, k in 1:nsteps+1]
-    global ap_corr_tfl = [uwreal(real.(sum(ap_corr_tfl_mc,dims = 3)[:,i,1,k]./Ns),runame) for i in 1:T,k in 1:Nfl]
+    global pptilde_corr_t = [uwreal(real.(sum(pptilde_t_mc,dims = 3)[:,i,1,k]./Nsff),runame) for i in 1:T,k in 1:nsteps+1]
-    global Sigma = [uwreal(real.(sum(Sigma_mc,dims = 3)[:,i,1,k]./Ns),runame) for i in 1:T, k in 1:Nfl]
+    global pp_corr_tfl = [uwreal(sum(pp_corr_tfl_mc,dims = 3)[:,i,1,k]./Nsbf,runame) for i in 1:T, k in 1:Nfl]
-    global Sigma_cfl = [uwreal(real.(sum(Sigma_cfl_mc,dims = 3)[:,i,1,k]./Ns),runame) for i in 1:T, k in 1:Nfl]
+    global ap_corr_tfl = [uwreal(real.(sum(ap_corr_tfl_mc,dims = 3)[:,i,1,k]./Nsbf),runame) for i in 1:T,k in 1:Nfl]
-    global Sigma2 = [uwreal(real.(sum(Sigma2_mc,dims = 3)[:,i,1,k]./Ns),runame) for i in 1:T, k in 1:Nfl]
+    global Sigma = [uwreal(real.(sum(Sigma_mc,dims = 3)[:,i,1,k]./Nsbf),runame) for i in 1:T, k in 1:Nfl]
-    global Sigma2_cfl = [uwreal(real.(sum(Sigma2_cfl_mc,dims = 3)[:,i,1,k]./Ns),runame) for i in 1:T, k in 1:Nfl]
+    global Sigma_cfl = [uwreal(real.(sum(Sigma_cfl_mc,dims = 3)[:,i,1,k]./Nsbf),runame) for i in 1:T, k in 1:Nfl]
-    global ChiDchi = [uwreal(real.(sum(ChiDchi_mc,dims = 3)[:,i,1,k]./Ns),runame) for i in 1:T, k in 1:Nfl]
+    global Sigma2 = [uwreal(real.(sum(Sigma2_mc,dims = 3)[:,i,1,k]./Nsbf),runame) for i in 1:T, k in 1:Nfl]
-    global ChiDchi_cfl = [uwreal(real.(sum(ChiDchi_cfl_mc,dims = 3)[:,i,1,k]./Ns),runame) for i in 1:T, k in 1:Nfl]
+    global Sigma2_cfl = [uwreal(real.(sum(Sigma2_cfl_mc,dims = 3)[:,i,1,k]./Nsbf),runame) for i in 1:T, k in 1:Nfl]
-    global Phat_t = [uwreal(sum(Phat_t_mc,dims = 2)[:,1,i]./Ns,runame) for i in 1:Nfl]
+    global ChiDchi = [uwreal(real.(sum(ChiDchi_mc,dims = 3)[:,i,1,k]./Nsbf),runame) for i in 1:T, k in 1:Nfl]
+    global ChiDchi_cfl = [uwreal(real.(sum(ChiDchi_cfl_mc,dims = 3)[:,i,1,k]./Nsbf),runame) for i in 1:T, k in 1:Nfl]
    return nothing
 end