diff --git a/doc/sphinx/source/data/dataset-naming-convention.rst b/doc/sphinx/source/data/dataset-naming-convention.rst
index 51f4036b82..f0bee2dab7 100644
--- a/doc/sphinx/source/data/dataset-naming-convention.rst
+++ b/doc/sphinx/source/data/dataset-naming-convention.rst
@@ -13,11 +13,11 @@ constructed following this [Backus–Naur form]::
 
   <valid dataset name> ::= <set name> "_" <observable name>
 
-  <experiment> ::= "ATLAS" | "BCDMS" | "CDF" | "CHORUS" | "CMS" | "D0" | "DYE605" | "DYE866" 
+  <experiment> ::= "ATLAS" | "BCDMS" | "CDF" | "CHORUS" | "CMS" | "D0" | "DYE605" | "DYE866"
                 | "DYE906" | "EMC" | "H1" | "HERA" | "LHCB" | "NMC" | "NNPDF" | "NUTEV" | "SLAC"
                 | "ZEUS"
 
-  <process> ::= "1JET" | "2JET" | "CC" | "DY" | "INTEG" | "NC" | "PH" | "POS" | "SINGLETOP" 
+  <process> ::= "1JET" | "2JET" | "CC" | "DY" | "INTEG" | "NC" | "PH" | "POS" | "SHP" | "SINGLETOP"
               | "TTBAR" | "WCHARM" | "WJ" | "WPWM" | "Z0" | "Z0J"
 
   <energy> ::= <integer> <unit> | <integer> "P" <integer> <unit>  | "NOTFIXED"
@@ -91,6 +91,7 @@ Processes
 - `NC`: DIS neutral-current
 - `POS`: auxiliary dataset for positivity constraints; only valid for
   `NNPDF` experiment
+- `SHP`: single hadron production
 - `TTBAR`: top–anti-top production
 - `WM`: production of a single negatively-charged lepton (charged current
   off-shell Drell–Yan)
