diff --git a/.pylintrc b/.pylintrc
index ef3ba76d..0ae1ec5d 100644
--- a/.pylintrc
+++ b/.pylintrc
@@ -130,7 +130,8 @@ disable=
 	global-statement,
 	too-many-public-method,
 	too-many-ancestors,
-	too-many-positional-arguments
+	too-many-positional-arguments,
+	too-many-statements
 
 # Enable the message, report, category or checker with the given id(s). You can
 # either give multiple identifier separated by comma (,) or put this option
diff --git a/docs/source/fitting_code/running.md b/docs/source/fitting_code/running.md
index dafc4ec3..59c36345 100644
--- a/docs/source/fitting_code/running.md
+++ b/docs/source/fitting_code/running.md
@@ -41,13 +41,14 @@ theory_path:
 
 ```
 ### Theory specifications
-The perturbative order of the QCD theory prediction (LO or NLO) should be specified using ``order``.
-``use_quad`` should be set to ``True`` for a fit with quadratic corrections, ``use_t0`` controls the use
+The default perturbative order of the theory prediction is set by the key ``default_order``. Orders may also be specified
+per datset, see [here](./example.html#datasets-to-consider-and-coefficients-to-fit) for more details.
+The order in the EFT expansion should be specified by setting ``use_quad`` to either ``True`` or ``False`` to include quadratic or only linear corrections respectively. The option ``use_t0`` controls the use
 of the ``t0`` prescription and ``use_theory_covmat`` specifies whether or not to use the theory covariance matrix
 which can be specified in the theory files.
 
 ```yaml
-order: NLO
+default_order: LO
 use_quad: False
 use_t0: False
 use_theory_covmat: True
@@ -91,22 +92,32 @@ n_samples: 1000 # number of the required samples of the posterior distribution
 
 ### Datasets to consider and coefficients to fit
 The datasets and Wilson coefficients to be included in the analysis must be listed under ``datasets``
