diff --git a/src/qibocal/protocols/state_tomography.py b/src/qibocal/protocols/state_tomography.py
index 6094de710..70fbcb74e 100644
--- a/src/qibocal/protocols/state_tomography.py
+++ b/src/qibocal/protocols/state_tomography.py
@@ -20,8 +20,8 @@
 
 BASIS = ["X", "Y", "Z"]
 """Single qubit measurement basis."""
-SIMULATED_DENSITY_MATRIX = "ideal"
-"""Filename for simulated density matrix."""
+CIRCUIT_PATH = "circuit.json"
+"""Path where circuit is stored."""
 
 
 @dataclass
@@ -53,25 +53,30 @@ def __post_init__(self):
 class StateTomographyData(Data):
     """Tomography data"""
 
-    ideal: dict[tuple[QubitId, str], np.float64] = field(default_factory=dict)
-    """Ideal samples measurements."""
+    targets: dict[QubitId, int]
+    """Store targets order."""
+    circuit: Circuit
+    """Circuit where tomography will be executed."""
     data: dict[tuple[QubitId, str], np.int64] = field(default_factory=dict)
     """Hardware measurements."""
 
+    @property
+    def params(self) -> dict:
+        """Convert non-arrays attributes into dict."""
+        params = super().params
+        params.pop("circuit")
+        return params
+
     def save(self, path):
-        self._to_npz(path, DATAFILE)
-        np.savez(
-            path / f"{SIMULATED_DENSITY_MATRIX}.npz",
-            **{json.dumps(i): self.ideal[i] for i in self.ideal},
-        )
+        super().save(path)
+        (path / CIRCUIT_PATH).write_text(json.dumps(self.circuit.raw))
 
     @classmethod
     def load(cls, path):
-        instance = cls()
-        instance.data = super().load_data(path, DATAFILE)
-        instance.ideal = super().load_data(path, SIMULATED_DENSITY_MATRIX)
-
-        return instance
+        circuit = Circuit.from_dict(json.loads((path / CIRCUIT_PATH).read_text()))
+        data = super().load_data(path, DATAFILE)
+        params = super().load_params(path, DATAFILE)
+        return cls(data=data, circuit=circuit, targets=params["targets"])
 
 
 @dataclass
@@ -94,7 +99,6 @@ def _acquisition(
     params: StateTomographyParameters, platform: Platform, targets: list[QubitId]
 ) -> StateTomographyData:
     """Acquisition protocol for single qubit state tomography experiment."""
-
     if params.circuit is None:
         params.circuit = Circuit(len(targets))
 
@@ -102,7 +106,9 @@ def _acquisition(
     backend.platform = platform
     transpiler = dummy_transpiler(backend)
 
-    data = StateTomographyData()
+    data = StateTomographyData(
+        circuit=params.circuit, targets={target: i for i, target in enumerate(targets)}
+    )
 
     for basis in BASIS:
         basis_circuit = deepcopy(params.circuit)
@@ -111,7 +117,7 @@ def _acquisition(
             for i in range(len(targets)):
                 basis_circuit.add(getattr(gates, basis)(i).basis_rotation())
 
-        basis_circuit.add(gates.M(i) for i in range(len(targets)))
+        basis_circuit.add(gates.M(*range(len(targets))))
         _, results = execute_transpiled_circuit(
             basis_circuit,
             targets,
@@ -127,9 +133,6 @@ def _acquisition(
                     samples=np.array(results.samples()).T[i],
                 ),
             )
-            data.ideal[target, basis] = np.array(
-                NumpyBackend().execute_circuit(basis_circuit, nshots=10000).samples()
-            ).T[i]
     return data
 
 
@@ -140,30 +143,27 @@ def _fit(data: StateTomographyData) -> StateTomographyResults:
     target_density_matrix_real = {}
     target_density_matrix_imag = {}
     fid = {}
-    for qubit in data.qubits:
-        x_exp, y_exp, z_exp = (
-            1 - 2 * np.mean(data[qubit, basis].samples) for basis in BASIS
+    circuit = data.circuit
+    circuit.density_matrix = True
+    total_density_matrix = NumpyBackend().execute_circuit(circuit=circuit).state()
+    for i, qubit in enumerate(data.targets):
+        traced_qubits = [q for q in range(len(data.qubits)) if q != i]
+        target_density_matrix = NumpyBackend().partial_trace_density_matrix(
+            total_density_matrix, traced_qubits, len(data.qubits)
         )
-        density_matrix = 0.5 * (
+        x_exp = 1 - 2 * np.mean(data[qubit, "X"].samples)
+        y_exp = 1 - 2 * np.mean(data[qubit, "Y"].samples)
+        z_exp = 1 - 2 * np.mean(data[qubit, "Z"].samples)
+        measured_density_matrix = 0.5 * (
             matrices.I + matrices.X * x_exp + matrices.Y * y_exp + matrices.Z * z_exp
         )
-        measured_density_matrix_real[qubit] = np.real(density_matrix).tolist()
-        measured_density_matrix_imag[qubit] = np.imag(density_matrix).tolist()
+        measured_density_matrix_real[qubit] = np.real(measured_density_matrix).tolist()
+        measured_density_matrix_imag[qubit] = np.imag(measured_density_matrix).tolist()
 
-        x_theory, y_theory, z_theory = (
-            1 - 2 * np.mean(data.ideal[qubit, basis]) for basis in BASIS
-        )
-        target_density_matrix = 0.5 * (
-            matrices.I
-            + matrices.X * x_theory
-            + matrices.Y * y_theory
-            + matrices.Z * z_theory
-        )
         target_density_matrix_real[qubit] = np.real(target_density_matrix).tolist()
         target_density_matrix_imag[qubit] = np.imag(target_density_matrix).tolist()
         fid[qubit] = fidelity(
-            np.array(measured_density_matrix_real[qubit])
-            + 1.0j * np.array(measured_density_matrix_imag[qubit]),
+            measured_density_matrix,
             target_density_matrix,
         )