diff --git a/nnpdf_data/nnpdf_data/commondata/PHENIX-2009_SHP_200GEV_PI0/data.yaml b/nnpdf_data/nnpdf_data/commondata/PHENIX-2009_SHP_200GEV_PI0/data.yaml
new file mode 100644
index 0000000000..e8daaf85ab
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/PHENIX-2009_SHP_200GEV_PI0/data.yaml
@@ -0,0 +1,13 @@
+data_central:
+- 0.00051
+- 0.00096
+- 0.00039
+- -0.00023
+- 0.0006
+- 0.0002
+- 0.0013
+- 0.0026
+- -0.0039
+- 0.0096
+- 0.008
+- 0.061
diff --git a/nnpdf_data/nnpdf_data/commondata/PHENIX-2009_SHP_200GEV_PI0/filter.py b/nnpdf_data/nnpdf_data/commondata/PHENIX-2009_SHP_200GEV_PI0/filter.py
new file mode 100644
index 0000000000..642445b7bc
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/PHENIX-2009_SHP_200GEV_PI0/filter.py
@@ -0,0 +1,93 @@
+from pathlib import Path
+
+import pandas as pd
+import yaml
+
+
+def read_data(path_rawdata: str) -> pd.DataFrame:
+    data = yaml.safe_load(Path(path_rawdata).read_text())
+
+    ptvals = data["independent_variables"][0]["values"]
+    asym_obs = data["dependent_variables"][0]["values"]
+
+    concatenated_table = []
+    for i in range(len(ptvals)):
+        # Compute pT mid-value if not given
+        pt_mid = (float(ptvals[i]["high"]) + float(ptvals[i]["low"])) / 2
+        pt_midval = ptvals[i].get("value", pt_mid)
+
+        concatenated_table.append(
+            pd.DataFrame(
+                {
+                    "pT_low": [ptvals[i]["low"]],
+                    "pT": [pt_midval],
+                    "pT_high": [ptvals[i]["high"]],
+                    "asym": [asym_obs[i]["value"]],
+                    "stat_err": [asym_obs[i]["errors"][0]["symerror"]],
+                    "lumi_err": [asym_obs[i]["errors"][1]["symerror"]],
+                    "spol_err": [asym_obs[i]["errors"][2]["symerror"].removesuffix("%")],
+                }
+            )
+        )
+
+    return pd.concat(concatenated_table, ignore_index=True)
+
+
+def dump_data(df_table: pd.DataFrame) -> None:
+    # Dump central data into Yaml file
+    data_central = []
+    for i in range(len(df_table["asym"])):
+        data_central.append(float(df_table.loc[i, "asym"]))
+
+    with open("data.yaml", "w") as file:
+        yaml.dump({"data_central": data_central}, file, sort_keys=False)
+
+    # Dump the kinematics into Yaml file
+    kinematics = []
+    for i in range(len(df_table["asym"])):
+        kin_value = {
+            "pT": {
+                "min": float(df_table.loc[i, "pT_low"]),
+                "mid": float(df_table.loc[i, "pT"]),
+                "max": float(df_table.loc[i, "pT_high"]),
+            },
+            "eta": {"min": -0.35, "mid": 0.0, "max": 0.35},
+        }
+        kinematics.append(kin_value)
+
+    with open("kinematics.yaml", "w") as file:
+        yaml.dump({"bins": kinematics}, file, sort_keys=False)
+
+    # Dump the uncertainties into Yaml file
+    errors = []
+    for i in range(len(df_table)):
+        error_per_bin = {
+            "stat": float(df_table.loc[i, "stat_err"]),
+            "sys_lumi": float(df_table.loc[i, "lumi_err"]),
+            "sys_pol": abs(data_central[i]) * float(df_table.loc[i, "spol_err"]) / 100.0,
+        }
+        errors.append(error_per_bin)
+
+    error_definition = {
+        "stat": {"description": "Statistical uncertainty", "treatment": "ADD", "type": "UNCORR"},
+        "sys_lumi": {
+            "description": "Systematic uncertainties due to luminosity",
+            "treatment": "MULT",
+            "type": "CORR",
+        },
+        "sys_pol": {
+            "description": "Systematic uncertainties due to polarization",
+            "treatment": "MULT",
+            "type": "CORR",
+        },
+    }
+
+    with open("uncertainties.yaml", "w") as file:
+        yaml.dump({"definitions": error_definition, "bins": errors}, file, sort_keys=False)
+
+    return
+
+
+if __name__ == "__main__":
+    df_table = read_data("./rawdata/HEPData-ins1282448-v1-Table_7.yaml")
+    dump_data(df_table=df_table)
diff --git a/nnpdf_data/nnpdf_data/commondata/PHENIX-2009_SHP_200GEV_PI0/kinematics.yaml b/nnpdf_data/nnpdf_data/commondata/PHENIX-2009_SHP_200GEV_PI0/kinematics.yaml
new file mode 100644
index 0000000000..c3a6860af0
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/PHENIX-2009_SHP_200GEV_PI0/kinematics.yaml
@@ -0,0 +1,97 @@
+bins:
+- pT:
+    min: 1.0
+    mid: 1.3
+    max: 1.5
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 1.5
+    mid: 1.75
+    max: 2.0
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 2.0
+    mid: 2.23
+    max: 2.5
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 2.5
+    mid: 2.72
+    max: 3.0
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 3.0
+    mid: 3.22
+    max: 3.5
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 3.5
+    mid: 3.72
+    max: 4.0
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 4.0
+    mid: 4.39
+    max: 5.0
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 5.0
+    mid: 5.4
+    max: 6.0
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 6.0
+    mid: 6.41
+    max: 7.0
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 7.0
+    mid: 7.74
+    max: 9.0
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 9.0
+    mid: 10.0
+    max: 12.0
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 12.0
+    mid: 13.1
+    max: 15.0
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
diff --git a/nnpdf_data/nnpdf_data/commondata/PHENIX-2009_SHP_200GEV_PI0/metadata.yaml b/nnpdf_data/nnpdf_data/commondata/PHENIX-2009_SHP_200GEV_PI0/metadata.yaml
new file mode 100644
index 0000000000..d5e0e4f1e3
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/PHENIX-2009_SHP_200GEV_PI0/metadata.yaml
@@ -0,0 +1,45 @@
+# Generalia
+setname: "PHENIX-2009_SHP_200GEV_PI0"
+
+version: 1
+version_comment: "Initial implementation"
+
+# References
+iNSPIRE:
+  url: "https://inspirehep.net/literature/1282448"
+hepdata:
+  url: "https://www.hepdata.net/record/ins1282448"
+  version: 1
+
+nnpdf_metadata:
+  nnpdf31_process: "SHP" # Single Hadron Production
+  experiment: "PHENIX"
+
+implemented_observables:
+  - observable_name: "ALL"
+    observable:
+      description: "Double helicity asymmetry in inclusive pi^0 production in polarized p+p collisions at sqrt(s)=200 GeV"
+      label: "$A_{LL}$"
+      units: ""
+    process_type: "SHP_ASY"
+    ndata: 12
+    tables: [7]
+    npoints: [12] # List of datapoints per table
+
+    # Plotting information
+    plotting:
+      kinematics_override: identity # TODO
+      dataset_label: "PHENIX ALL"
+      y_label: "$A_{LL}(p_T)$"
+      plot_x: pT
+    kinematic_coverage: [pT, eta]
+
+    kinematics:
+      variables:
+        pT: { description: "Transverse momentum", label: "$p_T$", units: "GeV" }
+        eta: { description: "pi^0 pseudorapidity", label: r"$\eta$", units: "" }
+      file: kinematics.yaml
+
+    data_central: data.yaml
+    data_uncertainties:
+      - uncertainties.yaml
diff --git a/nnpdf_data/nnpdf_data/commondata/PHENIX-2009_SHP_200GEV_PI0/rawdata/HEPData-ins1282448-v1-Table_7.yaml b/nnpdf_data/nnpdf_data/commondata/PHENIX-2009_SHP_200GEV_PI0/rawdata/HEPData-ins1282448-v1-Table_7.yaml
new file mode 100644
index 0000000000..f8ab996936
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/PHENIX-2009_SHP_200GEV_PI0/rawdata/HEPData-ins1282448-v1-Table_7.yaml
@@ -0,0 +1,81 @@
+dependent_variables:
+- header: {name: ASYM(LL)}
+  qualifiers:
+  - {name: RE, value: P P --> PI0 < GAMMA GAMMA > X}
+  - {name: SQRT(S), units: GeV, value: '200.0'}
+  values:
+  - errors:
+    - {label: stat, symerror: 0.00085}
+    - {label: 'sys,rel.lumi.', symerror: 0.00035}
+    - {label: 'sys,pol.', symerror: 3.4%}
+    value: 0.00051
+  - errors:
+    - {label: stat, symerror: 0.00055}
+    - {label: 'sys,rel.lumi.', symerror: 0.00033}
+    - {label: 'sys,pol.', symerror: 3.5%}
+    value: 0.00096
+  - errors:
+    - {label: stat, symerror: 0.00058}
+    - {label: 'sys,rel.lumi.', symerror: 0.00034}
+    - {label: 'sys,pol.', symerror: 3.5%}
+    value: 0.00039
+  - errors:
+    - {label: stat, symerror: 0.00074}
+    - {label: 'sys,rel.lumi.', symerror: 0.00036}
+    - {label: 'sys,pol.', symerror: 3.4%}
+    value: -0.00023
+  - errors:
+    - {label: stat, symerror: 0.0011}
+    - {label: 'sys,rel.lumi.', symerror: 0.0004}
+    - {label: 'sys,pol.', symerror: 3.2%}
+    value: 0.0006
+  - errors:
+    - {label: stat, symerror: 0.0015}
+    - {label: 'sys,rel.lumi.', symerror: 0.00041}
+    - {label: 'sys,pol.', symerror: 3.1%}
+    value: 0.0002
+  - errors:
+    - {label: stat, symerror: 0.0018}
+    - {label: 'sys,rel.lumi.', symerror: 0.00043}
+    - {label: 'sys,pol.', symerror: 3.1%}
+    value: 0.0013
+  - errors:
+    - {label: stat, symerror: 0.0035}
+    - {label: 'sys,rel.lumi.', symerror: 0.00045}
+    - {label: 'sys,pol.', symerror: 3.0%}
+    value: 0.0026
+  - errors:
+    - {label: stat, symerror: 0.0061}
+    - {label: 'sys,rel.lumi.', symerror: 0.00045}
+    - {label: 'sys,pol.', symerror: 2.9%}
+    value: -0.0039
+  - errors:
+    - {label: stat, symerror: 0.0085}
+    - {label: 'sys,rel.lumi.', symerror: 0.00045}
+    - {label: 'sys,pol.', symerror: 2.9%}
+    value: 0.0096
+  - errors:
+    - {label: stat, symerror: 0.018}
+    - {label: 'sys,rel.lumi.', symerror: 0.00058}
+    - {label: 'sys,pol.', symerror: 3.3%}
+    value: 0.008
+  - errors:
+    - {label: stat, symerror: 0.069}
+    - {label: 'sys,rel.lumi.', symerror: 0.001}
+    - {label: 'sys,pol.', symerror: 3.0%}
+    value: 0.061
+independent_variables:
+- header: {name: PT(PI0), units: GEV}
+  values:
+  - {high: 1.5, low: 1.0, value: 1.3}
+  - {high: 2.0, low: 1.5}
+  - {high: 2.5, low: 2.0, value: 2.23}
+  - {high: 3.0, low: 2.5, value: 2.72}
+  - {high: 3.5, low: 3.0, value: 3.22}
+  - {high: 4.0, low: 3.5, value: 3.72}
+  - {high: 5.0, low: 4.0, value: 4.39}
+  - {high: 6.0, low: 5.0, value: 5.4}
+  - {high: 7.0, low: 6.0, value: 6.41}
+  - {high: 9.0, low: 7.0, value: 7.74}
+  - {high: 12.0, low: 9.0, value: 10.0}
+  - {high: 15.0, low: 12.0, value: 13.1}
diff --git a/nnpdf_data/nnpdf_data/commondata/PHENIX-2009_SHP_200GEV_PI0/uncertainties.yaml b/nnpdf_data/nnpdf_data/commondata/PHENIX-2009_SHP_200GEV_PI0/uncertainties.yaml
new file mode 100644
index 0000000000..c5cc374156
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/PHENIX-2009_SHP_200GEV_PI0/uncertainties.yaml
@@ -0,0 +1,50 @@
+definitions:
+  stat:
+    description: Statistical uncertainty
+    treatment: ADD
+    type: UNCORR
+  sys_lumi:
+    description: Systematic uncertainties due to luminosity
+    treatment: MULT
+    type: CORR
+  sys_pol:
+    description: Systematic uncertainties due to polarization
+    treatment: MULT
+    type: CORR
+bins:
+- stat: 0.00085
+  sys_lumi: 0.00035
+  sys_pol: 1.734e-05
+- stat: 0.00055
+  sys_lumi: 0.00033
+  sys_pol: 3.3600000000000004e-05
+- stat: 0.00058
+  sys_lumi: 0.00034
+  sys_pol: 1.365e-05
+- stat: 0.00074
+  sys_lumi: 0.00036
+  sys_pol: 7.82e-06
+- stat: 0.0011
+  sys_lumi: 0.0004
+  sys_pol: 1.92e-05
+- stat: 0.0015
+  sys_lumi: 0.00041
+  sys_pol: 6.2e-06
+- stat: 0.0018
+  sys_lumi: 0.00043
+  sys_pol: 4.03e-05
+- stat: 0.0035
+  sys_lumi: 0.00045
+  sys_pol: 7.8e-05
+- stat: 0.0061
+  sys_lumi: 0.00045
+  sys_pol: 0.00011309999999999998
+- stat: 0.0085
+  sys_lumi: 0.00045
+  sys_pol: 0.0002784
+- stat: 0.018
+  sys_lumi: 0.00058
+  sys_pol: 0.000264
+- stat: 0.069
+  sys_lumi: 0.001
+  sys_pol: 0.00183
diff --git a/nnpdf_data/nnpdf_data/commondata/PHENIX-2013_SHP_510GEV_PI0/data.yaml b/nnpdf_data/nnpdf_data/commondata/PHENIX-2013_SHP_510GEV_PI0/data.yaml
new file mode 100644
index 0000000000..78c85297e6
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/PHENIX-2013_SHP_510GEV_PI0/data.yaml
@@ -0,0 +1,15 @@
+data_central:
+- -6.0e-05
+- -0.0009
+- 0.00035
+- -0.00047
+- 0.0011
+- 0.0008
+- 0.00039
+- 0.00158
+- 0.0013
+- 0.00303
+- 0.0064
+- 0.0078
+- 0.0001
+- 0.0086
diff --git a/nnpdf_data/nnpdf_data/commondata/PHENIX-2013_SHP_510GEV_PI0/filter.py b/nnpdf_data/nnpdf_data/commondata/PHENIX-2013_SHP_510GEV_PI0/filter.py
new file mode 100644
index 0000000000..a92cb51d5d
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/PHENIX-2013_SHP_510GEV_PI0/filter.py
@@ -0,0 +1,88 @@
+from pathlib import Path
+
+import pandas as pd
+import yaml
+
+
+def read_data(path_rawdata: str) -> pd.DataFrame:
+    data = yaml.safe_load(Path(path_rawdata).read_text())
+
+    ptvals = data["independent_variables"][0]["values"]
+    pt_midval = data["dependent_variables"][0]["values"]
+    asym_obs = data["dependent_variables"][1]["values"]
+
+    concatenated_table = []
+    for i in range(len(ptvals)):
+        concatenated_table.append(
+            pd.DataFrame(
+                {
+                    "pT_low": [ptvals[i]["low"]],
+                    "pT": [pt_midval[i]["value"]],
+                    "pT_high": [ptvals[i]["high"]],
+                    "asym": [asym_obs[i]["value"]],
+                    "stat_err": [asym_obs[i]["errors"][0]["symerror"]],
+                }
+            )
+        )
+
+    return pd.concat(concatenated_table, ignore_index=True)
+
+
+def dump_data(df_table: pd.DataFrame) -> None:
+    # Dump central data into Yaml file
+    data_central = []
+    for i in range(len(df_table["asym"])):
+        data_central.append(float(df_table.loc[i, "asym"]))
+
+    with open("data.yaml", "w") as file:
+        yaml.dump({"data_central": data_central}, file, sort_keys=False)
+
+    # Dump the kinematics into Yaml file
+    kinematics = []
+    for i in range(len(df_table["asym"])):
+        kin_value = {
+            "pT": {
+                "min": float(df_table.loc[i, "pT_low"]),
+                "mid": float(df_table.loc[i, "pT"]),
+                "max": float(df_table.loc[i, "pT_high"]),
+            },
+            "eta": {"min": -0.35, "mid": 0.0, "max": 0.35},
+        }
+        kinematics.append(kin_value)
+
+    with open("kinematics.yaml", "w") as file:
+        yaml.dump({"bins": kinematics}, file, sort_keys=False)
+
+    # Dump the uncertainties into Yaml file
+    errors = []
+    for i in range(len(df_table)):
+        error_per_bin = {
+            "stat": float(df_table.loc[i, "stat_err"]),
+            "sys_lumi": 3.6e-4,
+            "sys_pol": (abs(data_central[i]) * 6.5) / 100.0,
+        }
+        errors.append(error_per_bin)
+
+    error_definition = {
+        "stat": {"description": "Statistical uncertainty", "treatment": "ADD", "type": "UNCORR"},
+        "sys_lumi": {
+            "description": "Systematic uncertainties due to luminosity",
+            "treatment": "MULT",
+            "type": "CORR",
+        },
+        "sys_pol": {
+            "description": "Systematic uncertainties due to beam polarization",
+            "treatment": "MULT",
+            "type": "PHENIX2013POL",
+        },
+    }
+
+    with open("uncertainties.yaml", "w") as file:
+        yaml.dump({"definitions": error_definition, "bins": errors}, file, sort_keys=False)
+
+    return
+
+
+if __name__ == "__main__":
+    df_table = read_data("./rawdata/HEPData-ins1396712-v1-Figure_3.yaml")
+    dump_data(df_table=df_table)
diff --git a/nnpdf_data/nnpdf_data/commondata/PHENIX-2013_SHP_510GEV_PI0/kinematics.yaml b/nnpdf_data/nnpdf_data/commondata/PHENIX-2013_SHP_510GEV_PI0/kinematics.yaml
new file mode 100644
index 0000000000..ac6893b670
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/PHENIX-2013_SHP_510GEV_PI0/kinematics.yaml
@@ -0,0 +1,113 @@
+bins:
+- pT:
+    min: 2.0
+    mid: 2.28
+    max: 2.5
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 2.5
+    mid: 2.76
+    max: 3.0
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 3.0
+    mid: 3.25
+    max: 3.5
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 3.5
+    mid: 3.74
+    max: 4.0
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 4.0
+    mid: 4.24
+    max: 4.5
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 4.5
+    mid: 4.74
+    max: 5.0
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 5.0
+    mid: 5.45
+    max: 6.0
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 6.0
+    mid: 6.45
+    max: 7.0
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 7.0
+    mid: 7.45
+    max: 8.0
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 8.0
+    mid: 8.45
+    max: 9.0
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 9.0
+    mid: 9.45
+    max: 10.0
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 10.0
+    mid: 10.82
+    max: 12.0
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 12.0
+    mid: 13.14
+    max: 14.0
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
+- pT:
+    min: 14.0
+    mid: 16.62
+    max: 18.0
+  eta:
+    min: -0.35
+    mid: 0.0
+    max: 0.35
diff --git a/nnpdf_data/nnpdf_data/commondata/PHENIX-2013_SHP_510GEV_PI0/metadata.yaml b/nnpdf_data/nnpdf_data/commondata/PHENIX-2013_SHP_510GEV_PI0/metadata.yaml
new file mode 100644
index 0000000000..79d290b485
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/PHENIX-2013_SHP_510GEV_PI0/metadata.yaml
@@ -0,0 +1,42 @@
+# Generalia
+setname: "PHENIX-2013_SHP_510GEV_PI0"
+
+version: 1
+version_comment: "Initial implementation"
+
+# References
+iNSPIRE:
+  url: "https://inspirehep.net/literature/1396712"
+hepdata:
+  url: "https://www.hepdata.net/record/ins1396712"
+  version: 1
+
+nnpdf_metadata:
+  nnpdf31_process: "SHP" # Single Hadron Production
+  experiment: "PHENIX"
+
+implemented_observables:
+  - observable_name: "ALL"
+    observable:
+      description: "Inclusive cross section and double-helicity asymmetry for π0 production at midrapidity in collisions at 200 GeV"
+      label: "$A_{LL}$"
+      units: ""
+    process_type: "SHP_ASY"
+    ndata: 14
+    tables: [2]
+    npoints: [14] # List of datapoints per table
+    plotting:
+      kinematics_override: identity
+      dataset_label: "PHENIX ALL"
+      y_label: "$A_{LL}(p_T)$"
+      plot_x: pT
+    kinematic_coverage: [pT, eta]
+    kinematics:
+      variables:
+        pT: { description: "Transverse momentum", label: "$p_T$", units: "GeV" }
+        eta: { description: "π0 pseudorapidity", label: r"$\eta$", units: ""}
+      file: kinematics.yaml
+
+    data_central: data.yaml
+    data_uncertainties:
+      - uncertainties.yaml
diff --git a/nnpdf_data/nnpdf_data/commondata/PHENIX-2013_SHP_510GEV_PI0/rawdata/HEPData-ins1396712-v1-Figure_3.yaml b/nnpdf_data/nnpdf_data/commondata/PHENIX-2013_SHP_510GEV_PI0/rawdata/HEPData-ins1396712-v1-Figure_3.yaml
new file mode 100644
index 0000000000..d482e846b2
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/PHENIX-2013_SHP_510GEV_PI0/rawdata/HEPData-ins1396712-v1-Figure_3.yaml
@@ -0,0 +1,82 @@
+independent_variables:
+- header: {name: '$p_T$ (GeV/$c$)'}
+  values:
+  - {low: 2.00, high: 2.50}
+  - {low: 2.50, high: 3.00}
+  - {low: 3.00, high: 3.50}
+  - {low: 3.50, high: 4.00}
+  - {low: 4.00, high: 4.50}
+  - {low: 4.50, high: 5.00}
+  - {low: 5.00, high: 6.00}
+  - {low: 6.00, high: 7.00}
+  - {low: 7.00, high: 8.00}
+  - {low: 8.00, high: 9.00}
+  - {low: 9.00, high: 10.00}
+  - {low: 10.00, high: 12.00}
+  - {low: 12.00, high: 14.00}
+  - {low: 14.00, high: 18.00}
+dependent_variables:
+- header: {name: '$\langle p_T \rangle$ (GeV/$c$)'}
+  qualifiers:
+  - {name: 'particle', value: '$\pi^0$'}
+  values:
+  - value: 2.280
+  - value: 2.760
+  - value: 3.250
+  - value: 3.740
+  - value: 4.240
+  - value: 4.740
+  - value: 5.450
+  - value: 6.450
+  - value: 7.450
+  - value: 8.450
+  - value: 9.450
+  - value: 10.820
+  - value: 13.140
+  - value: 16.620
+- header: {name: '$A_{LL}$'}
+  qualifiers:
+  - {name: 'particle', value: '$\pi^0$'}
+  values:
+  - value: -6.00e-5
+    errors:
+    - {symerror: 0.0012000, label: 'tot'}
+  - value: -9.00e-4
+    errors:
+    - {symerror: 8.60e-4, label: 'tot'}
+  - value: 3.50e-4
+    errors:
+    - {symerror: 7.60e-4, label: 'tot'}
+  - value: -4.70e-4
+    errors:
+    - {symerror: 7.60e-4, label: 'tot'}
+  - value: 0.00110
+    errors:
+    - {symerror: 8.2e-4, label: 'tot'}
+  - value: 8.0e-4
+    errors:
+    - {symerror: 9.3e-4, label: 'tot'}
+  - value: 3.90e-4
+    errors:
+    - {symerror: 8.30e-4, label: 'tot'}
+  - value: 0.00158
+    errors:
+    - {symerror: 0.00117, label: 'tot'}
+  - value: 0.00130
+    errors:
+    - {symerror: 0.00165, label: 'tot'}
+  - value: 0.00303
+    errors:
+    - {symerror: 0.00231, label: 'tot'}
+  - value: 0.0064
+    errors:
+    - {symerror: 0.0032, label: 'tot'}
+  - value: 0.0078
+    errors:
+    - {symerror: 0.0035, label: 'tot'}
+  - value: 1.00e-4
+    errors:
+    - {symerror: 0.005500, label: 'tot'}
+  - value: 0.0086
+    errors:
+    - {symerror: 0.0110, label: 'tot'}
diff --git a/nnpdf_data/nnpdf_data/commondata/PHENIX-2013_SHP_510GEV_PI0/uncertainties.yaml b/nnpdf_data/nnpdf_data/commondata/PHENIX-2013_SHP_510GEV_PI0/uncertainties.yaml
new file mode 100644
index 0000000000..8f849fa568
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/PHENIX-2013_SHP_510GEV_PI0/uncertainties.yaml
@@ -0,0 +1,56 @@
+definitions:
+  stat:
+    description: Statistical uncertainty
+    treatment: ADD
+    type: UNCORR
+  sys_lumi:
+    description: Systematic uncertainties due to luminosity
+    treatment: MULT
+    type: CORR
+  sys_pol:
+    description: Systematic uncertainties due to beam polarization
+    treatment: MULT
+    type: PHENIX2013POL
+bins:
+- stat: 0.0012
+  sys_lumi: 0.00036
+  sys_pol: 3.9e-06
+- stat: 0.00086
+  sys_lumi: 0.00036
+  sys_pol: 5.85e-05
+- stat: 0.00076
+  sys_lumi: 0.00036
+  sys_pol: 2.275e-05
+- stat: 0.00076
+  sys_lumi: 0.00036
+  sys_pol: 3.055e-05
+- stat: 0.00082
+  sys_lumi: 0.00036
+  sys_pol: 7.15e-05
+- stat: 0.00093
+  sys_lumi: 0.00036
+  sys_pol: 5.2000000000000004e-05
+- stat: 0.00083
+  sys_lumi: 0.00036
+  sys_pol: 2.535e-05
+- stat: 0.00117
+  sys_lumi: 0.00036
+  sys_pol: 0.0001027
+- stat: 0.00165
+  sys_lumi: 0.00036
+  sys_pol: 8.45e-05
+- stat: 0.00231
+  sys_lumi: 0.00036
+  sys_pol: 0.00019695000000000002
+- stat: 0.0032
+  sys_lumi: 0.00036
+  sys_pol: 0.00041600000000000003
+- stat: 0.0035
+  sys_lumi: 0.00036
+  sys_pol: 0.000507
+- stat: 0.0055
+  sys_lumi: 0.00036
+  sys_pol: 6.5000000000000004e-06
+- stat: 0.011
+  sys_lumi: 0.00036
+  sys_pol: 0.000559
diff --git a/nnpdf_data/nnpdf_data/commondata/STAR-2006_SHP_200GEV_PI0/data.yaml b/nnpdf_data/nnpdf_data/commondata/STAR-2006_SHP_200GEV_PI0/data.yaml
new file mode 100644
index 0000000000..1d04db9fd2
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/STAR-2006_SHP_200GEV_PI0/data.yaml
@@ -0,0 +1,7 @@
+data_central:
+- -0.023
+- 0.004
+- 0.048
+- -0.001
+- 0.038
+- -0.045
diff --git a/nnpdf_data/nnpdf_data/commondata/STAR-2006_SHP_200GEV_PI0/filter.py b/nnpdf_data/nnpdf_data/commondata/STAR-2006_SHP_200GEV_PI0/filter.py
new file mode 100644
index 0000000000..fdd8c754fd
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/STAR-2006_SHP_200GEV_PI0/filter.py
@@ -0,0 +1,82 @@
+from pathlib import Path
+
+import pandas as pd
+import yaml
+
+
+def read_data(path_rawdata: str) -> pd.DataFrame:
+    data = yaml.safe_load(Path(path_rawdata).read_text())
+
+    ptvals = data["independent_variables"][0]["values"]
+    asym_obs = data["dependent_variables"][0]["values"]
+
+    concatenated_table = []
+    for i in range(len(ptvals)):
+        concatenated_table.append(
+            pd.DataFrame(
+                {
+                    "pT": [ptvals[i]["value"]],
+                    "asym": [asym_obs[i]["value"]],
+                    "stat_err": [asym_obs[i]["errors"][0]["symerror"]],
+                    "syst_err": [asym_obs[i]["errors"][1]["symerror"]],
+                }
+            )
+        )
+
+    return pd.concat(concatenated_table, ignore_index=True)
+
+
+def dump_data(df_table: pd.DataFrame) -> None:
+    # Dump central data into Yaml file
+    data_central = []
+    for i in range(len(df_table["asym"])):
+        data_central.append(float(df_table.loc[i, "asym"]))
+
+    with open("data.yaml", "w") as file:
+        yaml.dump({"data_central": data_central}, file, sort_keys=False)
+
+    # Dump the kinematics into Yaml file
+    kinematics = []
+    for i in range(len(df_table["asym"])):
+        kin_value = {
+            "pT": {"min": None, "mid": float(df_table.loc[i, "pT"]), "max": None},
+            "eta": {"min": 0.8, "mid": 1.4, "max": 2.0},
+        }
+        kinematics.append(kin_value)
+
+    with open("kinematics.yaml", "w") as file:
+        yaml.dump({"bins": kinematics}, file, sort_keys=False)
+
+    # Dump the uncertainties into Yaml file
+    errors = []
+    for i in range(len(df_table)):
+        error_per_bin = {
+            "stat": float(df_table.loc[i, "stat_err"]),
+            "syst": float(df_table.loc[i, "syst_err"]),
+            "sys_pol": abs(data_central[i]) * 6 / 100.0,
+        }
+        errors.append(error_per_bin)
+
+    error_definition = {
+        "stat": {"description": "Statistical uncertainty", "treatment": "ADD", "type": "UNCORR"},
+        "syst": {
+            "description": "Total systematic uncertainties",
+            "treatment": "MULT",
+            "type": "CORR",
+        },
+        "sys_pol": {
+            "description": "Systematic uncertainties due to beam polarization",
+            "treatment": "MULT",
+            "type": "STAR2006POL",
+        },
+    }
+
+    with open("uncertainties.yaml", "w") as file:
+        yaml.dump({"definitions": error_definition, "bins": errors}, file, sort_keys=False)
+
+    return
+
+
+if __name__ == "__main__":
+    df_table = read_data("./rawdata/HEPData-ins1253360-v1-Figure_7_Data.yaml")
+    dump_data(df_table=df_table)
diff --git a/nnpdf_data/nnpdf_data/commondata/STAR-2006_SHP_200GEV_PI0/kinematics.yaml b/nnpdf_data/nnpdf_data/commondata/STAR-2006_SHP_200GEV_PI0/kinematics.yaml
new file mode 100644
index 0000000000..87f13208a4
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/STAR-2006_SHP_200GEV_PI0/kinematics.yaml
@@ -0,0 +1,49 @@
+bins:
+- pT:
+    min: null
+    mid: 5.57692
+    max: null
+  eta:
+    min: 0.8
+    mid: 1.4
+    max: 2.0
+- pT:
+    min: null
+    mid: 6.51397
+    max: null
+  eta:
+    min: 0.8
+    mid: 1.4
+    max: 2.0
+- pT:
+    min: null
+    mid: 7.48737
+    max: null
+  eta:
+    min: 0.8
+    mid: 1.4
+    max: 2.0
+- pT:
+    min: null
+    mid: 8.4765
+    max: null
+  eta:
+    min: 0.8
+    mid: 1.4
+    max: 2.0
+- pT:
+    min: null
+    mid: 9.44489
+    max: null
+  eta:
+    min: 0.8
+    mid: 1.4
+    max: 2.0
+- pT:
+    min: null
+    mid: 10.7949
+    max: null
+  eta:
+    min: 0.8
+    mid: 1.4
+    max: 2.0
diff --git a/nnpdf_data/nnpdf_data/commondata/STAR-2006_SHP_200GEV_PI0/metadata.yaml b/nnpdf_data/nnpdf_data/commondata/STAR-2006_SHP_200GEV_PI0/metadata.yaml
new file mode 100644
index 0000000000..46abcc9659
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/STAR-2006_SHP_200GEV_PI0/metadata.yaml
@@ -0,0 +1,45 @@
+# Generalia
+setname: "STAR-2006_SHP_200GEV_PI0"
+
+version: 1
+version_comment: "Initial implementation"
+
+# References
+iNSPIRE:
+  url: "https://inspirehep.net/literature/1253360"
+hepdata:
+  url: "https://www.hepdata.net/record/ins1253360"
+  version: 1
+
+nnpdf_metadata:
+  nnpdf31_process: "SHP" # Single Hadron Production
+  experiment: "STAR"
+
+implemented_observables:
+  - observable_name: "ALL"
+    observable:
+      description: "Neutral pion spin asymmetries at intermediate pseudorapidity in polarized proton collisions sqrt(s)=200 GeV"
+      label: "$A_{LL}$"
+      units: ""
+    process_type: "SHP_ASY"
+    ndata: 6
+    tables: [7]
+    npoints: [6] # List of datapoints per table
+
+    # Plotting information
+    plotting:
+      kinematics_override: identity # TODO
+      dataset_label: "STAR ALL"
+      y_label: "$A_{LL}(p_T)$"
+      plot_x: pT
+    kinematic_coverage: [pT, eta]
+
+    kinematics:
+      variables:
+        pT: { description: "Transverse momentum", label: "$p_T$", units: "GeV" }
+        eta: { description: "Pseudorapidity", label: r"$\eta$", units: "" }
+      file: kinematics.yaml
+
+    data_central: data.yaml
+    data_uncertainties:
+      - uncertainties.yaml
diff --git a/nnpdf_data/nnpdf_data/commondata/STAR-2006_SHP_200GEV_PI0/rawdata/HEPData-ins1253360-v1-Figure_7_Data.yaml b/nnpdf_data/nnpdf_data/commondata/STAR-2006_SHP_200GEV_PI0/rawdata/HEPData-ins1253360-v1-Figure_7_Data.yaml
new file mode 100644
index 0000000000..73fccd7b0a
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/STAR-2006_SHP_200GEV_PI0/rawdata/HEPData-ins1253360-v1-Figure_7_Data.yaml
@@ -0,0 +1,38 @@
+independent_variables:
+- header: {name: '$p_{T} [GeV/c]$'}
+  values:
+  - {value: 5.57692}
+  - {value: 6.51397}
+  - {value: 7.48737}
+  - {value: 8.4765}
+  - {value: 9.44489}
+  - {value: 10.7949}
+dependent_variables:
+- header: {name: '$A_{LL}$'}
+  qualifiers:
+  - {name: '-', value: 'STAR Data'}
+  values:
+  - value: -0.023
+    errors:
+    - {symerror: 0.019, label: 'stat error'}
+    - {symerror: 0.002, label: 'sys error'}
+  - value: 0.004
+    errors:
+    - {symerror: 0.021, label: 'stat error'}
+    - {symerror: 0.004, label: 'sys error'}
+  - value: 0.048
+    errors:
+    - {symerror: 0.028, label: 'stat error'}
+    - {symerror: 0.006, label: 'sys error'}
+  - value: -0.001
+    errors:
+    - {symerror: 0.045, label: 'stat error'}
+    - {symerror: 0.007, label: 'sys error'}
+  - value: 0.038
+    errors:
+    - {symerror: 0.060, label: 'stat error'}
+    - {symerror: 0.009, label: 'sys error'}
+  - value: -0.045
+    errors:
+    - {symerror: 0.072, label: 'stat error'}
+    - {symerror: 0.009, label: 'sys error'}
diff --git a/nnpdf_data/nnpdf_data/commondata/STAR-2006_SHP_200GEV_PI0/uncertainties.yaml b/nnpdf_data/nnpdf_data/commondata/STAR-2006_SHP_200GEV_PI0/uncertainties.yaml
new file mode 100644
index 0000000000..051f353edd
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/STAR-2006_SHP_200GEV_PI0/uncertainties.yaml
@@ -0,0 +1,32 @@
+definitions:
+  stat:
+    description: Statistical uncertainty
+    treatment: ADD
+    type: UNCORR
+  syst:
+    description: Total systematic uncertainties
+    treatment: MULT
+    type: CORR
+  sys_pol:
+    description: Systematic uncertainties due to beam polarization
+    treatment: MULT
+    type: STAR2006POL
+bins:
+- stat: 0.019
+  syst: 0.002
+  sys_pol: 0.0013800000000000002
+- stat: 0.021
+  syst: 0.004
+  sys_pol: 0.00024
+- stat: 0.028
+  syst: 0.006
+  sys_pol: 0.00288
+- stat: 0.045
+  syst: 0.007
+  sys_pol: 6.0e-05
+- stat: 0.06
+  syst: 0.009
+  sys_pol: 0.00228
+- stat: 0.072
+  syst: 0.009
+  sys_pol: 0.0027
diff --git a/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/data_highrap.yaml b/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/data_highrap.yaml
new file mode 100644
index 0000000000..e3541a51a4
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/data_highrap.yaml
@@ -0,0 +1,5 @@
+data_central:
+- -0.00191
+- -0.00085
+- -0.00177
+- 0.00015
diff --git a/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/data_lowrap.yaml b/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/data_lowrap.yaml
new file mode 100644
index 0000000000..c6388bec27
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/data_lowrap.yaml
@@ -0,0 +1,5 @@
+data_central:
+- -0.00154
+- 0.0021
+- -0.0014
+- 0.00076
diff --git a/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/filter.py b/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/filter.py
new file mode 100644
index 0000000000..1d1f8e6017
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/filter.py
@@ -0,0 +1,89 @@
+from pathlib import Path
+
+import pandas as pd
+import yaml
+
+
+def read_data(path_rawdata: str) -> pd.DataFrame:
+    data = yaml.safe_load(Path(path_rawdata).read_text())
+
+    ptvals = data["independent_variables"][0]["values"]
+    asym_obs = data["dependent_variables"][0]["values"]
+
+    concatenated_table = []
+    for i in range(len(ptvals)):
+        concatenated_table.append(
+            pd.DataFrame(
+                {
+                    "pT_low": [ptvals[i]["low"]],
+                    "pT_high": [ptvals[i]["high"]],
+                    "asym": [asym_obs[i]["value"]],
+                    "stat_err": [asym_obs[i]["errors"][0]["symerror"]],
+                    "syst_err": [asym_obs[i]["errors"][1]["symerror"]],
+                }
+            )
+        )
+
+    return pd.concat(concatenated_table, ignore_index=True)
+
+
+def dump_data(df_table: pd.DataFrame, tableid: int) -> None:
+    suffix = "lowrap.yaml" if tableid == 1 else "highrap.yaml"
+    # Dump central data into Yaml file
+    data_central = []
+    for i in range(len(df_table["asym"])):
+        data_central.append(float(df_table.loc[i, "asym"]))
+
+    with open(f"data_{suffix}", "w") as file:
+        yaml.dump({"data_central": data_central}, file, sort_keys=False)
+
+    # Dump the kinematics into Yaml file
+    kinematics = []
+    for i in range(len(df_table["asym"])):
+        kin_value = {
+            "pT": {
+                "min": float(df_table.loc[i, "pT_low"]),
+                "mid": (float(df_table.loc[i, "pT_high"]) + float(df_table.loc[i, "pT_low"])) / 2,
+                "max": float(df_table.loc[i, "pT_high"]),
+            },
+            "eta": {"min": 3.15, "mid": 3.525, "max": 3.90},
+        }
+        kinematics.append(kin_value)
+
+    with open(f"kinematics_{suffix}", "w") as file:
+        yaml.dump({"bins": kinematics}, file, sort_keys=False)
+
+    # Dump the uncertainties into Yaml file
+    errors = []
+    for i in range(len(df_table)):
+        error_per_bin = {
+            "stat": float(df_table.loc[i, "stat_err"]),
+            "syst": float(df_table.loc[i, "syst_err"]),
+            "sys_pol": abs(data_central[i]) * 6.7 / 100.0,
+        }
+        errors.append(error_per_bin)
+
+    error_definition = {
+        "stat": {"description": "Statistical uncertainty", "treatment": "ADD", "type": "UNCORR"},
+        "syst": {
+            "description": "Total systematic uncertainties",
+            "treatment": "MULT",
+            "type": "CORR",
+        },
+        "sys_pol": {
+            "description": "Systematic uncertainties due to beam polarization",
+            "treatment": "MULT",
+            "type": "STAR2013POL",
+        },
+    }
+
+    with open(f"uncertainties_{suffix}", "w") as file:
+        yaml.dump({"definitions": error_definition, "bins": errors}, file, sort_keys=False)
+
+    return
+
+
+if __name__ == "__main__":
+    for tabid in [1, 2]:
+        df_table = read_data(f"./rawdata/HEPData-ins1674826-v1-Table_{tabid}.yaml")
+        dump_data(df_table=df_table, tableid=tabid)
diff --git a/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/kinematics_highrap.yaml b/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/kinematics_highrap.yaml
new file mode 100644
index 0000000000..0da835fd68
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/kinematics_highrap.yaml
@@ -0,0 +1,33 @@
+bins:
+- pT:
+    min: 3.7
+    mid: 3.91
+    max: 4.12
+  eta:
+    min: 3.15
+    mid: 3.525
+    max: 3.9
+- pT:
+    min: 4.48
+    mid: 4.73
+    max: 4.98
+  eta:
+    min: 3.15
+    mid: 3.525
+    max: 3.9
+- pT:
+    min: 5.33
+    mid: 5.62
+    max: 5.91
+  eta:
+    min: 3.15
+    mid: 3.525
+    max: 3.9
+- pT:
+    min: 6.71
+    mid: 7.08
+    max: 7.45
+  eta:
+    min: 3.15
+    mid: 3.525
+    max: 3.9
diff --git a/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/kinematics_lowrap.yaml b/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/kinematics_lowrap.yaml
new file mode 100644
index 0000000000..6bc318481f
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/kinematics_lowrap.yaml
@@ -0,0 +1,33 @@
+bins:
+- pT:
+    min: 2.37
+    mid: 2.5
+    max: 2.63
+  eta:
+    min: 3.15
+    mid: 3.525
+    max: 3.9
+- pT:
+    min: 3.15
+    mid: 3.33
+    max: 3.51
+  eta:
+    min: 3.15
+    mid: 3.525
+    max: 3.9
+- pT:
+    min: 3.89
+    mid: 4.11
+    max: 4.33
+  eta:
+    min: 3.15
+    mid: 3.525
+    max: 3.9
+- pT:
+    min: 5.08
+    mid: 5.37
+    max: 5.66
+  eta:
+    min: 3.15
+    mid: 3.525
+    max: 3.9
diff --git a/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/metadata.yaml b/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/metadata.yaml
new file mode 100644
index 0000000000..001a969d1f
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/metadata.yaml
@@ -0,0 +1,66 @@
+# Generalia
+setname: "STAR-2013_SHP_510GEV_PI0"
+
+version: 1
+version_comment: "Initial implementation"
+
+# References
+iNSPIRE:
+  url: "https://inspirehep.net/literature/1674826"
+hepdata:
+  url: "https://www.hepdata.net/record/ins1674826"
+  version: 1
+
+nnpdf_metadata:
+  nnpdf31_process: "SHP" # Single Hadron Production
+  experiment: "STAR"
+
+implemented_observables:
+  - observable_name: "ALL-LOWRAP"
+    observable:
+      description: "Neutral pion spin asymmetries at low pseudorapidity in polarized proton collisions sqrt(s)=510 GeV"
+      label: "$A_{LL}$"
+      units: ""
+    process_type: "SHP_ASY"
+    ndata: 4
+    tables: [2]
+    npoints: [4] # List of datapoints per table
+    # Plotting information
+    plotting:
+      kinematics_override: identity
+      dataset_label: "STAR ALL"
+      y_label: "$A_{LL}(p_T)$"
+      plot_x: pT
+    kinematic_coverage: [pT, eta]
+    kinematics:
+      variables:
+        pT: { description: "Transverse momentum", label: "$p_T$", units: "GeV" }
+        eta: { description: "Pseudorapidity", label: r"$\eta$", units: "" }
+      file: kinematics_lowrap.yaml
+    data_central: data_lowrap.yaml
+    data_uncertainties:
+      - uncertainties_lowrap.yaml
+  - observable_name: "ALL-HIGHRAP"
+    observable:
+      description: "Neutral pion spin asymmetries at high pseudorapidity in polarized proton collisions sqrt(s)=510 GeV"
+      label: "$A_{LL}$"
+      units: ""
+    process_type: "SHP_ASY"
+    ndata: 4
+    tables: [1]
+    npoints: [4] # List of datapoints per table
+    # Plotting information
+    plotting:
+      kinematics_override: identity
+      dataset_label: "STAR ALL"
+      y_label: "$A_{LL}(p_T)$"
+      plot_x: pT
+    kinematic_coverage: [pT, eta]
+    kinematics:
+      variables:
+        pT: { description: "Transverse momentum", label: "$p_T$", units: "GeV" }
+        eta: { description: "Pseudorapidity", label: r"$\eta$", units: "" }
+      file: kinematics_highrap.yaml
+    data_central: data_highrap.yaml
+    data_uncertainties:
+      - uncertainties_highrap.yaml
diff --git a/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/rawdata/HEPData-ins1674826-v1-Table_1.yaml b/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/rawdata/HEPData-ins1674826-v1-Table_1.yaml
new file mode 100644
index 0000000000..14d5dfd6d1
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/rawdata/HEPData-ins1674826-v1-Table_1.yaml
@@ -0,0 +1,28 @@
+independent_variables:
+- header: {name: '$p_T$ range (GeV/$c$)'}
+  values:
+  - {low: 2.37, high: 2.63}
+  - {low: 3.15, high: 3.51}
+  - {low: 3.89, high: 4.33}
+  - {low: 5.08, high: 5.66}
+dependent_variables:
+- header: {name: '$A_{LL}(\pi^0)$'}
+  qualifiers:
+  - {name: '$p_T$ range (GeV/$c$)', value: 'Pseudorapidity range $3.15 < \eta < 3.90$'}
+  values:
+  - value: -0.00154
+    errors:
+    - {symerror: 0.0022, label: 'Stat. error (mb GeV$^{-2}c^3$)'}
+    - {symerror: 0.00021, label: 'Syst. error (mb GeV$^{-2}c^3$)'}
+  - value: 0.00210
+    errors:
+    - {symerror: 0.0021, label: 'Stat. error (mb GeV$^{-2}c^3$)'}
+    - {symerror: 0.00021, label: 'Syst. error (mb GeV$^{-2}c^3$)'}
+  - value: -0.00140
+    errors:
+    - {symerror: 0.0022, label: 'Stat. error (mb GeV$^{-2}c^3$)'}
+    - {symerror: 0.00021, label: 'Syst. error (mb GeV$^{-2}c^3$)'}
+  - value: 0.00076
+    errors:
+    - {symerror: 0.0023, label: 'Stat. error (mb GeV$^{-2}c^3$)'}
+    - {symerror: 0.00021, label: 'Syst. error (mb GeV$^{-2}c^3$)'}
diff --git a/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/rawdata/HEPData-ins1674826-v1-Table_2.yaml b/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/rawdata/HEPData-ins1674826-v1-Table_2.yaml
new file mode 100644
index 0000000000..13ed77e1fb
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/rawdata/HEPData-ins1674826-v1-Table_2.yaml
@@ -0,0 +1,28 @@
+independent_variables:
+- header: {name: '$p_T$ range (GeV/$c$)'}
+  values:
+  - {low: 3.70, high: 4.12}
+  - {low: 4.48, high: 4.98}
+  - {low: 5.33, high: 5.91}
+  - {low: 6.71, high: 7.45}
+dependent_variables:
+- header: {name: '$A_{LL}(\pi^0)$'}
+  qualifiers:
+  - {name: '$p_T$ range (GeV/$c$)', value: 'Pseudorapidity range $2.65 < \eta < 3.15$'}
+  values:
+  - value: -0.00191
+    errors:
+    - {symerror: 0.0022, label: 'Stat. error (mb GeV$^{-2}c^3$)'}
+    - {symerror: 0.00034, label: 'Syst. error (mb GeV$^{-2}c^3$)'}
+  - value: -0.00085
+    errors:
+    - {symerror: 0.0024, label: 'Stat. error (mb GeV$^{-2}c^3$)'}
+    - {symerror: 0.00032, label: 'Syst. error (mb GeV$^{-2}c^3$)'}
+  - value: -0.00177
+    errors:
+    - {symerror: 0.0025, label: 'Stat. error (mb GeV$^{-2}c^3$)'}
+    - {symerror: 0.00031, label: 'Syst. error (mb GeV$^{-2}c^3$)'}
+  - value: 0.00015
+    errors:
+    - {symerror: 0.0024, label: 'Stat. error (mb GeV$^{-2}c^3$)'}
+    - {symerror: 0.00029, label: 'Syst. error (mb GeV$^{-2}c^3$)'}
diff --git a/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/uncertainties_highrap.yaml b/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/uncertainties_highrap.yaml
new file mode 100644
index 0000000000..432e675bc9
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/uncertainties_highrap.yaml
@@ -0,0 +1,26 @@
+definitions:
+  stat:
+    description: Statistical uncertainty
+    treatment: ADD
+    type: UNCORR
+  syst:
+    description: Total systematic uncertainties
+    treatment: MULT
+    type: CORR
+  sys_pol:
+    description: Systematic uncertainties due to beam polarization
+    treatment: MULT
+    type: STAR2013POL
+bins:
+- stat: 0.0022
+  syst: 0.00034
+  sys_pol: 0.00012797
+- stat: 0.0024
+  syst: 0.00032
+  sys_pol: 5.6949999999999995e-05
+- stat: 0.0025
+  syst: 0.00031
+  sys_pol: 0.00011859000000000001
+- stat: 0.0024
+  syst: 0.00029
+  sys_pol: 1.005e-05
diff --git a/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/uncertainties_lowrap.yaml b/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/uncertainties_lowrap.yaml
new file mode 100644
index 0000000000..96adb7298c
--- /dev/null
+++ b/nnpdf_data/nnpdf_data/commondata/STAR-2013_SHP_510GEV_PI0/uncertainties_lowrap.yaml
@@ -0,0 +1,26 @@
+definitions:
+  stat:
+    description: Statistical uncertainty
+    treatment: ADD
+    type: UNCORR
+  syst:
+    description: Total systematic uncertainties
+    treatment: MULT
+    type: CORR
+  sys_pol:
+    description: Systematic uncertainties due to beam polarization
+    treatment: MULT
+    type: STAR2013POL
+bins:
+- stat: 0.0022
+  syst: 0.00021
+  sys_pol: 0.00010317999999999998
+- stat: 0.0021
+  syst: 0.00021
+  sys_pol: 0.0001407
+- stat: 0.0022
+  syst: 0.00021
+  sys_pol: 9.379999999999999e-05
+- stat: 0.0023
+  syst: 0.00021
+  sys_pol: 5.092e-05
diff --git a/validphys2/src/validphys/commondataparser.py b/validphys2/src/validphys/commondataparser.py
index 60440869e7..d64d3c0d65 100644
--- a/validphys2/src/validphys/commondataparser.py
+++ b/validphys2/src/validphys/commondataparser.py
@@ -109,6 +109,7 @@ def _quick_yaml_load(filepath):
     "JET": ("$\\eta$", "$p_T^2 (GeV^2)$", "$\\sqrt{s} (GeV)$"),
     "PHT": ("$\\eta_\\gamma$", "$E_{T,\\gamma}^2 (GeV^2)$", "$\\sqrt{s} (GeV)$"),
     "SIA": ("$z$", "$Q^2 (GeV^2)$", "$y$"),
