diff --git a/pyerrors/input/__init__.py b/pyerrors/input/__init__.py
index 7a835029..5c7496e9 100644
--- a/pyerrors/input/__init__.py
+++ b/pyerrors/input/__init__.py
@@ -1,2 +1,3 @@
 from .input import *
 from . import bdio
+from . import hadrons
diff --git a/pyerrors/input/hadrons.py b/pyerrors/input/hadrons.py
new file mode 100644
index 00000000..b63f385d
--- /dev/null
+++ b/pyerrors/input/hadrons.py
@@ -0,0 +1,56 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+import os
+import h5py
+import numpy as np
+from ..pyerrors import Obs
+from ..correlators import Corr
+
+def read_meson_hd5(path, filestem, ens_id, meson='meson_0'):
+    """Read hadrons meson hdf5 file and extract the meson labeled 'meson'
+
+    Parameters
+    -----------------
+    path -- path to the files to read
+    filestem -- namestem of the files to read
+    ens_id -- name of the ensemble, required for internal bookkeeping
+    meson -- label of the meson to be extracted, standard value meson_0 which corresponds to the pseudoscalar pseudoscalar two-point function.
+    """
+    ls = []
+    for (dirpath, dirnames, filenames) in os.walk(path):
+        ls.extend(filenames)
+        break
+
+    # Clean up file list
+    files = []
+    for l in ls:
+        if l.startswith(filestem):
+            files.append(l)
+
+    # Sort according to configuration number
+    get_cnfg_number = lambda x : int(x[len(filestem) + 1:-3])
+    files.sort(key=get_cnfg_number)
+
+    # Check that configurations are evenly spaced
+    cnfg_numbers = []
+    for l in files:
+        cnfg_numbers.append(get_cnfg_number(l))
+
+    if not all(np.diff(cnfg_numbers)==np.diff(cnfg_numbers)[0]):
+        raise Exception('Configurations are not evenly spaced.')
+
+    corr_data = []
+    for hd5_file in files:
+        file = h5py.File(path + '/' + hd5_file, "r")
+
+        raw_data = list(file['meson/' + meson + '/corr'])
+        real_data = [o[0] for o in raw_data]
+        corr_data.append(real_data)
+    corr_data = np.array(corr_data)
+
+    l_obs = []
+    for c in corr_data.T:
+        l_obs.append(Obs([c], [ens_id]))
+
+    return Corr(l_obs)