diff --git a/docs/make.jl b/docs/make.jl
index 45403ed2..b665c973 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -10,6 +10,7 @@ makedocs(;
     pages = [
         "Home" => "index.md"
         "Library" => "lib/lib.md"
+        "CFT Data" => "cft.md"
         "Finalizers" => "finalizers.md"
         "References" => "references.md"
     ],
diff --git a/docs/src/cft.md b/docs/src/cft.md
new file mode 100644
index 00000000..c78e460f
--- /dev/null
+++ b/docs/src/cft.md
@@ -0,0 +1,32 @@
+# Conformal Field Theory Data
+TNRKit provides extensive tools for calculating conformal field theory data. Details about the implementation can be found in the TNRKit paper ([arxiv/2604.06922](https://arxiv.org/abs/2604.06922)).
+
+The core idea behind calculating the central charge, scaling dimensions, and conformal spins, is to calculate the spectrum of the fixed point tensor on a tube geometry. There are different ways to put fixed point tensors on a tube and the geometry of this tube is characterised by 3 parameters:
+$$[h, L, x]$$
+Where $h$ is the height of the tube, $L$ is the circumference, and $x$ is the horizontal shift. The higher the ratio $\frac{L}{h}$, the higher the resolution but also the more expensive the calculation.
+
+To calculate cft data we provide the `CFTData` struct which can be used in the following ways:
+
+```julia
+CFTData(scheme; shape=[h, L, x])
+CFTData(T::TensorMap; kwargs...) # 1 fixed point tensor
+CFTData(TA::TensorMap, TB::TensorMap; kwargs...) # 2x2 checkerboard unitcell
+```
+
+The shapes we provide are: $[1, 1, 0]$, $[\sqrt{2}, 2\sqrt{2}, 0]$, $[1, 4, 1]$, $[1, 8, 1]$, $[\frac{4}{\sqrt{10}}, 2 \sqrt{10}, \frac{2}{\sqrt{10}}]$
+
+The last two of which require intermediate truncation steps, the parameters of which can be tuned by:
+```julia
+CFTData(scheme; shape=[1, 8, 1], trunc = trunc1, truncentanglement=trunc2)
+```
+
+# CFTData struct
+The `CFTData` struct has two fields:
+- `central_charge`
+- `scaling_dimensions`
+
+The `central_charge` can be either `missing` (when using the $[1, 1, 0]$ shape), or a number.
+The `scaling_dimensions` field is a `SectorVector` from TensorKit.jl.
+
+The `scaling_dimensions` can be indexed like an `AbstractVector` (i.e. with scalars, slices, ...), or with sectors (e.g. `Z2Irrep(0)`), which will provide the scaling dimensions associated with that sector/charge.
+To check which sectors you can index the `scaling_dimensions` with you can use `keys(scaling_dimensions)`.
\ No newline at end of file
diff --git a/src/TNRKit.jl b/src/TNRKit.jl
index 4ed21251..1779ce99 100644
--- a/src/TNRKit.jl
+++ b/src/TNRKit.jl
@@ -120,14 +120,14 @@ include("utility/free_energy.jl")
 export free_energy
 
 include("utility/cft.jl")
-export cft_data, central_charge
+export CFTData, central_charge
 
 include("utility/gs_degeneracy.jl")
 export ground_state_degeneracy, gu_wen_ratio
 
 include("utility/finalize.jl")
 export Finalizer, two_by_two_Finalizer, finalize!, finalize_two_by_two!, finalize_cftdata!, finalize_central_charge!,
-    finalize_groundstatedegeneracy!, GSDegeneracy_Finalizer, guwenratio_Finalizer
+    finalize_groundstatedegeneracy!, CFT_Finalizer, GSDegeneracy_Finalizer, guwenratio_Finalizer
 
 include("utility/cdl.jl")
 export cdl_tensor
diff --git a/src/schemes/ctm/ctm_hotrg.jl b/src/schemes/ctm/ctm_hotrg.jl
index 170214b4..8db97479 100644
--- a/src/schemes/ctm/ctm_hotrg.jl
+++ b/src/schemes/ctm/ctm_hotrg.jl
@@ -7,7 +7,7 @@ Corner Transfer Matrix environment + Higher-Order Tensor Renormalization Group
     $(FUNCTIONNAME)(T, χenv[, ctm_iter=maxiter(2.0e4), ctm_tol=trivial_convcrit(1.0e-9), ctm_obc=false, χenv_ini=2])
 
 ### Running the algorithm
-    run!(::ctm_HOTRG, trunc::TruncationStrategy, criterion::stopcrit[, sweep=30, return_cft=false, inv=false, conv_criterion=1.0e-12])
+    run!(::ctm_HOTRG, trunc::TruncationStrategy, criterion::stopcrit[, sweep=30, inv=false, conv_criterion=1.0e-12])
 
 ### Fields
 
@@ -129,29 +129,21 @@ end
 
 function run!(
         scheme::ctm_HOTRG, trunc::TruncationStrategy, criterion::stopcrit;
-        sweep = 30, return_cft = false, inv = false, conv_criterion = 1.0e-12
+        sweep = 30, inv = false, conv_criterion = 1.0e-12
     )
     area = 1
     lnz = 0.0
-    cft = []
 
     for i in 1:(criterion.n)
         area *= 4.0
         tr_norm = step!(scheme, trunc; sweep = sweep, inv = inv)
-        if return_cft
-            push!(cft, cft_data(scheme; unitcell = 2))
-        end
         lnz += log(tr_norm) / area
         if abs(log(abs(tr_norm)) / area) <= conv_criterion
             @info "CTM-HOTRG converged after $i iterations!"
             break
         end
     end
-    if return_cft
-        return lnz, cft
-    else
-        return lnz
-    end
+    return lnz
 end
 
 function Base.show(io::IO, scheme::ctm_HOTRG)
diff --git a/src/schemes/ctm/ctm_trg.jl b/src/schemes/ctm/ctm_trg.jl
index 9070a5ab..be3ee929 100644
--- a/src/schemes/ctm/ctm_trg.jl
+++ b/src/schemes/ctm/ctm_trg.jl
@@ -7,7 +7,7 @@ Corner Transfer Matrix environment + Tensor Renormalization Group
     $(FUNCTIONNAME)(T, χenv[, ctm_iter=2.0e4, ctm_tol=1.0e-9])
 
 ### Running the algorithm
-    run!(::ctm_TRG, trunc::TruncationStrategy, criterion::maxiter[, sweep=30, enlarge=true, return_cft=false, inv=false, conv_criterion=1.0e-12, modified=true])
+    run!(::ctm_TRG, trunc::TruncationStrategy, criterion::maxiter[, sweep=30, enlarge=true, inv=false, conv_criterion=1.0e-12, modified=true])
 
 ### Fields
 
@@ -167,14 +167,12 @@ function run!(
         criterion::maxiter;
         sweep = 30,
         enlarge = true,
-        return_cft = false,
         inv = false,
         conv_criterion = 1.0e-12,
         modified = true,
     )
     area = 1
     lnz = 0.0
-    cft = []
 
     steps = 0
     crit = true
@@ -182,9 +180,6 @@ function run!(
         area *= 4.0
         tr_norm = step!(scheme, trunc; sweep = sweep, enlarge = enlarge, inv = inv, modified)
         lnz += log(tr_norm) / area
-        if return_cft
-            push!(cft, cft_data(scheme; unitcell = 2))
-        end
         if abs(log(abs(tr_norm)) / area) <= conv_criterion
             @info "CTM-TRG converged after $steps iterations"
             break
@@ -192,11 +187,7 @@ function run!(
         steps += 1
         crit = criterion(steps, nothing)
     end
-    if return_cft
-        return lnz, cft
-    else
-        return lnz
-    end
+    return lnz
 end
 
 function Base.show(io::IO, scheme::ctm_TRG)
diff --git a/src/utility/cft.jl b/src/utility/cft.jl
index 07af97b4..4dd5afbd 100644
--- a/src/utility/cft.jl
+++ b/src/utility/cft.jl
@@ -1,75 +1,128 @@
-function next_τ(τ)
-    return (τ - 1) / (τ + 1)
+# Extensions on top of TensorKit.SectorVector
+Base.broadcasted(f, v::TensorKit.SectorVector) = TensorKit.SectorVector(broadcast(f, parent(v)), v.structure)
+Base.broadcasted(f, v::TensorKit.SectorVector, a) = TensorKit.SectorVector(broadcast(f, parent(v), a), v.structure)
+Base.broadcasted(f, a, v::TensorKit.SectorVector) = TensorKit.SectorVector(broadcast(f, a, parent(v)), v.structure)
+function Base.broadcasted(f, v1::TensorKit.SectorVector, v2::TensorKit.SectorVector)
+    if v1.structure != v2.structure
+        throw(ArgumentError("Cannot broadcast two SectorVectors with different structures"))
+    end
+    return TensorKit.SectorVector(broadcast(f, parent(v1), parent(v2)), v1.structure)
 end
 
-function cft_data(scheme::TNRScheme; v = 1, unitcell = 1, is_real = true)
-    # make the indices
-    indices = [[i, -i, -(i + unitcell), i + 1] for i in 1:unitcell]
-    indices[end][4] = 1
+function Base.filter(f, v::TensorKit.SectorVector)
+    data = copy(parent(v))
+    structure = copy(v.structure)
 
-    T = ncon(fill(scheme.T, unitcell), indices)
+    kept_inds = findall(f, parent(v))
+    sectors = keys(structure)
+    for (i, sector) in enumerate(sectors)
+        structure[sector] = i == 1 ? (1:findlast(x -> x <= structure[sector].stop, kept_inds)) : ((structure[sectors[i - 1]].stop + 1):findlast(x -> x <= structure[sector].stop, kept_inds))
+    end
+    data = data[kept_inds]
+    return TensorKit.SectorVector(data, structure)
+end
 
-    outinds = Tuple(collect(1:unitcell))
-    ininds = Tuple(collect((unitcell + 1):(2unitcell)))
+function Base.sort(v::TensorKit.SectorVector; kwargs...)
+    # sort within the sectors, then concatenate the data and update the structure
+    # Ideally this would sort the total data, but sectorvectors are only compatible with
+    # structures that contain unitranges, we cannot interweave the data of different sectors.
+    data = copy(parent(v))
+    newdata = similar(data)
+    structure = copy(v.structure)
+    sectors = keys(structure)
+    for sector in sectors
+        newdata[structure[sector]] = sort(data[structure[sector]]; kwargs...)
+    end
+    return TensorKit.SectorVector(newdata, structure)
+end
 
-    T = permute(T, (outinds, ininds))
-    D, _ = eig_full(T)
+Base.:*(a::Number, v::TensorKit.SectorVector) = scale(v, a)
+Base.:*(v::TensorKit.SectorVector, a::Number) = scale(v, a)
 
-    data = zeros(ComplexF64, dim(space(D, 1)))
+"""
+    CFTData{E, I} where {E, I}
 
-    i = 1
-    for (_, b) in blocks(D)
-        for I in LinearAlgebra.diagind(b)
-            data[i] = b[I]
-            i += 1
-        end
-    end
+A struct to hold conformal data extracted from a TNR scheme.
 
-    data = sort(data; by = x -> abs(x), rev = true) # sorting by magnitude
-    data = filter(x -> real(x) > 0, data) # filtering out negative real values
-    data = filter(x -> abs(x) > 1.0e-12, data) # filtering out small values
+### Constructors
 
-    if is_real
-        data = real(data)
+
+### Fields
+    - `central_charge::Union{E, Missing}`: The central charge of the CFT. Will be `nothing` if not calculated.
+    - `scaling_dimensions::TensorKit.SectorVector{E, I}`: The scaling dimensions of the CFT, organized in a `TensorKit.SectorVector` where the sectors correspond to different spin sectors (or other quantum numbers) and the data contains the scaling dimensions within those sectors
+
+"""
+struct CFTData{E, I}
+    "Central charge of the CFT. Will be `nothing` if not calculated."
+    central_charge::Union{E, Missing}
+
+    "Scaling dimensions of the CFT."
+    scaling_dimensions::TensorKit.SectorVector{E, I}
+end
+
+function Base.show(io::IO, data::CFTData)
+    println(io, "CFTData")
+    println(io, "  * central charge: $(data.central_charge)")
+    println(io, "  * scaling dimensions: $(data.scaling_dimensions)")
+    return nothing
+end
+
+function CFTData(T::TensorMap{E, S, 2, 2}; shape = [sqrt(2), 2 * sqrt(2), 0], kwargs...) where {E, S}
+    if shape == [1, 1, 0] # trivial implementation
+        return CFTData(missing, _scaling_dimensions(T))
+    else
+        CFTData(T, T; shape, kwargs...)
     end
+end
+CFTData(scheme::TNRScheme; kwargs...) = CFTData(scheme.T; kwargs...) # simple 1-site unitcell schemes
+CFTData(scheme::LoopTNR; kwargs...) = CFTData(scheme.TA, scheme.TB; kwargs...) # simple 1-site unitcell schemes
+function CFTData(scheme::BTRG; kwargs...) # merge bond tensors into central tensor
+    @tensor T_unit[-1 -2; -3 -4] := scheme.T[1 2; -3 -4] * scheme.S1[-2; 2] *
+        scheme.S2[-1; 1]
+    return CFTData(T_unit; kwargs...)
+end
 
-    return unitcell * (1 / (2π * v)) * log.(data[1] ./ data)
+# Main implementation, two-site unitcell
+function CFTData(TA::TensorMap{E, S, 2, 2}, TB::TensorMap{E, S, 2, 2}; shape = [sqrt(2), 2 * sqrt(2), 0], trunc = truncrank(16), truncentanglement = trunctol(; rtol = 1.0e-14)) where {E, S}
+    if shape == [1, 1, 0]
+        throw(ArgumentError("The shape [1, 1, 0] is not compatible with a two-site unit cell."))
+    elseif (shape ≈ [sqrt(2), 2 * sqrt(2), 0]) || (shape == [1, 4, 1]) # these shapes need no truncation
+        norm_const = area_term(TA, TB)^(1 / 4) # canonical normalisation constant
+        return spec(TA / norm_const, TB / norm_const, shape)
+    elseif (shape == [1, 8, 1]) || (shape ≈ [4 / sqrt(10), 2 * sqrt(10), 2 / sqrt(10)])
+
+        norm_const = area_term(TA, TB)^(1 / 4) # canonical normalisation constant
+
+        dl, ur, ul, dr = MPO_opt(
+            TA / norm_const, TB / norm_const, trunc, truncentanglement
+        )
+        T = reduced_MPO(dl, ur, ul, dr, trunc)
+        return spec(T, T, shape)
+    else
+        throw(ArgumentError("Shape $shape is not implemented."))
+    end
 end
 
-function cft_data(scheme::BTRG; v = 1, unitcell = 1, is_real = true)
-    # make the indices
+# Trivial diagonalisation of the transfer matrix. Currently the v and unitcell are not acessible from the outside.
+# The user should really be using the other shapes anyways.
+function _scaling_dimensions(T::TensorMap{E, S, 2, 2}; v = 1, unitcell = 1) where {E, S}
+    # stack unitcell copies of T and trace
     indices = [[i, -i, -(i + unitcell), i + 1] for i in 1:unitcell]
     indices[end][4] = 1
 
-    @tensor T_unit[-1 -2; -3 -4] := scheme.T[1 2; -3 -4] * scheme.S1[-2; 2] *
-        scheme.S2[-1; 1]
-    T = ncon(fill(T_unit, unitcell), indices)
+    T = ncon(fill(T, unitcell), indices)
 
     outinds = Tuple(collect(1:unitcell))
     ininds = Tuple(collect((unitcell + 1):(2unitcell)))
 
     T = permute(T, (outinds, ininds))
-    D, _ = eig_full(T)
-
-    data = zeros(ComplexF64, dim(space(D, 1)))
-
-    i = 1
-    for (_, b) in blocks(D)
-        for I in LinearAlgebra.diagind(b)
-            data[i] = b[I]
-            i += 1
-        end
-    end
 
+    data = eig_vals(T)
     data = sort(data; by = x -> abs(x), rev = true) # sorting by magnitude
     data = filter(x -> real(x) > 0, data) # filtering out negative real values
     data = filter(x -> abs(x) > 1.0e-12, data) # filtering out small values
 
-    if is_real
-        data = real(data)
-    end
-
-    return unitcell * (1 / (2π * v)) * log.(data[1] ./ data)
+    return unitcell * (1 / (2π * v)) .* log.(data[1] ./ data)
 end
 
 """
@@ -98,64 +151,7 @@ function area_term(A, B; is_real = true)
     end
 end
 
-function MPO_opt(
-        TA::TensorMap, TB::TensorMap, trunc::TruncationStrategy,
-        truncentanglement::TruncationStrategy
-    )
-    pretrunc = truncrank(2 * trunc.howmany)
-    dl, ur = SVD12(TA, pretrunc)
-    dr, ul = SVD12(transpose(TB, ((2, 4), (1, 3))), pretrunc)
-
-    transfer_MPO = [
-        transpose(dl, ((1,), (3, 2))), ur, transpose(ul, ((2,), (3, 1))),
-        transpose(dr, ((3,), (2, 1))),
-    ]
-
-    in_inds = [1, 1, 1, 1]
-    out_inds = [1, 2, 2, 1]
-    MPO_function(steps, data) = abs(data[end])
-    criterion = maxiter(10) & convcrit(1.0e-12, MPO_function)
-    PR_list, PL_list = find_projectors(
-        transfer_MPO, in_inds, out_inds, criterion,
-        trunc & truncentanglement
-    )
-
-    MPO_disentangled!(transfer_MPO, in_inds, out_inds, PR_list, PL_list)
-    return transfer_MPO
-end
-
-function reduced_MPO(
-        dl::TensorMap, ur::TensorMap, ul::TensorMap, dr::TensorMap,
-        trunc::TruncationStrategy
-    )
-    @plansor temp[-1 -2; -3 -4] := ur[-1; 1 4] *
-        ul[4; 3 -2] *
-        dr[-3; 2 1] * dl[2; -4 3]
-    D, U = SVD12(temp, trunc)
-    @plansor translate[-1 -2; -3 -4] := U[-2; 1 -4] * D[-1 1; -3]
-    return translate
-end
-
-function MPO_action_1x4(TA::TensorMap, TB::TensorMap, x::TensorMap)
-    @tensor TTTTx[-1 -2 -3 -4; -5] := x[1 2 3 4; -5] * TA[41 -1; 1 12] *
-        TB[12 -2; 2 23] *
-        TA[23 -3; 3 34] * TB[34 -4; 4 41]
-    return TTTTx
-end
-
-function MPO_action_1x4_twist(TA::TensorMap, TB::TensorMap, x::TensorMap)
-    TTTTx = MPO_action_1x4(TA, TB, x)
-    return permute(TTTTx, ((2, 3, 4, 1), (5,)))
-end
-
-# Fig.25 of https://arxiv.org/pdf/2311.18785. Firstly appear in Chenfeng Bao's thesis, see http://hdl.handle.net/10012/14674.
-function MPO_action_2gates(TA::TensorMap, TB::TensorMap, x::TensorMap)
-    @tensor fx[-1 -2 -3 -4; 5] := TB[-1 -2; 1 2] * x[1 2 3 4; 5] * TB[-3 -4; 3 4]
-    @tensor ffx[-1 -2 -3 -4; 5] := TA[-3 -4; 2 3] * fx[1 2 3 4; 5] *
-        TA[-1 -2; 4 1]
-    return permute(ffx, ((2, 3, 4, 1), (5,)))
-end
-
+# The case with spin is based on https://arxiv.org/pdf/1512.03846 and some private communications with Yingjie Wei and Atsushi Ueda
 function spec(TA::TensorMap, TB::TensorMap, shape::Array; Nh = 25)
     area = shape[1] * shape[2]
     Imτ = shape[1] / shape[2]
@@ -198,72 +194,88 @@ function spec(TA::TensorMap, TB::TensorMap, shape::Array; Nh = 25)
         end
     )
 
-    conformal_data = Dict()
-
     norm_const_0 = spec_sector[one(I)][1]
-    conformal_data["c"] = 6 / pi / (Imτ - area / 4) * log(norm_const_0)
+    central_charge = 6 / pi / (Imτ - area / 4) * log(norm_const_0)
 
+    # Construct a SectorVector from the data of the different sectors
+    data = ComplexF64[]
+    structure = TensorKit.SectorDict{sectortype(xspace), UnitRange{Int}}()
+    last_index = 1
     for charge in sectors(fuse(xspace))
         DeltaS = -1 / (2 * pi * Imτ) * log.(spec_sector[charge] / norm_const_0)
         if !(relative_shift ≈ 0)
-            conformal_data[charge] = real.(DeltaS) + imag.(DeltaS) / relative_shift * im
+            push!(data, (real.(DeltaS) + imag.(DeltaS) / relative_shift * im)...)
+            structure[charge] = last_index:(last_index + length(DeltaS) - 1)
         else
-            conformal_data[charge] = DeltaS
+            push!(data, real.(DeltaS)...)
+            structure[charge] = last_index:(last_index + length(DeltaS) - 1)
         end
+        last_index += length(DeltaS)
     end
-    return conformal_data
+
+    sv = TensorKit.SectorVector(data, structure)
+    sv = sort(sv; by = real)
+    sv = filter(x -> real(x) ≤ 1.0e16, sv)
+
+    return CFTData(central_charge, sv)
 end
 
-# The function to obtain central charge and conformal spectrum from the fixed-point tensor with G-symmetry. Here the conformal spectrum is obtained by different charge sectors.
-# The case with spin is based on https://arxiv.org/pdf/1512.03846 and some private communications with Yingjie Wei and Atsushi Ueda
-function cft_data(
-        scheme::LoopTNR, shape::Array,
-        trunc::TruncationStrategy,
+function MPO_opt(
+        TA::TensorMap, TB::TensorMap, trunc::TruncationStrategy,
         truncentanglement::TruncationStrategy
     )
-    if !(shape in [[1, 8, 1], [4 / sqrt(10), 2 * sqrt(10), 2 / sqrt(10)]])
-        throw(ArgumentError("The shape $shape is not correct."))
-    end
+    pretrunc = truncrank(2 * trunc.howmany)
+    dl, ur = SVD12(TA, pretrunc)
+    dr, ul = SVD12(transpose(TB, ((2, 4), (1, 3))), pretrunc)
+
+    transfer_MPO = [
+        transpose(dl, ((1,), (3, 2))), ur, transpose(ul, ((2,), (3, 1))),
+        transpose(dr, ((3,), (2, 1))),
+    ]
 
-    @infov 2 "CFT data calculating"
-    norm_const = area_term(scheme.TA, scheme.TB)^(1 / 4)
-    dl, ur, ul, dr = MPO_opt(
-        scheme.TA / norm_const, scheme.TB / norm_const, trunc, truncentanglement
+    in_inds = [1, 1, 1, 1]
+    out_inds = [1, 2, 2, 1]
+    MPO_function(steps, data) = abs(data[end])
+    criterion = maxiter(10) & convcrit(1.0e-12, MPO_function)
+    PR_list, PL_list = find_projectors(
+        transfer_MPO, in_inds, out_inds, criterion,
+        trunc & truncentanglement
     )
-    T = reduced_MPO(dl, ur, ul, dr, trunc)
 
-    # Calculate conformal data with spin from -4 to 4. Most error is introduced in the second step of the SVD.
-    conformal_data = spec(T, T, shape)
-    return conformal_data
+    MPO_disentangled!(transfer_MPO, in_inds, out_inds, PR_list, PL_list)
+    return transfer_MPO
 end
 
-function cft_data(scheme::LoopTNR, shape::Array)
-    if !(shape in [[1, 4, 1], [sqrt(2), 2 * sqrt(2), 0]])
-        throw(ArgumentError("The shape $shape is not correct."))
-    end
-
-    @infov 2 "CFT data calculating"
-    norm_const = area_term(scheme.TA, scheme.TB)^(1 / 4)
-    conformal_data = spec(scheme.TA / norm_const, scheme.TB / norm_const, shape)
-    return conformal_data
+# Apply functions for diagonalising different shapes of transfer matrices
+# =======================================================================
+# Fig.25 of https://arxiv.org/pdf/2311.18785. Firstly appear in Chenfeng Bao's thesis, see http://hdl.handle.net/10012/14674.
+function MPO_action_2gates(TA::TensorMap, TB::TensorMap, x::TensorMap)
+    @tensor fx[-1 -2 -3 -4; 5] := TB[-1 -2; 1 2] * x[1 2 3 4; 5] * TB[-3 -4; 3 4]
+    @tensor ffx[-1 -2 -3 -4; 5] := TA[-3 -4; 2 3] * fx[1 2 3 4; 5] *
+        TA[-1 -2; 4 1]
+    return permute(ffx, ((2, 3, 4, 1), (5,)))
 end
 
-"""
-    central_charge(scheme::TNRScheme, n::Number)
+function MPO_action_1x4(TA::TensorMap, TB::TensorMap, x::TensorMap)
+    @tensor TTTTx[-1 -2 -3 -4; -5] := x[1 2 3 4; -5] * TA[41 -1; 1 12] *
+        TB[12 -2; 2 23] *
+        TA[23 -3; 3 34] * TB[34 -4; 4 41]
+    return TTTTx
+end
 
-Get the central charge given the current state of a `TNRScheme` and the previous normalization factor `n`
-"""
-function central_charge(scheme::TNRScheme, n::Number)
-    @tensor M[-1; -2] := (scheme.T / n)[1 -1; -2 1]
-    _, S, _ = svd_compact(M)
-    return log(S.data[1]) * 6 / (π)
+function MPO_action_1x4_twist(TA::TensorMap, TB::TensorMap, x::TensorMap)
+    TTTTx = MPO_action_1x4(TA, TB, x)
+    return permute(TTTTx, ((2, 3, 4, 1), (5,)))
 end
 
-function central_charge(scheme::BTRG, n::Number)
-    @tensor M[-1; -2] := (
-        (scheme.T)[1 -1; 3 2] * scheme.S1[3; -2] *
-            scheme.S2[2; 1]
-    ) / n
-    _, S, _ = svd_compact(M)
-    return log(S.data[1]) * 6 / (π)
+function reduced_MPO(
+        dl::TensorMap, ur::TensorMap, ul::TensorMap, dr::TensorMap,
+        trunc::TruncationStrategy
+    )
+    @plansor temp[-1 -2; -3 -4] := ur[-1; 1 4] *
+        ul[4; 3 -2] *
+        dr[-3; 2 1] * dl[2; -4 3]
+    D, U = SVD12(temp, trunc)
+    @plansor translate[-1 -2; -3 -4] := U[-2; 1 -4] * D[-1 1; -3]
+    return translate
 end
diff --git a/src/utility/finalize.jl b/src/utility/finalize.jl
index ec2dff2c..1f9901fb 100644
--- a/src/utility/finalize.jl
+++ b/src/utility/finalize.jl
@@ -124,12 +124,14 @@ function finalize_phase23!(scheme::CorrelationHOTRG)
 end
 
 # cft data finalize
-function finalize_cftdata!(scheme::LoopTNR)
+function finalize_cftdata!(scheme::TNRScheme)
     finalize!(scheme)
-    return cft_data(scheme; is_real = true)
+    return CFTData(scheme)
 end
 
-function finalize_cft!(scheme::SLoopTNR)
+CFT_Finalizer = Finalizer(finalize_cftdata!, CFTData)
+
+function finalize_cftdata!(scheme::SLoopTNR) # TODO: remove this
     tr_norm = trnorm_2x2(scheme.T)
     scheme.T /= tr_norm^0.25
     Tflip = flip(scheme.T, (1, 2, 3, 4))
@@ -141,16 +143,9 @@ function finalize_cft!(scheme::SLoopTNR)
     return data
 end
 
-# central charge finalize
-function finalize_central_charge!(scheme::TNRScheme)
-    n = finalize!(scheme)
-    return central_charge(scheme, n)
-end
-
 # TODO: add Finalizers for CFT and central charge
 two_by_two_Finalizer = Finalizer(finalize_two_by_two!, Float64)
 
-
 # Finalizer for ground state degeneracy
 function finalize_groundstatedegeneracy!(scheme::TNRScheme)
     finalize!(scheme)
diff --git a/test/models.jl b/test/models.jl
index 37b98afd..3f9e3c5e 100644
--- a/test/models.jl
+++ b/test/models.jl
@@ -44,8 +44,8 @@ end
     T_KT = classical_XY(U1Irrep, XY_βc + 0.1, 8)
     scheme = LoopTNR(T_KT)
     data = run!(scheme, truncrank(16), maxiter(20))
-    cft = cft_data(scheme, [sqrt(2), 2 * sqrt(2), 0])
-    central_charge = cft["c"]
+    cft = CFTData(scheme)
+    central_charge = cft.central_charge
     @test central_charge ≈ 1.0 atol = 1.0e-2
     @info "Obtained central charge:\n$central_charge."
 
@@ -53,8 +53,8 @@ end
     T_sym = classical_XY(U1Irrep, XY_βc - 0.1, 8)
     scheme = LoopTNR(T_sym)
     data = run!(scheme, truncrank(16), maxiter(20))
-    cft = cft_data(scheme, [sqrt(2), 2 * sqrt(2), 0])
-    central_charge = cft["c"]
+    cft = CFTData(scheme)
+    central_charge = cft.central_charge
     @test central_charge ≈ 0.0 atol = 1.0e-13
     @info "Obtained central charge:\n$central_charge."
 
@@ -62,8 +62,8 @@ end
     T_KT = classical_XY(CU1Irrep, XY_βc + 0.1, 8)
     scheme = LoopTNR(T_KT)
     data = run!(scheme, truncrank(16), maxiter(20))
-    cft = cft_data(scheme, [sqrt(2), 2 * sqrt(2), 0])
-    central_charge = cft["c"]
+    cft = CFTData(scheme)
+    central_charge = cft.central_charge
     @test central_charge ≈ 1.0 atol = 1.0e-2
     @info "Obtained central charge:\n$central_charge."
 
@@ -71,8 +71,8 @@ end
     T_sym = classical_XY(CU1Irrep, XY_βc - 0.1, 8)
     scheme = LoopTNR(T_sym)
     data = run!(scheme, truncrank(16), maxiter(20))
-    cft = cft_data(scheme, [sqrt(2), 2 * sqrt(2), 0])
-    central_charge = cft["c"]
+    cft = CFTData(scheme)
+    central_charge = cft.central_charge
     @test central_charge ≈ 0.0 atol = 1.0e-13
     @info "Obtained central charge:\n$central_charge."
 end
diff --git a/test/schemes.jl b/test/schemes.jl
index 37b71c23..2bfe4240 100644
--- a/test/schemes.jl
+++ b/test/schemes.jl
@@ -11,7 +11,7 @@ const f_benchmark3D = -3.507
 
 function cft_finalize!(scheme)
     finalize!(scheme)
-    return cft_data(scheme)
+    return CFTData(scheme)
 end
 
 # TRG
@@ -26,7 +26,7 @@ end
     scheme = TRG(T)
     run!(scheme, truncrank(24), maxiter(10))
 
-    cft = cft_data(scheme)[2:end]
+    cft = sort(CFTData(scheme; shape = [1, 1, 0]).scaling_dimensions[2:end]; by = abs) .|> real
 
     @test cft[1] ≈ ising_cft_exact[1] rtol = 2.0e-4
     @test cft[2] ≈ ising_cft_exact[2] rtol = 1.0e-2
@@ -65,7 +65,7 @@ end
     scheme = BTRG(T)
     run!(scheme, truncrank(24), maxiter(10))
 
-    cft = cft_data(scheme)[2:end]
+    cft = sort(CFTData(scheme; shape = [1, 1, 0]).scaling_dimensions[2:end]; by = abs) .|> real
 
     @test cft[1] ≈ ising_cft_exact[1] rtol = 3.0e-4
     @test cft[2] ≈ ising_cft_exact[2] rtol = 2.0e-2
@@ -102,7 +102,7 @@ end
     scheme = HOTRG(T)
     run!(scheme, truncrank(16), maxiter(4))
 
-    cft = cft_data(scheme)[2:end]
+    cft = sort(CFTData(scheme; shape = [1, 1, 0]).scaling_dimensions[2:end]; by = abs) .|> real
 
     @test cft[1] ≈ ising_cft_exact[1] rtol = 6.0e-4
     @test cft[2] ≈ ising_cft_exact[2] rtol = 1.0e-2
@@ -139,7 +139,7 @@ end
     scheme = ATRG(T)
     run!(scheme, truncrank(24), maxiter(3))
 
-    cft = cft_data(scheme)[2:end]
+    cft = sort(CFTData(scheme; shape = [1, 1, 0]).scaling_dimensions[2:end]; by = abs) .|> real
 
     @test cft[1] ≈ ising_cft_exact[1] rtol = 1.0e-2
     @test cft[2] ≈ ising_cft_exact[2] rtol = 1.0e-2
@@ -185,7 +185,7 @@ end
     run!(scheme, truncrank(12), maxiter(10))
 
     for shape in [[1, 4, 1], [sqrt(2), 2 * sqrt(2), 0]]
-        cft = cft_data(scheme, shape)
+        cft = sort(CFTData(scheme; shape = shape).scaling_dimensions; by = abs)
         d1, d2 = real(cft[Z2Irrep(1)][1]), real(cft[Z2Irrep(0)][2])
         @info "Obtained lowest scaling dimensions:\n$(d1), $(d2)."
         @test d1 ≈ ising_cft_exact[1] rtol = 5.0e-4
@@ -193,7 +193,7 @@ end
     end
 
     for shape in [[1, 8, 1], [4 / sqrt(10), 2 * sqrt(10), 2 / sqrt(10)]]
-        cft = cft_data(scheme, shape, truncrank(12), trunctol(atol = 1.0e-10))
+        cft = sort(CFTData(scheme; shape = shape, trunc = truncrank(16), truncentanglement = trunctol(atol = 1.0e-10)).scaling_dimensions; by = real)
         d1, d2 = real(cft[Z2Irrep(1)][1]), real(cft[Z2Irrep(0)][2])
         @info "Obtained lowest scaling dimensions:\n$(d1), $(d2)."
         @test d1 ≈ ising_cft_exact[1] rtol = 1.0e-3