+    "SHP_ASY": ("$\\eta$", "$p_T (GeV)$", "$\\sqrt{s} (GeV)$"),
     "JET_POL": ("$\\eta$", "$p_T^2 (GeV^2)$", "$\\sqrt{s} (GeV)$"),
     "DIJET_POL": ("$\\m_{1,2} (GeV)", "$\\eta_1$", "$\\eta_2$"),
 }
@@ -137,6 +138,7 @@ def _quick_yaml_load(filepath):
     "DYP": "Fixed-Target Drell-Yan",
     "JET_POL": "Inclusive Jet longitudinal double-spin asymmetry",
     "DIJET_POL": "Dijets longitudinal double-spin asymmetry",
+    "SHP_ASY": "double spin asymmetry in single hadron production",
 }
 
 
diff --git a/validphys2/src/validphys/filters.py b/validphys2/src/validphys/filters.py
index 780f00eb09..0f47bd08db 100644
--- a/validphys2/src/validphys/filters.py
+++ b/validphys2/src/validphys/filters.py
@@ -27,6 +27,7 @@
     "DIJET": ("eta", "m_12", "sqrts"),
     "PHT": ("eta_gamma", "E_{T,gamma)2", "sqrts"),
     "INC": ("0", "mu2", "sqrts"),
+    "SHP_ASY": ("eta", "pT", "sqrts"),
     "EWK_RAP": ("etay", "M2", "sqrts"),
     "EWK_RAP_ASY": ("etay", "M2", "sqrts"),
     "EWK_PT": ("p_T", "M2", "sqrts"),
