diff --git a/eval.py b/eval.py
index 9d1fe64bfed8ec5988a0b0c66e15e2e22dc6b35b..6c9b88021663870cf5191b4b327a74d3c1f6d615 100755
--- a/eval.py
+++ b/eval.py
@@ -6,6 +6,7 @@ import copy
 import datetime
 import fnmatch
 import os
+import signal
 import subprocess
 import sys
 from pathlib import Path
@@ -93,6 +94,21 @@ PERF_EXE = 'perf'
 PERF_EVENT_SELECTION = '-dd'
 
 
+def perf_get_energy_events() -> str:
+    """Return the energy events supported by this system"""
+    if not perf_get_energy_events.events:  # type: ignore
+        perf_get_energy_events.events = ','.join([  # type: ignore
+            l.split()[0] for l in subprocess.check_output(
+                'perf list'.split(), text=True).splitlines() if 'energy' in l
+        ])
+    return perf_get_energy_events.events  # type: ignore
+
+
+perf_get_energy_events.events = ''  # type: ignore
+
+MEASURE_ENERGY_CMD = f'{PERF_EXE} stat -a -x; -e {{energy_events}} -o {{out}}'
+
+
 def main(args):
     """Run an evaluation"""
     for target, target_conf in TARGETS.items():
@@ -124,6 +140,13 @@ def main(args):
         if args.worker_count:
             target_env['EMPER_WORKER_COUNT'] = str(args.worker_count)
 
+        measure_energy_proc = None
+        if args.measure_energy:
+            measure_energy_out = RESULT_DIR / f'{target}.energy.stats'
+            measure_energy_cmd = MEASURE_ENERGY_CMD.format(
+                energy_events=perf_get_energy_events(), out=measure_energy_out)
+            measure_energy_proc = subprocess.Popen(measure_energy_cmd.split())
+
         out_path = RESULT_DIR / f'{target}.out'
         err_path = RESULT_DIR / f'{target}.err'
         with open(out_path, 'w', encoding='utf-8') as out_file, open(
@@ -134,6 +157,10 @@ def main(args):
                            stderr=err_file,
                            env=prepare_env(target_env))
 
+        if measure_energy_proc:
+            measure_energy_proc.send_signal(signal.SIGINT)
+            measure_energy_proc.wait()
+
         # delete empty files
         if not os.path.getsize(err_path):
             os.remove(err_path)
@@ -196,6 +223,9 @@ if __name__ == '__main__':
     parser.add_argument('--perf-record',
                         help='use perf to record a profile',
                         action='store_true')
+    parser.add_argument('--measure-energy',
+                        help='Use perf to measure the used energy',
+                        action='store_true')
     parser.add_argument('--flamegraph',
                         help='generate flamegraphs',
                         action='store_true')
diff --git a/gen_boxplots.py b/gen_boxplots.py
index 23e80838b76b706c362ff1a1549ec6343f0ffa0c..97388f68411dc93f2c94279e5fdcc3ec50b4c040 100755
--- a/gen_boxplots.py
+++ b/gen_boxplots.py
@@ -132,7 +132,8 @@ def generate_boxplot(data: Data) -> str:
 
     plots = ''
 
-    for variant, stats in data.items():
+    for variant, variant_stats in data.items():
+        stats = variant_stats['times']
         outlier_cords = ''
         for outlier in stats['outliers']:
             outlier_cords += f'(0, {outlier}) '
diff --git a/summarize.py b/summarize.py
index 010a2454fef091b1ba4233bbfa113a8220924dd9..0648b052738de913885d93f4e5a78c76dd9ef848 100755
--- a/summarize.py
+++ b/summarize.py
@@ -11,24 +11,37 @@ import yaml
 
 Variant = str
 Measurements = list[int]
-Data = dict[Variant, Measurements]
+VariantData = dict[str, T.Union[Measurements, float]]
+Data = dict[Variant, VariantData]
 
 
 def collect(result_dir) -> Data:
     """Collect the data in result_dir and return calculated averages"""
     result_dir = Path(result_dir)
 
-    data = {}
+    data: Data = {}
     for result_file_path in result_dir.iterdir():
-        if result_file_path.suffix != '.out':
-            continue
-
         variant = result_file_path.name.split('.')[0]