-and ``coefficients`` respectively.
+and ``coefficients`` respectively. The default order for each dataset is taken from  ``default_order``. However, it is
+possible to specify specific orders per dataset. To do this, add the key ``order`` to the dataset entry as follows.
 
 ```yaml
 datasets:
 
-  - ATLAS_tt_8TeV_ljets_Mtt
-  - ATLAS_tt_8TeV_dilep_Mtt
-  - CMS_tt_8TeV_ljets_Ytt
-  - CMS_tt2D_8TeV_dilep_MttYtt
-  - CMS_tt_13TeV_ljets_2015_Mtt
-  - CMS_tt_13TeV_dilep_2015_Mtt
-  - CMS_tt_13TeV_ljets_2016_Mtt
-  - CMS_tt_13TeV_dilep_2016_Mtt
-  - ATLAS_tt_13TeV_ljets_2016_Mtt
-  - ATLAS_CMS_tt_AC_8TeV
-  - ATLAS_tt_AC_13TeV
+  - name: ATLAS_tt_8TeV_ljets_Mtt
+  - name: ATLAS_tt_8TeV_dilep_Mtt
+    order: NLO_QCD
+  - name: CMS_tt_8TeV_ljets_Ytt
+    order: NLO_QCD
+  - name: CMS_tt2D_8TeV_dilep_MttYtt
+    order: NLO_QCD
+  - name: CMS_tt_13TeV_ljets_2015_Mtt
+    order: NLO_QCD
+  - name: CMS_tt_13TeV_dilep_2015_Mtt
+    order: NLO_QCD
+  - name: CMS_tt_13TeV_ljets_2016_Mtt
+    order: NLO_QCD
+  - name: CMS_tt_13TeV_dilep_2016_Mtt
+    order: NLO_QCD
+  - name: ATLAS_tt_13TeV_ljets_2016_Mtt
+    order: NLO_QCD
+  - name: ATLAS_CMS_tt_AC_8TeV
+    order: NLO_QCD
+  - name: ATLAS_tt_AC_13TeV
   ...
   ...
 
diff --git a/src/smefit/analyze/pca.py b/src/smefit/analyze/pca.py
index bda167c1..c10d5d8b 100644
--- a/src/smefit/analyze/pca.py
+++ b/src/smefit/analyze/pca.py
@@ -48,11 +48,11 @@ def from_dict(cls, config):
             config["data_path"],
             config["datasets"],
             config["coefficients"],
-            config["order"],
             config["use_quad"],
             config["use_theory_covmat"],
             config["use_t0"],
             config.get("use_multiplicative_prescription", False),
+            config.get("default_order", "LO"),
             config.get("theory_path", None),
             config.get("rot_to_fit_basis", None),
             config.get("uv_couplings", False),
diff --git a/src/smefit/analyze/summary.py b/src/smefit/analyze/summary.py
index 6ddefcd0..e10ec053 100644
--- a/src/smefit/analyze/summary.py
+++ b/src/smefit/analyze/summary.py
@@ -69,7 +69,6 @@ def fit_settings(self):
             summary_dict["EFT order"] = (
                 "Qudratic" if fit.config["use_quad"] else "Linear"
             )
-            summary_dict["pQCD"] = fit.config["order"]
             summary_dict["Replicas"] = fit.n_replica
             label = fit.label.replace(r"\ ", "").replace(r"\rm", "")
             summaries[label] = summary_dict
diff --git a/src/smefit/chi2.py b/src/smefit/chi2.py
index 87df6d2a..e12b3b25 100644
--- a/src/smefit/chi2.py
+++ b/src/smefit/chi2.py
@@ -85,11 +85,11 @@ def __init__(self, run_card, n_replica):
             run_card["data_path"],
             run_card["datasets"],
             run_card["coefficients"],
-            run_card["order"],
             run_card["use_quad"],
             run_card["use_theory_covmat"],
             False,
             self.use_multiplicative_prescription,
+            run_card.get("default_order", "LO"),
             run_card.get("theory_path", None),
             run_card.get("rot_to_fit_basis", None),
             run_card.get("uv_couplings", False),
diff --git a/src/smefit/fit_manager.py b/src/smefit/fit_manager.py
index 594f88d3..49853f17 100644
--- a/src/smefit/fit_manager.py
+++ b/src/smefit/fit_manager.py
@@ -117,11 +117,11 @@ def load_datasets(self):
             self.config["data_path"],
             self.config["datasets"],
             self.config["coefficients"],
-            self.config["order"],
             self.config["use_quad"],
             self.config["use_theory_covmat"],
             False,  # t0 is not used here because in the report we look at the experimental chi2
             self.config.get("use_multiplicative_prescription", False),
+            self.config.get("default_order", "LO"),
             self.config.get("theory_path", None),
             self.config.get("rot_to_fit_basis", None),
             self.config.get("uv_couplings", False),
diff --git a/src/smefit/loader.py b/src/smefit/loader.py
index 234f08b1..05b66804 100644
--- a/src/smefit/loader.py
+++ b/src/smefit/loader.py
@@ -545,11 +545,11 @@ def load_datasets(
     commondata_path,
     datasets,
     operators_to_keep,
-    order,
     use_quad,
     use_theory_covmat,
     use_t0,
     use_multiplicative_prescription,
+    default_order="LO",
     theory_path=None,
     rot_to_fit_basis=None,
     has_uv_couplings=False,
@@ -565,11 +565,11 @@ def load_datasets(
         commondata_path : str, pathlib.Path
             path to commondata folder, commondata excluded
         datasets : list
-            list of datasets to be loaded
+            List of datasets to be loaded
         operators_to_keep: list
             list of operators for which corrections are loaded
-        order: "LO", "NLO"
-            EFT perturbative order
+        default_order: str
+            Default perturbative order of the theory predictions
         use_quad: bool
             if True loads also |HO| corrections
         use_theory_covmat: bool
@@ -602,18 +602,16 @@ def load_datasets(
     th_cov = []
 
     Loader.commondata_path = pathlib.Path(commondata_path)
-    if theory_path is not None:
-        Loader.theory_path = pathlib.Path(theory_path)
-    else:
-        Loader.theory_path = pathlib.Path(commondata_path)
+    Loader.theory_path = pathlib.Path(theory_path or commondata_path)
 
     _logger.info(f"Applying cutoff scale: {cutoff_scale} GeV.")
-    for sset in np.unique(datasets):
+    for sset in datasets:
+        dataset_name = sset.get("name")
 
         dataset = Loader(
-            sset,
+            dataset_name,
             operators_to_keep,
-            order,
+            sset.get("order", default_order),
             use_quad,
             use_theory_covmat,
             use_multiplicative_prescription,
@@ -625,7 +623,7 @@ def load_datasets(
         if np.all(~dataset.mask):
             continue
 
-        exp_name.append(sset)
+        exp_name.append(dataset_name)
         n_data_exp.append(dataset.n_data)
         lumi_exp.append(dataset.lumi)
         exp_data.extend(dataset.central_values)
diff --git a/src/smefit/optimize/analytic.py b/src/smefit/optimize/analytic.py
index 064a1ae6..da6ea4b3 100644
--- a/src/smefit/optimize/analytic.py
+++ b/src/smefit/optimize/analytic.py
@@ -82,11 +82,11 @@ def from_dict(cls, config):
             config["data_path"],
             config["datasets"],
             config["coefficients"],
-            config["order"],
             False,
             config["use_theory_covmat"],
             config["use_t0"],
             False,
+            config.get("default_order", "LO"),
             config.get("theory_path", None),
             config.get("rot_to_fit_basis", None),
             config.get("uv_couplings", False),
diff --git a/src/smefit/optimize/mc.py b/src/smefit/optimize/mc.py
index 18aa4878..67a65bae 100644
--- a/src/smefit/optimize/mc.py
+++ b/src/smefit/optimize/mc.py
@@ -121,11 +121,11 @@ def from_dict(cls, config):
             config["data_path"],
             config["datasets"],
             config["coefficients"],
-            config["order"],
             config["use_quad"],
             config["use_theory_covmat"],
             config["use_t0"],
             config.get("use_multiplicative_prescription", False),
+            config.get("default_order", "LO"),
             config.get("theory_path", None),
             config.get("rot_to_fit_basis", None),
             config.get("uv_couplings", False),
diff --git a/src/smefit/optimize/ultranest.py b/src/smefit/optimize/ultranest.py
index c241d775..48757066 100644
--- a/src/smefit/optimize/ultranest.py
+++ b/src/smefit/optimize/ultranest.py
@@ -168,11 +168,11 @@ def from_dict(cls, config):
                 config["data_path"],
                 config["datasets"],
                 operators_to_keep,
-                config["order"],
                 config["use_quad"],
                 config["use_theory_covmat"],
                 config["use_t0"],
                 config.get("use_multiplicative_prescription", False),
+                config.get("default_order", "LO"),
                 config.get("theory_path", None),
                 config.get("rot_to_fit_basis", None),
                 config.get("uv_couplings", False),
diff --git a/src/smefit/prefit/__init__.py b/src/smefit/prefit/__init__.py
index 3372430c..47f067fd 100644
--- a/src/smefit/prefit/__init__.py
+++ b/src/smefit/prefit/__init__.py
@@ -19,11 +19,11 @@ def __init__(self, config):
             config["data_path"],
             config["datasets"],
             config["coefficients"],
-            config["order"],
             config["use_quad"],
             config["use_theory_covmat"],
             config["use_t0"],
             False,
+            config.get("default_order", "LO"),
             config.get("theory_path", None),
             config.get("rot_to_fit_basis", None),
             config.get("uv_couplings", False),
diff --git a/src/smefit/projections/__init__.py b/src/smefit/projections/__init__.py
index 0e10c9f1..209fb30c 100644
--- a/src/smefit/projections/__init__.py
+++ b/src/smefit/projections/__init__.py
@@ -19,10 +19,10 @@ def __init__(
         self,
         commondata_path,
         theory_path,
-        dataset_names,
+        datasets,
         projections_path,
         coefficients,
-        order,
+        default_order,
         use_quad,
         use_theory_covmat,
         rot_to_fit_basis,
@@ -32,10 +32,10 @@ def __init__(
     ):
         self.commondata_path = commondata_path
         self.theory_path = theory_path
-        self.dataset_names = dataset_names
+        self.datasets = datasets
         self.projections_path = projections_path
         self.coefficients = coefficients
-        self.order = order
+        self.default_order = default_order
         self.use_quad = use_quad
         self.use_theory_covmat = use_theory_covmat
         self.rot_to_fit_basis = rot_to_fit_basis
@@ -45,13 +45,13 @@ def __init__(
 
         self.datasets = load_datasets(
             self.commondata_path,
-            self.dataset_names,
+            self.datasets,
             self.coefficients,
-            self.order,
             self.use_quad,
             self.use_theory_covmat,
             self.use_t0,
             False,
+            self.default_order,
             theory_path=self.theory_path,
         )
 
@@ -82,10 +82,10 @@ def from_config(cls, projection_card):
         projections_path = pathlib.Path(
             projection_config["projections_path"]
         ).absolute()
-        dataset_names = projection_config["datasets"]
+        datasets = projection_config["datasets"]
 
         coefficients = projection_config.get("coefficients", [])
-        order = projection_config.get("order", "LO")
+        default_order = projection_config.get("default_order", "LO")
         use_quad = projection_config.get("use_quad", False)
         use_theory_covmat = projection_config.get("use_theory_covmat", True)
         rot_to_fit_basis = projection_config.get("rot_to_fit_basis", None)
@@ -98,10 +98,10 @@ def from_config(cls, projection_card):
         return cls(
             commondata_path,
             theory_path,
-            dataset_names,
+            datasets,
             projections_path,
             coefficients,
-            order,
+            default_order,
             use_quad,
             use_theory_covmat,
             rot_to_fit_basis,
diff --git a/src/smefit/rge.py b/src/smefit/rge.py
index ee309c57..ee9b87bd 100644
--- a/src/smefit/rge.py
+++ b/src/smefit/rge.py
@@ -379,11 +379,12 @@ def load_scales(datasets, theory_path, default_scale=1e3, cutoff_scale=None):
         list of scales for the datasets
     """
     scales = []
-    for dataset in np.unique(datasets):
+    for dataset in datasets:
+
         Loader.theory_path = pathlib.Path(theory_path)
         # dummy call just to get the scales
         _, _, _, _, dataset_scales = Loader.load_theory(
-            dataset,
+            dataset.get("name"),
             operators_to_keep={},
             order="LO",
             use_quad=False,
@@ -397,7 +398,7 @@ def load_scales(datasets, theory_path, default_scale=1e3, cutoff_scale=None):
         else:
             scales.extend([default_scale] * len(dataset_scales))
 
-        _logger.info(f"Loaded scales for dataset {dataset}: {dataset_scales}")
+        _logger.info(f"Loaded scales for dataset {dataset['name']}: {dataset_scales}")
 
     if cutoff_scale is not None:
         scales = [scale for scale in scales if scale < cutoff_scale]
diff --git a/tests/fake_results/fake_results.yaml b/tests/fake_results/fake_results.yaml
index a5c400b3..2fd1ff0e 100644
--- a/tests/fake_results/fake_results.yaml
+++ b/tests/fake_results/fake_results.yaml
@@ -10,14 +10,14 @@ data_path:  ./tests/fake_data
 theory_path:  ./tests/fake_data
 
 use_quad: False
-order: NLO
 use_theory_covmat: True
 
 
 # Datasets to include
 datasets:
 
-  - data_test5
+  - name: data_test5
+    order: "NLO"
 
 # Coefficients to fit
 coefficients:
diff --git a/tests/test_fisher.py b/tests/test_fisher.py
index dded1795..56694b15 100644
--- a/tests/test_fisher.py
+++ b/tests/test_fisher.py
@@ -71,9 +71,8 @@ def test_fisher():
         operators_to_keep = np.array(["Op1", "Op2", "Op3"])
         dataset = load_datasets(
             commondata_path,
-            datasets=["data_test5"],
+            datasets=[{"name": "data_test5", "order": "NLO"}],
             operators_to_keep=operators_to_keep,
-            order="NLO",
             use_quad=use_quad,
             use_theory_covmat=True,
             use_t0=False,
@@ -86,9 +85,8 @@ def test_fisher():
             "result_path": None,
             "result_ID": None,
             "data_path": commondata_path,
-            "datasets": ["data_test5"],
+            "datasets": [{"name": "data_test5", "order": "NLO"}],
             "coefficients": coefficients_dict,
-            "order": "NLO",
             "use_theory_covmat": True,
             "theory_path": commondata_path,
             "use_multiplicative_prescription": True,
diff --git a/tests/test_loader.py b/tests/test_loader.py
index 0e3c9149..f5d7c167 100644
--- a/tests/test_loader.py
+++ b/tests/test_loader.py
@@ -72,15 +72,17 @@
 
 def test_load_datasets():
     operators_to_keep = np.array(["Op1", "Op2", "Op4"])
-    datasets = ["data_test1", "data_test2"]
 
     for use_quad in [True, False]:
         for order in ["LO", "NLO"]:
+            datasets = [
+                {"name": "data_test1", "order": f"{order}"},
+                {"name": "data_test2", "order": f"{order}"},
+            ]
             loaded_tuple = load_datasets(
                 commondata_path,
                 datasets=datasets,
                 operators_to_keep=operators_to_keep,
-                order=order,
                 use_quad=use_quad,
                 use_theory_covmat=True,
                 use_t0=False,
@@ -117,6 +119,7 @@ def test_load_datasets():
             np.testing.assert_equal(loaded_tuple.OperatorsNames, operators_to_keep)
             np.testing.assert_equal(loaded_tuple.ExpNames, ["data_test1", "data_test2"])
             np.testing.assert_equal(loaded_tuple.NdataExp, [2, 4])
+
             np.testing.assert_equal(loaded_tuple.LinearCorrections, lin_corr)
 
             if use_quad:
@@ -152,9 +155,8 @@ def test_operator_correction_sorted():
     for order in ["LO", "NLO"]:
         loaded_tuple = load_datasets(
             commondata_path,
-            datasets=["data_test1"],
+            datasets=[{"name": "data_test1", "order": f"{order}"}],
             operators_to_keep=operators_to_keep,
-            order=order,
             use_quad=True,
             use_theory_covmat=True,
             use_t0=False,
diff --git a/tests/test_optimize.py b/tests/test_optimize.py
index a5b5f91f..b5022455 100644
--- a/tests/test_optimize.py
+++ b/tests/test_optimize.py
@@ -170,14 +170,16 @@
 path_abs = pathlib.Path(__file__).parent.resolve()
 
 
-datasets_no_corr = ["data_test1", "data_test2"]
+datasets_no_corr = [
+    {"name": "data_test1", "order": "LO"},
+    {"name": "data_test2", "order": "LO"},
+]
 config_no_corr = {}
 config_no_corr["data_path"] = commondata_path
 config_no_corr["coefficients"] = coeffs_dict
 config_no_corr["result_path"] = commondata_path
 config_no_corr["result_ID"] = "test"
 config_no_corr["datasets"] = datasets_no_corr
-config_no_corr["order"] = "LO"
 config_no_corr["use_quad"] = True
 config_no_corr["use_theory_covmat"] = True
 config_no_corr["use_t0"] = False
@@ -270,14 +272,16 @@
 chi2_corr_t0_ext = chi2_corr_t0 + chi2_ext
 
 
-datasets_corr = ["data_test3", "data_test4"]
+datasets_corr = [
+    {"name": "data_test3", "order": "LO"},
+    {"name": "data_test4", "order": "LO"},
+]
 config_corr = {}
 config_corr["data_path"] = commondata_path
 config_corr["coefficients"] = coeffs_dict
 config_corr["result_path"] = commondata_path
 config_corr["result_ID"] = "test"
 config_corr["datasets"] = datasets_corr
-config_corr["order"] = "LO"
 config_corr["use_quad"] = True
 config_corr["use_theory_covmat"] = True
 config_corr["theory_path"] = commondata_path
@@ -315,9 +319,12 @@ class TestOptimize_NS:
 
     def test_init(self):
         assert self.test_opt.results_path == commondata_path / "test"
+
         np.testing.assert_equal(
-            self.test_opt.loaded_datasets.ExpNames, datasets_no_corr
+            self.test_opt.loaded_datasets.ExpNames,
+            [dataset.get("name") for dataset in datasets_no_corr],
         )
+
         np.testing.assert_equal(
             self.test_opt.coefficients.name, ["Op1", "Op2", "Op3", "Op4"]
         )
@@ -388,7 +395,8 @@ class TestOptimize_MC:
     def test_init(self):
         assert self.test_opt.results_path == commondata_path / "test"
         np.testing.assert_equal(
-            self.test_opt.loaded_datasets.ExpNames, datasets_no_corr
+            self.test_opt.loaded_datasets.ExpNames,
+            [dataset.get("name") for dataset in datasets_no_corr],
         )
         np.testing.assert_equal(
             self.test_opt.coefficients.name, ["Op1", "Op2", "Op3", "Op4"]
@@ -487,9 +495,11 @@ def test_from_dict(self):
             opt.analytic.ALOptimizer.from_dict(config_quad)
 
     def test_init(self):
+
         assert self.test_opt.results_path == commondata_path / "test"
         np.testing.assert_equal(
-            self.test_opt.loaded_datasets.ExpNames, datasets_no_corr
+            self.test_opt.loaded_datasets.ExpNames,
+            [dataset.get("name") for dataset in datasets_no_corr],
         )
         np.testing.assert_equal(
             self.test_opt.coefficients.name, ["Op1", "Op2", "Op3", "Op4"]
diff --git a/tests/test_pca.py b/tests/test_pca.py
index af973708..169bc568 100644
--- a/tests/test_pca.py
+++ b/tests/test_pca.py
@@ -17,9 +17,8 @@
 
 dataset = load_datasets(
     commondata_path,
-    datasets=["data_test5"],
+    datasets=[{"name": "data_test5", "order": "NLO"}],
     operators_to_keep=operators_to_keep,
-    order="NLO",
     use_quad=True,
     use_theory_covmat=True,
     use_t0=False,
@@ -176,9 +175,8 @@ def test_pca_id(self):
         pca_coeffs = CoefficientManager.from_dict(pca_coeffs_dict)
         rotated_datasets = load_datasets(
             commondata_path,
-            datasets=["data_test5"],
+            datasets=[{"name": "data_test5", "order": "NLO"}],
             operators_to_keep=["PC00", "PC01", "PC02", "Op3"],
-            order="NLO",
             use_quad=True,
             use_theory_covmat=True,
             use_t0=False,