diff --git a/validphys2/src/validphys/process_options.py b/validphys2/src/validphys/process_options.py
index 151d359103..ec3cd7da5c 100644
--- a/validphys2/src/validphys/process_options.py
+++ b/validphys2/src/validphys/process_options.py
@@ -132,7 +132,7 @@ def _dis_xq2map(kin_info):
     """
     x = kin_info.get_one_of("k1", _Vars.x)
     if "k2" in kin_info._kins:
-        q2 = kin_info.get_one_of("k2")**2
+        q2 = kin_info.get_one_of("k2") ** 2
     else:
         q2 = kin_info.get_one_of(_Vars.Q2)
     return x, q2
@@ -150,6 +150,18 @@ def _jets_xq2map(kin_info):
     return np.clip(x, a_min=None, a_max=1, out=x), np.concatenate((q2, q2))
 
 
+def _shp_xq2map(kin_info):
+    # Then compute x, Q2
+    pT = kin_info[_Vars.pT]
+    ratio = pT / kin_info[_Vars.sqrts]
+    rap = kin_info[_Vars.eta]
+    x1 = 2 * ratio * np.exp(rap)
+    x2 = 2 * ratio * np.exp(-rap)
+    q2 = pT * pT
+    x = np.concatenate((x1, x2))
+    return np.clip(x, a_min=None, a_max=1, out=x), np.concatenate((q2, q2))
+
+
 def _dijets_xq2map(kin_info):
     # Here we can have either ystar or ydiff, but in either case we need to do the same
     ylab_1 = kin_info.get_one_of(_Vars.ystar, _Vars.ydiff, _Vars.eta_1, _Vars.abs_eta_1)
@@ -258,6 +270,13 @@ def _dybosonpt_xq2map(kin_dict):
     xq2map_function=_jets_xq2map,
 )
 
+SHP = _Process(
+    "SHP",
+    "Single Hadron Production",
+    accepted_variables=(_Vars.eta, _Vars.pT, _Vars.sqrts),
+    xq2map_function=_shp_xq2map,
+)
+
 DIJET = _Process(
     "DIJET",
     "DiJets production",
@@ -359,6 +378,7 @@ def _dybosonpt_xq2map(kin_dict):
     "DIS_POL": dataclasses.replace(DIS, name="DIS_POL"),
     "JET": JET,
     "DIJET": DIJET,
+    "SHP_ASY": SHP,
     "HQP_YQ": HQP_YQ,
     "HQP_YQQ": dataclasses.replace(HQP_YQ, name="HQP_YQQ"),
     "HQP_PTQ": HQP_PTQ,