-        with open(result_file_path, 'r', encoding='utf-8') as result_file:
-            # record times in us
-            times = [int(line) // 1000 for line in result_file.readlines()]
+        variant_data = data.get(variant, {})
+
+        if result_file_path.suffix == '.out':
+            with open(result_file_path, 'r', encoding='utf-8') as result_file:
+                # record times in us
+                variant_data['times'] = [
+                    int(line) // 1000 for line in result_file.readlines()
+                ]
+
+            data[variant] = variant_data
 
-        data[variant] = times
+        if result_file_path.suffixes == ['.energy', '.stats']:
+            with open(result_file_path, 'r', encoding='utf-8') as result_file:
+                # record times in us
+                variant_data.update({
+                    line.split(';')[2]: float(line.split(';')[0])
+                    for line in result_file.readlines()[2:]
+                })
+
+            data[variant] = variant_data
 
     return data
 
@@ -45,54 +58,65 @@ def calc_avgs(data) -> Avgs:
 
 Outliers = list[int]
 DescriptiveStats = dict[str, T.Union[float, Outliers]]
-Stats = dict[Variant, DescriptiveStats]
+VariantStats = dict[str, DescriptiveStats]
+Stats = dict[Variant, VariantStats]
 
 
 def calc_stats(data: Data) -> Stats:
     """Calculate and return descriptive stats of all measurements in data"""
     stats = {}
-    for variant, measurements in data.items():
-        variant_stats: DescriptiveStats = {}
+    for variant, variant_data in data.items():
+        variant_stats: VariantStats = {}
         stats[variant] = variant_stats
 
-        variant_stats['mean'] = numpy.mean(measurements)
-        variant_stats['std'] = numpy.std(measurements)
-
-        measurements.sort()
-        variant_stats['min'] = measurements[0]
-        variant_stats['max'] = measurements[-1]
-        variant_stats['median'] = float(numpy.median(measurements))
-        upper_quartile = float(numpy.percentile(measurements, 75))
-        variant_stats['upper_quartile'] = upper_quartile
-        lower_quartile = float(numpy.percentile(measurements, 25))
-        variant_stats['lower_quartile'] = lower_quartile
-        iqr = upper_quartile - lower_quartile
-
-        outliers = []
-        # find whiskers
-        i = 0
-        while measurements[i] < lower_quartile - 1.5 * iqr:
-            i += 1
-        variant_stats['lower_whisker'] = measurements[i]
-        outliers = measurements[:i]
-
-        i = len(measurements) - 1
-        while measurements[i] > upper_quartile + 1.5 * iqr:
-            i -= 1
-        variant_stats['upper_whisker'] = measurements[i]
-        outliers += measurements[i + 1:]
-        variant_stats['outliers'] = outliers
-
-        # convert everything to float to easily dump it using pyyaml
-        for key, value in variant_stats.items():
-            if isinstance(value, list):
+        for key, measurements in variant_data.items():
+            key_stats: DescriptiveStats = {}
+            variant_stats[key] = key_stats
+
+            # single data point
+            if isinstance(measurements, float):
+                key_stats['mean'] = measurements
                 continue
-            variant_stats[key] = float(value)
+
+            key_stats['mean'] = numpy.mean(measurements)
+            key_stats['std'] = numpy.std(measurements)
+
+            measurements.sort()
+            key_stats['min'] = measurements[0]
+            key_stats['max'] = measurements[-1]
+            key_stats['median'] = float(numpy.median(measurements))
+            upper_quartile = float(numpy.percentile(measurements, 75))
+            key_stats['upper_quartile'] = upper_quartile
+            lower_quartile = float(numpy.percentile(measurements, 25))
+            key_stats['lower_quartile'] = lower_quartile
+            iqr = upper_quartile - lower_quartile
+
+            outliers = []
+            # find whiskers
+            i = 0
+            while measurements[i] < lower_quartile - 1.5 * iqr:
+                i += 1
+            key_stats['lower_whisker'] = measurements[i]
+            outliers = measurements[:i]
+
+            i = len(measurements) - 1
+            while measurements[i] > upper_quartile + 1.5 * iqr:
+                i -= 1
+            key_stats['upper_whisker'] = measurements[i]
+            outliers += measurements[i + 1:]
+            key_stats['outliers'] = outliers
+
+            # convert everything to float to easily dump it using pyyaml
+            for stat, value in key_stats.items():
+                if isinstance(value, list):
+                    continue
+                key_stats[stat] = float(value)
     return stats
 
 
 def summarize(avgs: T.Optional[Avgs] = None,
               stats: T.Optional[Stats] = None,
+              keys: T.Optional[T.Sequence[str]] = None,
               desc_stats: T.Optional[T.Sequence[str]] = None) -> int:
     """Print a summary for each selected key of the collected stats"""
     data = avgs or stats
@@ -106,8 +130,10 @@ def summarize(avgs: T.Optional[Avgs] = None,
             print(f'{variant}: {avgs[variant]}')
         else:
             assert stats is not None
-            for stat in desc_stats or stats[variant].keys():
-                print(f'{variant}-{stat}: {stats[variant][stat]}')
+            for key in keys or stats[variant].keys():
+                for stat in desc_stats or stats[variant][key].keys():
+                    print(
+                        f'{variant}-{key}-{stat}: {stats[variant][key][stat]}')
 
     return 0