diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 85707c1f..c8bb3b19 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -6,6 +6,7 @@ CHANGELOG
 ================
 
 * We have decided internally to revert the versioning schema back to SemVer, starting at ``1.0.0``, in order to improve better backward and forward compatibility support for plugins.
+* Fix bug related to Probabilistic History Matching result reader. The shape of result was wrong.
 
 
 2024.2 (2024-09-10)
diff --git a/src/alfasim_sdk/result_reader/aggregator.py b/src/alfasim_sdk/result_reader/aggregator.py
index ef2339b5..b0c926a5 100644
--- a/src/alfasim_sdk/result_reader/aggregator.py
+++ b/src/alfasim_sdk/result_reader/aggregator.py
@@ -1965,11 +1965,14 @@ def read_history_matching_result(
     if hm_type not in ("HM-deterministic", "HM-probabilistic"):
         raise ValueError(f"history matching of type `{hm_type}` not supported")
 
+    slicer: Callable[[int], Any]
     if hm_type == "HM-deterministic":
         dataset_key = HISTORY_MATCHING_DETERMINISTIC_DSET_NAME
+        slicer = lambda index: index
     else:
         assert hm_type == "HM-probabilistic"
         dataset_key = HISTORY_MATCHING_PROBABILISTIC_DSET_NAME
+        slicer = lambda index: (slice(None), index)
 
     with open_result_file(metadata.result_directory) as result_file:
         if not result_file:
@@ -1980,11 +1983,11 @@ def read_history_matching_result(
         result_map = {}
         if hm_result_key is None:
             for key, meta in metadata.hm_items.items():
-                result_map[key] = result[meta.data_index]
+                result_map[key] = result[slicer(meta.data_index)]
         else:
             meta = metadata.hm_items.get(hm_result_key)
             if meta is not None:
-                result_map[hm_result_key] = result[meta.data_index]
+                result_map[hm_result_key] = result[slicer(meta.data_index)]
 
         return result_map
 
diff --git a/tests/conftest.py b/tests/conftest.py
index a9fe44d8..615bc102 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -463,7 +463,7 @@ def hm_probabilistic_results_dir(datadir: Path) -> Path:
 
     result_dir = datadir / "main-HM-probabilistic"
     probabilistic_result = np.array(
-        [[0.1, 0.22, 1.0, 0.8, 0.55], [3.0, 6.0, 5.1, 4.7, 6.3]]
+        [[0.1, 3.0], [0.22, 6.0], [1.0, 5.1], [0.8, 4.7], [0.55, 6.3]]
     )
     historic_data_curves = {
         "observed_curve_1": np.array([[0.1, 0.5, 0.9], [1.1, 2.2, 3.3]]),
diff --git a/tests/results/test_aggregator.py b/tests/results/test_aggregator.py
index 96d4f44e..2d037255 100644
--- a/tests/results/test_aggregator.py
+++ b/tests/results/test_aggregator.py
@@ -511,13 +511,18 @@ def test_read_history_matching_result_data(
     Check reading the result of both HM type analysis. Both results are available simultaneously by
     the means of the fixtures, but only one is used at a time.
     """
+    import numpy as np
+
     # Setup.
     if hm_type == "HM-probabilistic":
-        expected_results = ([0.1, 0.22, 1.0, 0.8, 0.55], [3.0, 6.0, 5.1, 4.7, 6.3])
+        expected_results = {
+            "parametric_var_1": np.array([0.1, 0.22, 1.0, 0.8, 0.55]),
+            "parametric_var_2": np.array([3.0, 6.0, 5.1, 4.7, 6.3]),
+        }
         results_dir = hm_probabilistic_results_dir
     else:
         assert hm_type == "HM-deterministic"
-        expected_results = (0.1, 3.2)
+        expected_results = {"parametric_var_1": 0.1, "parametric_var_2": 3.2}
         results_dir = hm_deterministic_results_dir
 
     metadata = read_history_matching_metadata(results_dir)
@@ -527,13 +532,20 @@ def test_read_history_matching_result_data(
         metadata, hm_type=hm_type, hm_result_key="parametric_var_1"
     )
     assert len(result) == 1
-    assert result["parametric_var_1"] == pytest.approx(expected_results[0])
+    assert len(result) == 1
+    assert result["parametric_var_1"] == pytest.approx(
+        expected_results["parametric_var_1"]
+    )
 
     # Read the result of all entries.
     result = read_history_matching_result(metadata, hm_type=hm_type)
     assert len(result) == 2
-    assert result["parametric_var_1"] == pytest.approx(expected_results[0])
-    assert result["parametric_var_2"] == pytest.approx(expected_results[1])
+    assert result["parametric_var_1"] == pytest.approx(
+        expected_results["parametric_var_1"]
+    )
+    assert result["parametric_var_2"] == pytest.approx(
+        expected_results["parametric_var_2"]
+    )
 
     # Unexistent result key, result should be empty.
     result = read_history_matching_result(