diff --git a/irrad_control/analysis/constants.py b/irrad_control/analysis/constants.py
index d8897785..cd4daec8 100644
--- a/irrad_control/analysis/constants.py
+++ b/irrad_control/analysis/constants.py
@@ -8,3 +8,12 @@
 
 # nano prefix
 nano = 1e-9
+
+# Conversion factor for MeV/g to Mrad, 1 eV = 1.602e-19 J, 1 rad = 0.01 J/kg
+# -> MeV / g = 1e6 * 1.602e-19 J / 1e-3 kg 
+# -> MeV / g = 1e9 * 1.602e-19 J / kg
+# -> MeV / g = 1e9 * 1.602e-19 * 1e2 rad 
+# -> MeV / g = 1e11 * 1.602e-19 rad
+# -> Mev / g = 1e5 * 1.602e-19 Mrad
+# -> Mev / g = 1e5 * elementary_charge * Mrad
+MEV_PER_GRAM_TO_MRAD = 1e5 * elementary_charge
diff --git a/irrad_control/analysis/damage.py b/irrad_control/analysis/damage.py
index ab1b8635..e46cf989 100644
--- a/irrad_control/analysis/damage.py
+++ b/irrad_control/analysis/damage.py
@@ -13,7 +13,7 @@ def analyse_radiation_damage(data, config=None):
     bins = (100, 100)
 
     # Dict that holds results and error maps; bin centers
-    results = {r: None for r in ('proton', 'neq', 'tid')}
+    results = {r: None for r in ('primary', 'neq', 'tid')}
     errors = {e: None for e in results}
     bin_centers = {'x': None, 'y': None}
 
@@ -21,6 +21,7 @@ def analyse_radiation_damage(data, config=None):
     if config is None:
         
         server = None  # Only allow files with exactly one server for multipart to avoid adding unrelated fluence maps
+        ion_name = None
 
         # Loop over generator and get partial data files
         for nfile, data_part, config_part, session_basename in data:
@@ -35,6 +36,7 @@ def analyse_radiation_damage(data, config=None):
             # Only allow one fixed server for multipart
             if server is None:
                 server = server_config['name']
+                ion_name = server_config['daq']['ion']
 
             if server not in data_part:
                 raise KeyError(f"Server '{server}' not present in file {session_basename}!")
@@ -42,7 +44,7 @@ def analyse_radiation_damage(data, config=None):
             # Initialize damage and error maps
             if nfile == 0:
 
-                results['proton'], errors['proton'], bin_centers['x'], bin_centers['y'] = fluence.generate_fluence_map(beam_data=data_part[server]['Beam'],
+                results['primary'], errors['primary'], bin_centers['x'], bin_centers['y'] = fluence.generate_fluence_map(beam_data=data_part[server]['Beam'],
                                                                                                                        scan_data=data_part[server]['Scan'],
                                                                                                                        beam_sigma=beam_sigma,
                                                                                                                        bins=bins)
@@ -50,13 +52,13 @@ def analyse_radiation_damage(data, config=None):
                 if server_config['daq']['kappa'] is None:
                     del results['neq']
                 else:
-                    results['neq'] = results['proton'] * server_config['daq']['kappa']['nominal']
+                    results['neq'] = results['primary'] * server_config['daq']['kappa']['nominal']
                 
                 print(server_config['daq']['stopping_power'], type(server_config['daq']['stopping_power']))
                 if server_config['daq']['stopping_power'] is None:
                     del results['tid']
                 else:
-                    results['tid'] = formulas.tid_scan(proton_fluence=results['proton'], stopping_power=server_config['daq']['stopping_power'])
+                    results['tid'] = formulas.tid_per_scan(primary_fluence=results['primary'], stopping_power=server_config['daq']['stopping_power'])
 
                 continue
 
@@ -65,23 +67,24 @@ def analyse_radiation_damage(data, config=None):
                                                                                           beam_sigma=beam_sigma,
                                                                                           bins=bins)
             # Add to overall map
-            results['proton'] += fluence_map_part
-            errors['proton'] = (errors['proton']**2 + fluence_map_part_error**2)**.5
+            results['primary'] += fluence_map_part
+            errors['primary'] = (errors['primary']**2 + fluence_map_part_error**2)**.5
             
             # Add to eqivalent fluence map
             if 'neq' in results:
-                results['neq'] += results['proton'] * server_config['daq']['kappa']['nominal']
-                errors['neq'] = ((server_config['daq']['kappa']['nominal'] * errors['proton'])**2 + (results['proton'] * server_config['daq']['kappa']['sigma'])**2)**0.5
+                results['neq'] += results['primary'] * server_config['daq']['kappa']['nominal']
+                errors['neq'] = ((server_config['daq']['kappa']['nominal'] * errors['primary'])**2 + (results['primary'] * server_config['daq']['kappa']['sigma'])**2)**0.5
             
             if 'tid' in results:
-                results['tid'] += formulas.tid_scan(proton_fluence=results['proton'], stopping_power=server_config['daq']['stopping_power'])
-                errors['tid'] = formulas.tid_scan(proton_fluence=errors['proton'], stopping_power=server_config['daq']['stopping_power'])
+                results['tid'] += formulas.tid_per_scan(primary_fluence=results['primary'], stopping_power=server_config['daq']['stopping_power'])
+                errors['tid'] = formulas.tid_per_scan(primary_fluence=errors['primary'], stopping_power=server_config['daq']['stopping_power'])
 
     else:
 
         server = config['name']
+        ion_name = config['daq']['ion']
                     
-        results['proton'], errors['proton'], bin_centers['x'], bin_centers['y'] = fluence.generate_fluence_map(beam_data=data[server]['Beam'],
+        results['primary'], errors['primary'], bin_centers['x'], bin_centers['y'] = fluence.generate_fluence_map(beam_data=data[server]['Beam'],
                                                                                                                scan_data=data[server]['Scan'],
                                                                                                                beam_sigma=beam_sigma,
                                                                                                                bins=bins)
@@ -89,14 +92,14 @@ def analyse_radiation_damage(data, config=None):
         if config['daq']['kappa'] is None:
             del results['neq']
         else:
-            results['neq'] = results['proton'] * config['daq']['kappa']['nominal']
-            errors['neq'] = ((config['daq']['kappa']['nominal'] * errors['proton'])**2 + (results['proton'] * config['daq']['kappa']['sigma'])**2)**.5
+            results['neq'] = results['primary'] * config['daq']['kappa']['nominal']
+            errors['neq'] = ((config['daq']['kappa']['nominal'] * errors['primary'])**2 + (results['primary'] * config['daq']['kappa']['sigma'])**2)**.5
         
         if config['daq']['stopping_power'] is None:
             del results['tid']
         else:
-            results['tid'] = formulas.tid_scan(proton_fluence=results['proton'], stopping_power=config['daq']['stopping_power'])
-            errors['tid'] = formulas.tid_scan(proton_fluence=errors['proton'], stopping_power=config['daq']['stopping_power'])
+            results['tid'] = formulas.tid_per_scan(primary_fluence=results['primary'], stopping_power=config['daq']['stopping_power'])
+            errors['tid'] = formulas.tid_per_scan(primary_fluence=errors['primary'], stopping_power=config['daq']['stopping_power'])
 
     if any(a is None for a in (list(bin_centers.values()) + list(results.values()))):
         raise ValueError('Uninitialized values! Something went wrong - maybe files not found?')
@@ -123,16 +126,16 @@ def analyse_radiation_damage(data, config=None):
 
             is_dut = damage_map.shape == dut_map.shape                
 
-            fig, _ = plotting.plot_damage_map_3d(damage_map=damage_map, map_centers_x=centers_x, map_centers_y=centers_y, contour=not is_dut, damage=damage, server=server, dut=is_dut)
+            fig, _ = plotting.plot_damage_map_3d(damage_map=damage_map, map_centers_x=centers_x, map_centers_y=centers_y, contour=not is_dut, damage=damage, ion_name=ion_name, server=server, dut=is_dut)
             figs.append(fig)
 
-            fig, _ = plotting.plot_damage_error_3d(damage_map=damage_map, error_map=errors[damage] if not is_dut else dut_error_map, map_centers_x=centers_x, map_centers_y=centers_y, contour=not is_dut, damage=damage, server=server, dut=is_dut)
+            fig, _ = plotting.plot_damage_error_3d(damage_map=damage_map, error_map=errors[damage] if not is_dut else dut_error_map, map_centers_x=centers_x, map_centers_y=centers_y, contour=not is_dut, damage=damage, ion_name=ion_name,  server=server, dut=is_dut)
             figs.append(fig)
 
-            fig, _ = plotting.plot_damage_map_2d(damage_map=damage_map, map_centers_x=centers_x, map_centers_y=centers_y, damage=damage, server=server, dut=is_dut)
+            fig, _ = plotting.plot_damage_map_2d(damage_map=damage_map, map_centers_x=centers_x, map_centers_y=centers_y, damage=damage, ion_name=ion_name, server=server, dut=is_dut)
             figs.append(fig)
 
-            fig, _ = plotting.plot_damage_map_contourf(damage_map=damage_map, map_centers_x=centers_x, map_centers_y=centers_y, damage=damage, server=server, dut=is_dut)
+            fig, _ = plotting.plot_damage_map_contourf(damage_map=damage_map, map_centers_x=centers_x, map_centers_y=centers_y, damage=damage, ion_name=ion_name, server=server, dut=is_dut)
             figs.append(fig)
 
     logging.info("Finished plotting.")
diff --git a/irrad_control/analysis/dtype.py b/irrad_control/analysis/dtype.py
index 0f0fcfb7..49fdf0f9 100644
--- a/irrad_control/analysis/dtype.py
+++ b/irrad_control/analysis/dtype.py
@@ -37,8 +37,8 @@
                ('row_stop_y', '<f4'),  # # y component of the stopping position of currently-scanned row [mm]
                ('row_mean_beam_current', '<f4'),  # Mean of the beam current during scanning current row [A]
                ('row_mean_beam_current_error', '<f4'),  # Error of the beam current; quadratic addition of std of beam current and measurement error [A]
-               ('row_proton_fluence', '<f8'),  # The proton fluence during scanning current row [protons/cm^2]
-               ('row_proton_fluence_error', '<f8'),  # Error of the proton fluence during scanning current row [protons/cm^2]
+               ('row_primary_fluence', '<f8'),  # The ion fluence during scanning current row [ions/cm^2]
+               ('row_primary_fluence_error', '<f8'),  # Error of the ion fluence during scanning current row [ions/cm^2]
                ('row_tid', '<f4'),  # The TID during scanning current row [Mrad]
                ('row_tid_error', '<f4'),  # Error of the tid [Mrad]
                ('row_scan_speed', '<f4'),  # Speed with which the sample is scanned [mm/s]
@@ -48,8 +48,8 @@
 _irrad_dtype = [('timestamp', '<f8'),  # Posix-timestamp of init [s]
                 ('row_separation', '<f4'),  # Row separation e.g. step size of scan, spacing in between scanned rows [mm]
                 ('n_rows', '<i2'),  # Number of total rows in scan
-                ('aim_damage', 'S4'),  # Either NIEL or TID
-                ('aim_value', '<f4'),  # Nominal value of damage to be induced, either in neq/cm^2 or Mrad
+                ('aim_damage', 'S8'),  # Either neq, tid or primary
+                ('aim_value', '<f4'),  # Nominal value of damage to be induced, either in neq/cm^2, Mrad or ions / cm^2
                 ('min_scan_current', '<f4'),  # Minimum current for scanning [A]
                 ('scan_origin_x', '<f4'),  # x component of the scan origin from where the rel. coord. system is constructed [mm]
                 ('scan_origin_y', '<f4'),  # y component of the scan origin from where the rel. coord. system is constructed [mm]
@@ -61,18 +61,18 @@
 # Damage data dtype; contains NIEL and TID damage data on a per-scan basis
 _damage_dtype = [('timestamp', '<f8'),  # Timestamp [s]
                  ('scan', '<i2'),  # Number of *completed* scans,
-                 ('scan_proton_fluence', '<f8'),  # Proton fluence delivered in this scan [protons/cm^2]
-                 ('scan_proton_fluence_error', '<f8'),  # Error of proton fluence delivered in this scan [protons/cm^2]
+                 ('scan_primary_fluence', '<f8'),  # ion fluence delivered in this scan [ions/cm^2]
+                 ('scan_primary_fluence_error', '<f8'),  # Error of ion fluence delivered in this scan [ions/cm^2]
                  ('scan_tid', '<f8'),  # Total-ionizing dose delivered in this scan [Mrad]
                  ('scan_tid_error', '<f8')]  # Error of total-ionizing dose delivered in this scan [Mrad]
 
 
-# Result data type: contains proton as well as neutron fluence and scaling factor
+# Result data type: contains ion as well as neutron fluence and scaling factor
 _result_dtype = [('timestamp', '<f8'),
-                 ('proton_fluence', '<f8'),
-                 ('proton_fluence_error', '<f8'),
-                 ('neutron_fluence', '<f8'),
-                 ('neutron_fluence_error', '<f8'),
+                 ('primary_fluence', '<f8'),
+                 ('primary_fluence_error', '<f8'),
+                 ('neq_fluence', '<f8'),
+                 ('neq_fluence_error', '<f8'),
                  ('tid', '<f4'),
                  ('tid_error', '<f4')]
 
diff --git a/irrad_control/analysis/fluence.py b/irrad_control/analysis/fluence.py
index dac07c3a..1a52915e 100644
--- a/irrad_control/analysis/fluence.py
+++ b/irrad_control/analysis/fluence.py
@@ -88,7 +88,7 @@ def generate_fluence_map(beam_data, scan_data, beam_sigma, bins=(100, 100)):
     # Take sqrt of error map squared
     fluence_map_error = np.sqrt(fluence_map_error)                                  
 
-    # Scale from protons / mm² (intrinsic unit) to protons / cm²
+    # Scale from ions / mm² (intrinsic unit) to ions / cm²
     fluence_map *= 100
     fluence_map_error *= 100
 
@@ -403,7 +403,7 @@ def _process_row_wait(row_data, wait_beam_data, fluence_map, fluence_map_error,
     wait_mu_y = row_data['row_start_y'] - scan_y_offset
 
     # Add variation to the uncertainty
-    wait_protons_std = np.std(wait_beam_data['beam_current'])
+    wait_ions_std = np.std(wait_beam_data['beam_current'])
     
     # Loop over currents and apply Gauss kernel at given position
     for i in range(wait_beam_data.shape[0] - 1):
@@ -415,16 +415,16 @@ def _process_row_wait(row_data, wait_beam_data, fluence_map, fluence_map_error,
         # Calculate how many seconds this current was present while waiting
         wait_interval = wait_beam_data[i+1]['timestamp'] - wait_beam_data[i]['timestamp']
 
-        # Integrate over *wait_interval* to obtain number of protons induced
-        wait_protons = wait_current * wait_interval / elementary_charge
-        wait_protons_error = wait_current_error * wait_interval / elementary_charge
-        wait_protons_error = (wait_protons_error**2 + wait_protons_std**2)**.5
+        # Integrate over *wait_interval* to obtain number of ions induced
+        wait_ions = wait_current * wait_interval / elementary_charge
+        wait_ions_error = wait_current_error * wait_interval / elementary_charge
+        wait_ions_error = (wait_ions_error**2 + wait_ions_std**2)**.5
 
-        # Apply Gaussian kernel for protons
+        # Apply Gaussian kernel for ions
         apply_gauss_2d_kernel(map_2d=fluence_map,
                               map_2d_error=fluence_map_error,
-                              amplitude=wait_protons,
-                              amplitude_error=wait_protons_error,
+                              amplitude=wait_ions,
+                              amplitude_error=wait_ions_error,
                               bin_centers_x=map_bin_centers_x,
                               bin_centers_y=map_bin_centers_y,
                               mu_x=wait_mu_x,
@@ -479,23 +479,23 @@ def _process_row_scan(row_data, row_beam_data, fluence_map, fluence_map_error, r
     row_bin_center_currents = np.interp(row_bin_center_timestamps, row_beam_data['timestamp'], row_beam_data['beam_current'])
     row_bin_center_current_errors = np.interp(row_bin_center_timestamps, row_beam_data['timestamp'], row_beam_data['beam_current_error'])
 
-    # Integrate the current measurements with the times spent in each bin to calculate the amount of protons in the bin
-    row_bin_center_protons = (row_bin_center_currents * row_bin_transit_times) / elementary_charge
-    row_bin_center_proton_errors = (row_bin_center_current_errors * row_bin_transit_times) / elementary_charge
-    row_bin_center_proton_errors = (row_bin_center_proton_errors**2 + np.std(row_bin_center_protons)**2)**.5
+    # Integrate the current measurements with the times spent in each bin to calculate the amount of ions in the bin
+    row_bin_center_ions = (row_bin_center_currents * row_bin_transit_times) / elementary_charge
+    row_bin_center_ion_errors = (row_bin_center_current_errors * row_bin_transit_times) / elementary_charge
+    row_bin_center_ion_errors = (row_bin_center_ion_errors**2 + np.std(row_bin_center_ions)**2)**.5
 
     # Loop over row times
-    for i in range(row_bin_center_protons.shape[0]):
+    for i in range(row_bin_center_ions.shape[0]):
         
         # Update mean location of the distribution
         mu_x = map_bin_centers_x[(-(i+1) if row_data['row'] % 2 else i)]
         mu_y = row_data['row_start_y'] - scan_y_offset
         
-        # Apply Gaussian kernel for protons
+        # Apply Gaussian kernel for ions
         apply_gauss_2d_kernel(map_2d=fluence_map,
                               map_2d_error=fluence_map_error,
-                              amplitude=row_bin_center_protons[i],
-                              amplitude_error=row_bin_center_proton_errors[i],
+                              amplitude=row_bin_center_ions[i],
+                              amplitude_error=row_bin_center_ion_errors[i],
                               bin_centers_x=map_bin_centers_x,
                               bin_centers_y=map_bin_centers_y,
                               mu_x=mu_x,
diff --git a/irrad_control/analysis/formulas.py b/irrad_control/analysis/formulas.py
index 888586fb..fbd93468 100644
--- a/irrad_control/analysis/formulas.py
+++ b/irrad_control/analysis/formulas.py
@@ -1,51 +1,51 @@
 """Collection of analysis functions"""
-import irrad_control.analysis.constants as irrad_consts
 import numpy as np
+import irrad_control.analysis.constants as irrad_consts
 
 
-def tid_scan(proton_fluence, stopping_power):
+def tid_per_scan(primary_fluence, stopping_power):
     """
 
     Parameters
     ----------
-    proton_fluence: float
-        Number of protons per square centimeter
+    primary_fluence: float
+        Number of ions per square centimeter
 
     stopping_power:
-        Total stopping power of the protons in MeV cm^2 / g
+        Total stopping power of the ions in MeV cm^2 / g
 
     Returns
     -------
     Total ionizing dose in Mrad
     """
-    return 1e5 * irrad_consts.elementary_charge * proton_fluence * stopping_power
+    return irrad_consts.MEV_PER_GRAM_TO_MRAD * primary_fluence * stopping_power
 
 
-def tid_rate_scan(proton_flux, stopping_power):
+def tid_rate(ion_rate, stopping_power):
     """
 
     Parameters
     ----------
-    proton_flux: float
-        Number of protons per square centimeter per second
+    ion_rate: float
+        Number of ions per second
 
     stopping_power:
-        Total stopping power of the protons in MeV cm^2 / g
+        Total stopping power of the ions in MeV cm^2 / g
 
     Returns
     -------
-    Total ionizing dose in Mrad
+    Total ionizing dose rate in Mrad/s
     """
-    return 1e5 * irrad_consts.elementary_charge * proton_flux * stopping_power
+    return irrad_consts.MEV_PER_GRAM_TO_MRAD * ion_rate * stopping_power
 
 
-def proton_fluence_scan(proton_current, scan_step, scan_speed):
+def fluence_per_scan(ion_current, ion_n_charge, scan_step, scan_speed):
     """
 
     Parameters
     ----------
-    proton_current: float
-        Proton beam current in A
+    current: float
+        Ion beam current in A
 
     scan_step: float
         Separation between scanned rows in mm
@@ -55,34 +55,28 @@ def proton_fluence_scan(proton_current, scan_step, scan_speed):
 
     Returns
     -------
-    Fluence in protons / cm^2 delivered.
+    Fluence in ions / cm^2 delivered.
 
     """
-    return proton_current / (irrad_consts.elementary_charge * scan_speed * scan_step * 1e-2)
+    return ion_current / (ion_n_charge * irrad_consts.elementary_charge * scan_speed * scan_step * 1e-2)
 
 
-def proton_flux_scan(proton_current, scan_step, scan_speed, scan_duration):
+def neq_rate(ion_rate, hardness_factor):
     """
 
     Parameters
     ----------
-    proton_current: float
-        Proton beam current in A
-
-    scan_step: float
-        Separation between scanned rows in mm / cm / m
-
-    scan_speed: float
-        Speed with which rows are scanned in mm/s / cm/s / m/s
+    ion_rate: float
+        Ion rate in particles / s
 
-    scan_duration:
-        Duration which is needed to scan from first row to last row
+    hardness_factor:
+        Hardness factor to scale damage to neutron eqivalents
 
     Returns
     -------
-    Proton flux in protons/cm^2/s
+    1 MeV neutron eqivalent damages / s
     """
-    return proton_fluence_scan(proton_current, scan_step, scan_speed) / scan_duration
+    return ion_rate * hardness_factor
 
 
 def time_scan(scan_area, scan_step, scan_speed):
diff --git a/irrad_control/analysis/plotting.py b/irrad_control/analysis/plotting.py
index 90b44517..d09e9136 100644
--- a/irrad_control/analysis/plotting.py
+++ b/irrad_control/analysis/plotting.py
@@ -8,15 +8,15 @@
 from irrad_control.analysis.formulas import lin_odr
 
 
-def _get_damage_label_unit_target(damage, dut=False):
-    damage_unit = r'n$_\mathrm{eq}$ cm$^{-2}$' if damage == 'neq' else r'p cm$^{-2}$' if damage == 'proton' else 'Mrad'
-    damage_label = 'Fluence' if damage in ('neq', 'proton') else 'Total Ionizing Dose'
+def _get_damage_label_unit_target(damage, ion_name, dut=False):
+    damage_unit = r'n$_\mathrm{eq}$ cm$^{-2}$' if damage == 'neq' else f'{ion_name}s' + r' cm$^{-2}$' if damage == 'primary' else 'Mrad'
+    damage_label = 'Fluence' if damage in ('neq', 'primary') else 'Total Ionizing Dose'
     damage_target = "DUT" if dut else "Scan"
     return damage_label, damage_unit, damage_target
 
-def _apply_labels_damage_plots(ax, damage, server, cbar=None, dut=False, damage_map=None, uncertainty_map=False):
+def _apply_labels_damage_plots(ax, damage, ion_name, server, cbar=None, dut=False, damage_map=None, uncertainty_map=False):
 
-    damage_label, damage_unit, damage_target = _get_damage_label_unit_target(damage=damage, dut=dut)
+    damage_label, damage_unit, damage_target = _get_damage_label_unit_target(damage=damage, ion_name=ion_name, dut=dut)
 
     ax.set_xlabel(f'{damage_target} area horizontal / mm')
     ax.set_ylabel(f'{damage_target} area vertical / mm')
@@ -24,7 +24,7 @@ def _apply_labels_damage_plots(ax, damage, server, cbar=None, dut=False, damage_
 
     # 3D plot
     if hasattr(ax, 'set_zlabel'):
-        ax.set_zlabel(f"{damage_label} / {damage_unit}")
+        ax.set_zlabel(f"{damage_unit}")
 
     if damage_map is not None and dut and not uncertainty_map:
         mean, std = damage_map.mean(), damage_map.std()
@@ -35,9 +35,9 @@ def _apply_labels_damage_plots(ax, damage, server, cbar=None, dut=False, damage_
         cbar_label = f"{damage_label} / {damage_unit}"
         cbar.set_label(cbar_label)
 
-def _make_cbar(fig, damage_map, damage, rel_error_lims=None):
+def _make_cbar(fig, damage_map, damage, ion_name, rel_error_lims=None):
 
-    damage_label, damage_unit, _ = _get_damage_label_unit_target(damage=damage, dut=False)
+    damage_label, damage_unit, _ = _get_damage_label_unit_target(damage=damage, ion_name=ion_name, dut=False)
 
     # Make axis for cbar
     cbar_axis = plt.axes([0.85, 0.1, 0.033, 0.8])
@@ -70,7 +70,7 @@ def plot_damage_error_3d(damage_map, error_map, map_centers_x, map_centers_y, vi
     # Relative errors
     rel_damage_map = error_map / damage_map * 100.0
 
-    _make_cbar(fig=fig, damage_map=surface_3d, damage=damage_label_kwargs.get('damage', 'neq'), rel_error_lims=(rel_damage_map.min(), rel_damage_map.max()))
+    _make_cbar(fig=fig, damage_map=surface_3d, damage=damage_label_kwargs.get('damage', 'neq'), ion_name=damage_label_kwargs['ion_name'], rel_error_lims=(rel_damage_map.min(), rel_damage_map.max()))
 
     # Apply labels
     _apply_labels_damage_plots(ax=ax, damage_map=damage_map, uncertainty_map=True, **damage_label_kwargs)
@@ -99,7 +99,7 @@ def plot_damage_map_3d(damage_map, map_centers_x, map_centers_y, view_angle=(25,
     ax.view_init(*view_angle)
     ax.set_ylim(ax.get_ylim()[::-1])  # Inverty y axis in order to set origin to upper left
 
-    _make_cbar(fig=fig, damage_map=surface_3d, damage=damage_label_kwargs.get('damage', 'neq'))
+    _make_cbar(fig=fig, damage_map=surface_3d, damage=damage_label_kwargs.get('damage', 'neq'), ion_name=damage_label_kwargs['ion_name'])
     
     # Apply labels
     _apply_labels_damage_plots(ax=ax, damage_map=damage_map, **damage_label_kwargs)
diff --git a/irrad_control/gui/tabs/control_tab.py b/irrad_control/gui/tabs/control_tab.py
index fff829b0..9a746f03 100644
--- a/irrad_control/gui/tabs/control_tab.py
+++ b/irrad_control/gui/tabs/control_tab.py
@@ -45,7 +45,7 @@ def _init_tab(self, server):
 
         # Get widgets
         motorstage_widget = ic_cntrl_wdgts.MotorStageControlWidget(server=server)
-        scan_widget = ic_cntrl_wdgts.ScanControlWidget(server=server)
+        scan_widget = ic_cntrl_wdgts.ScanControlWidget(server=server, daq_setup=self.setup[server]['daq'])
         daq_widget = ic_cntrl_wdgts.DAQControlWidget(server=server, ro_device=ro_device)
         status_widget = ic_cntrl_wdgts.StatusInfoWidget('Status')
 
@@ -120,6 +120,12 @@ def check_no_beam(self, server, beam_current):
                     if time.time() - self._beam_down_timer[server] > 1.0:
                         self.send_cmd(hostname=server, target='__scan__', cmd='handle_event', cmd_data={'kwargs': {'event': 'beam_ok'}})
                         self._beam_down[server] = False
+    
+    def check_finish(self, server, eta_n_scans):
+        
+        if eta_n_scans == 0 and self.tab_widgets[server]['scan'].auto_finish_scan:
+            self.send_cmd(hostname=server, target='__scan__', cmd='handle_event', cmd_data={'kwargs': {'event': 'finish'}})
+
 
     def scan_status(self, server, status='started'):
         read_only_state = status == 'started'
diff --git a/irrad_control/gui/widgets/control_widgets.py b/irrad_control/gui/widgets/control_widgets.py
index 173d53c4..d9800d02 100644
--- a/irrad_control/gui/widgets/control_widgets.py
+++ b/irrad_control/gui/widgets/control_widgets.py
@@ -394,22 +394,24 @@ class ScanControlWidget(ControlWidget):
 
     scanParamsUpdated = QtCore.pyqtSignal(dict)
 
-    def __init__(self, server, parent=None):
+    def __init__(self, server, daq_setup, parent=None):
         super(ScanControlWidget, self).__init__(name='Scan Control', parent=parent)
 
         # Store server hostname
         self.server = server
+        self.daq_setup = daq_setup
 
         self.scan_params = {'row_sep': 1.0,
                             'scan_speed': 70.0,
                             'min_current': 0.0,
-                            'aim_damage': 'niel',
-                            'aim_value': 2e15,
+                            'aim_damage': 'primary',
+                            'aim_value': 1e14,
                             'rel_start': [0.0, 0.0],
                             'rel_end': [0.0, 0.0]}
 
         self._after_scan_container = None
         self.n_rows = None
+        self.auto_finish_scan = True
 
         self._init_ui()
 
@@ -424,6 +426,33 @@ def update_scan_params(self, **kwargs):
         self.scan_params.update(kwargs)
         self.scanParamsUpdated.emit(self.scan_params)
 
+    def _damage_toggled(self, damage_buttons, sv, se):
+        
+        # Get active button
+        active = damage_buttons.checkedButton()
+        if active is None:
+            damage = 'primary'
+        else:
+            damage = active.text().lower()
+
+        if damage == 'primary':
+            se.setSuffix(f" {self.daq_setup['ion']} / cm^2")
+            se.setRange(3, 20)
+            sv.setValue(1)
+            se.setValue(14)
+        elif damage == 'neq':
+            se.setSuffix(f" neq / cm^2")
+            se.setRange(3, 20)
+            sv.setValue(1)
+            se.setValue(14)
+        else:
+            se.setRange(1, 6)
+            se.setSuffix(' Mrad')
+            sv.setValue(1)
+            se.setValue(2)
+
+        self.update_scan_params(aim_damage=damage)
+
     def _init_ui(self):
 
         # Step size
@@ -461,21 +490,36 @@ def _init_ui(self):
         # Fluence
         label_aim_damage = QtWidgets.QLabel('Aim damage:')
         label_aim_damage.setToolTip('Select type and quantity of damage to be introduced to DUT')
-        rbtn_niel = QtWidgets.QRadioButton('NIEL')
+        rbtn_neq = QtWidgets.QRadioButton('NEQ')
         rbtn_tid = QtWidgets.QRadioButton('TID')
+        rbtn_primary = QtWidgets.QRadioButton('Primary')
+        damage_buttons = QtWidgets.QButtonGroup()
+        damage_buttons.addButton(rbtn_neq)
+        damage_buttons.addButton(rbtn_tid)
+        damage_buttons.addButton(rbtn_primary)
+        layout_aim_damage = QtWidgets.QHBoxLayout()
+
+        # Add radio buttons for different types of damage
+        if self.daq_setup['kappa'] is None:
+            damage_buttons.removeButton(rbtn_neq)    
+        if self.daq_setup['stopping_power'] is None:
+            damage_buttons.removeButton(rbtn_tid)            
+
         spx_damage_val = QtWidgets.QDoubleSpinBox()
         spx_damage_val.setRange(1e-3, 10)
         spx_damage_val.setDecimals(3)
         spx_damage_exp = QtWidgets.QSpinBox()
         spx_damage_exp.setPrefix('e ')
-        rbtn_niel.toggled.connect(lambda toggled, sv=spx_damage_val, se=spx_damage_exp:
-                                  (se.setRange(3, 20), se.setSuffix(' neq / cm^2'), sv.setValue(2), se.setValue(15)))
-        rbtn_tid.toggled.connect(lambda toggled, sv=spx_damage_val, se=spx_damage_exp:
-                                 (se.setRange(1, 6), se.setSuffix(' Mrad'), sv.setValue(1), se.setValue(3)))
-        rbtn_niel.toggled.connect(lambda toggled: self.update_scan_params(aim_damage='niel' if toggled else 'tid'))
+
+        for btn in damage_buttons.buttons():
+            btn.toggled.connect(lambda _, bg=damage_buttons, sv=spx_damage_val, se=spx_damage_exp: self._damage_toggled(bg, sv, se))
+            layout_aim_damage.addWidget(btn)
+        
         spx_damage_val.valueChanged.connect(lambda v: self.update_scan_params(aim_value=float(f'{v}e{spx_damage_exp.value()}')))
         spx_damage_exp.valueChanged.connect(lambda v: self.update_scan_params(aim_value=float(f'{spx_damage_val.value()}e{v}')))
-        rbtn_niel.toggle()
+
+        # Toggle initially
+        rbtn_primary.toggle()
 
         # Start point
         label_start = QtWidgets.QLabel('Relative start point:')
@@ -511,6 +555,12 @@ def _init_ui(self):
         spx_end_x.valueChanged.connect(lambda v: self.update_scan_params(rel_end=[v, spx_end_y.value()]))
         spx_end_y.valueChanged.connect(lambda v: self.update_scan_params(rel_start=[spx_end_x.value(), v]))
 
+        # Auto finish scan
+        checkbox_auto_finish = QtWidgets.QCheckBox('Auto finish scan')
+        checkbox_auto_finish.setToolTip("Automatically finish scan procedure when target damage is reached.")
+        checkbox_auto_finish.stateChanged.connect(lambda state: setattr(self, 'auto_finish_scan', bool(state)))
+        checkbox_auto_finish.setChecked(True)
+
         # Scan
         btn_start = QtWidgets.QPushButton('START')
         btn_start.setToolTip("Start scan.")
@@ -553,12 +603,13 @@ def _init_ui(self):
         layout_scan.addWidget(btn_pause)
         layout_scan.addWidget(btn_finish)
         layout_scan.addWidget(btn_stop)
+        layout_scan.addWidget(checkbox_auto_finish)
 
         # Add to layout
         self.add_widget(widget=[label_row_sep, spx_row_sep])
         self.add_widget(widget=[label_scan_speed, spx_scan_speed])
         self.add_widget(widget=[label_min_current, spx_min_current])
-        self.add_widget(widget=[label_aim_damage, rbtn_niel, rbtn_tid])
+        self.add_widget(widget=[label_aim_damage, layout_aim_damage])
         self.add_widget(widget=[QtWidgets.QLabel(''), spx_damage_val, spx_damage_exp])
         self.add_widget(widget=[label_start, spx_start_x, spx_start_y])
         self.add_widget(widget=[label_end, spx_end_x, spx_end_y])
diff --git a/irrad_control/ions/__init__.py b/irrad_control/ions/__init__.py
index 80c0bc48..c0ea80dc 100644
--- a/irrad_control/ions/__init__.py
+++ b/irrad_control/ions/__init__.py
@@ -4,6 +4,8 @@
 import numpy as np
 from importlib import import_module
 
+from irrad_control.analysis.constants import elementary_charge
+
 
 class IrradIon(object):
 
@@ -75,6 +77,24 @@ def _select_data(self, data_type, at_energy=None, at_index=None, as_dict=False,
 
         return _data
 
+    def rate(self, current):
+        """
+        Returns the *rate* in particles / second, calculated from *current* in Ampere.
+        For IrradIons with n_charge = 1 current / elementary charge and rate are the same
+
+        Parameters
+        ----------
+        current : float
+            Ion beam current in Ampere
+
+        Returns
+        -------
+        ion rate
+            Number of ions per second
+        """
+        # Ions per second
+        return current / (self.n_charge * elementary_charge)
+
     def ekin_range(self):
         """
         Return kinetic energy range as a tuple in MeV
@@ -114,6 +134,14 @@ def hardness_factor(self, at_energy=None, at_index=None, as_dict=False, return_i
 
 # Generate all ions
 def get_ions():
+    """
+    Returns a dict with all available IrradIon.name, IrradIon key-value pairs
+
+    Returns
+    -------
+    dict
+        dict with IrradIon.names as keys and the respective IrradIon as value
+    """
     ions = []
     for ion in os.listdir(os.path.dirname(__file__)):
         try:
diff --git a/irrad_control/ions/alpha/dut_energy.dat b/irrad_control/ions/alpha/dut_energy.dat
new file mode 100644
index 00000000..9de353f3
--- /dev/null
+++ b/irrad_control/ions/alpha/dut_energy.dat
@@ -0,0 +1,18 @@
+### This file contains beam monitor calibration data for alphas ###
+# Initial energy / MeV, energy at DUT mean / MeV,  energy at DUT sigma / MeV
+# [('initial_energy', '<f4'), ('final_energy', '<f4'), ('final_energy_std', '<f4')]
+2.800000000000000000e+01,1.801889991760253906e+01,1.835670620203018188e-01
+3.000000000000000000e+01,2.074927902221679688e+01,1.777266561985015869e-01
+3.200000000000000000e+01,2.335564804077148438e+01,1.738293319940567017e-01
+3.400000000000000000e+01,2.586944770812988281e+01,1.708793491125106812e-01
+3.600000000000000000e+01,2.831258583068847656e+01,1.686975806951522827e-01
+3.800000000000000000e+01,3.070089340209960938e+01,1.670033037662506104e-01
+4.000000000000000000e+01,3.304644393920898438e+01,1.656836122274398804e-01
+4.200000000000000000e+01,3.536841201782226562e+01,1.645111888647079468e-01
+4.400000000000000000e+01,3.764707565307617188e+01,1.635769158601760864e-01
+4.600000000000000000e+01,3.990035247802734375e+01,1.628875136375427246e-01
+4.800000000000000000e+01,4.213177871704101562e+01,1.622377187013626099e-01
+5.000000000000000000e+01,4.434376144409179688e+01,1.617864817380905151e-01
+5.200000000000000000e+01,4.653865814208984375e+01,1.613458693027496338e-01
+5.400000000000000000e+01,4.871849060058593750e+01,1.609955877065658569e-01
+5.600000000000000000e+01,5.088407897949218750e+01,1.606363803148269653e-01
diff --git a/irrad_control/ions/deuteron/dut_energy.dat b/irrad_control/ions/deuteron/dut_energy.dat
new file mode 100644
index 00000000..c7506145
--- /dev/null
+++ b/irrad_control/ions/deuteron/dut_energy.dat
@@ -0,0 +1,18 @@
+### This file contains beam monitor calibration data for deuterons ###
+# Initial energy / MeV, energy at DUT mean / MeV,  energy at DUT sigma / MeV
+# [('initial_energy', '<f4'), ('final_energy', '<f4'), ('final_energy_std', '<f4')]
+1.400000000000000000e+01,1.169845676422119141e+01,6.416837126016616821e-02
+1.500000000000000000e+01,1.283718967437744141e+01,6.333419680595397949e-02
+1.600000000000000000e+01,1.395740222930908203e+01,6.274832785129547119e-02
+1.700000000000000000e+01,1.506295490264892578e+01,6.221706792712211609e-02
+1.800000000000000000e+01,1.615665817260742188e+01,6.180695444345474243e-02
+1.900000000000000000e+01,1.724040222167968750e+01,6.100346893072128296e-02
+2.000000000000000000e+01,1.831780815124511719e+01,6.055201217532157898e-02
+2.100000000000000000e+01,1.938663101196289062e+01,6.036590412259101868e-02
+2.200000000000000000e+01,2.044962501525878906e+01,6.011645495891571045e-02
+2.300000000000000000e+01,2.150743865966796875e+01,5.995881184935569763e-02
+2.400000000000000000e+01,2.256084442138671875e+01,5.978257209062576294e-02
+2.500000000000000000e+01,2.360998725891113281e+01,5.964975431561470032e-02
+2.600000000000000000e+01,2.465555000305175781e+01,5.949187651276588440e-02
+2.700000000000000000e+01,2.569784164428710938e+01,5.942077189683914185e-02
+2.800000000000000000e+01,2.673683357238769531e+01,5.928938835859298706e-02
diff --git a/irrad_control/ions/proton/dut_energy.dat b/irrad_control/ions/proton/dut_energy.dat
index a3d96fc6..509baa5b 100644
--- a/irrad_control/ions/proton/dut_energy.dat
+++ b/irrad_control/ions/proton/dut_energy.dat
@@ -1,5 +1,21 @@
 ### This file contains beam monitor calibration data for protons ###
 # Initial energy / MeV, energy at DUT mean / MeV,  energy at DUT sigma / MeV
-7, 3, 0.5
-13, 12.1, 0.2
-14, 12.8, 0.15
\ No newline at end of file
+# [('initial_energy', '<f4'), ('final_energy', '<f4'), ('final_energy_std', '<f4')]
+7.000000000000000000e+00,4.485149383544921875e+00,7.274931669235229492e-02
+7.500000000000000000e+00,5.169362068176269531e+00,7.003546506166458130e-02
+8.000000000000000000e+00,5.822299003601074219e+00,6.818123161792755127e-02
+8.500000000000000000e+00,6.451987266540527344e+00,6.674090027809143066e-02
+9.000000000000000000e+00,7.063783645629882812e+00,6.558969616889953613e-02
+9.500000000000000000e+00,7.661623477935791016e+00,6.474756449460983276e-02
+1.000000000000000000e+01,8.248719215393066406e+00,6.371647864580154419e-02
+1.050000000000000000e+01,8.826905250549316406e+00,6.297903507947921753e-02
+1.100000000000000000e+01,9.397494316101074219e+00,6.240972504019737244e-02
+1.150000000000000000e+01,9.961688995361328125e+00,6.204782426357269287e-02
+1.200000000000000000e+01,1.052024650573730469e+01,6.163158640265464783e-02
+1.250000000000000000e+01,1.107397270202636719e+01,6.129443272948265076e-02
+1.300000000000000000e+01,1.162335300445556641e+01,6.104218214750289917e-02
+1.319999980926513672e+01,1.184196662902832031e+01,6.093817576766014099e-02
+1.339999961853027344e+01,1.205999088287353516e+01,6.086225062608718872e-02
+1.360000038146972656e+01,1.227746486663818359e+01,6.077045574784278870e-02
+1.380000019073486328e+01,1.249438858032226562e+01,6.073839962482452393e-02
+1.400000000000000000e+01,1.271081352233886719e+01,6.061572954058647156e-02
diff --git a/irrad_control/processes/converter.py b/irrad_control/processes/converter.py
index cfae02bd..20de4bdd 100644
--- a/irrad_control/processes/converter.py
+++ b/irrad_control/processes/converter.py
@@ -10,6 +10,7 @@
 import irrad_control.analysis as analysis
 import irrad_control.devices.readout as ro
 from irrad_control.processes.daq import DAQProcess
+from irrad_control.ions import get_ions
 
 
 class IrradConverter(DAQProcess):
@@ -29,6 +30,7 @@ def __init__(self, name=None):
 
         self.dtypes = analysis.dtype.IrradDtypes()
         self.hists = analysis.dtype.IrradHists()
+        self.ions = get_ions()
 
         # Call init of super class
         super(IrradConverter, self).__init__(name=name)
@@ -231,6 +233,7 @@ def _setup_daq_parameters(self, server, server_setup):
         
         daq_setup = server_setup['daq']
 
+        self._daq_params[server]['ion'] = self.ions[daq_setup['ion']]
         self._daq_params[server]['stopping_power'] = daq_setup['stopping_power'] or np.nan
         self._daq_params[server]['kappa'] = (np.nan, np.nan) if daq_setup['kappa'] is None else (daq_setup['kappa']['nominal'], daq_setup['kappa']['sigma'])
         self._daq_params[server]['lambda'] = (np.nan, np.nan) if daq_setup['lambda'] is None else (daq_setup['lambda']['nominal'], daq_setup['lambda']['sigma'])
@@ -488,7 +491,8 @@ def _interpret_beam_data(self, server, data, meta):
                                 extracted_current = self.data_arrays[server]['beam']['beam_current'][0] - blm_current
                                 logging.warning("Correcting extracted beam current from {:.2E} A to {:.2E} A".format(self.data_arrays[server]['beam']['beam_current'][0],
                                                                                                                 extracted_current))
-                                self.data_arrays[server]['beam']['beam_current'] = extracted_current
+
+                                self.data_arrays[server]['beam']['beam_current'] = beam_data['data']['current']['beam_current'] = extracted_current
 
                         except ZeroDivisionError:
                             pass
@@ -573,51 +577,52 @@ def _interpret_scan_data(self, server, data, meta):
             # Calculate mean row fluence and error
             row_mean_beam_current, row_mean_beam_current_err = self._calc_mean_and_error(data=self._scan_currents[server])
 
-            row_proton_fluence = analysis.formulas.proton_fluence_scan(proton_current=ufloat(row_mean_beam_current, row_mean_beam_current_err),
-                                                                       scan_step=self.data_arrays[server]['irrad']['row_separation'][0],
-                                                                       scan_speed=self.data_arrays[server]['scan']['row_scan_speed'][0])
+            row_primary_fluence = analysis.formulas.fluence_per_scan(ion_current=ufloat(row_mean_beam_current, row_mean_beam_current_err),
+                                                                     ion_n_charge=self._daq_params[server]['ion'].n_charge,
+                                                                     scan_step=self.data_arrays[server]['irrad']['row_separation'][0],
+                                                                     scan_speed=self.data_arrays[server]['scan']['row_scan_speed'][0])
 
-            row_proton_tid = analysis.formulas.tid_scan(proton_fluence=row_proton_fluence, stopping_power=self._daq_params[server]['stopping_power'])
+            row_tid = analysis.formulas.tid_per_scan(primary_fluence=row_primary_fluence, stopping_power=self._daq_params[server]['stopping_power'])
 
             self.data_arrays[server]['scan']['row_stop_timestamp'] = meta['timestamp']
             self.data_arrays[server]['scan']['row_stop_x'] = data['x_stop']
             self.data_arrays[server]['scan']['row_stop_y'] = data['y_stop']
             self.data_arrays[server]['scan']['row_mean_beam_current'] = row_mean_beam_current
             self.data_arrays[server]['scan']['row_mean_beam_current_error'] = row_mean_beam_current_err
-            self.data_arrays[server]['scan']['row_proton_fluence'] = row_proton_fluence.n
-            self.data_arrays[server]['scan']['row_proton_fluence_error'] = row_proton_fluence.s
-            self.data_arrays[server]['scan']['row_tid'] = row_proton_tid.n
-            self.data_arrays[server]['scan']['row_tid_error'] = row_proton_tid.s
-
-            # Log
-            logging.info("Row {} of scan {}: ({:.2E} +- {:.2E}) protons / cm^2 and ({:.2E} +- {:.2E}) Mrad".format(self.data_arrays[server]['scan']['row'][0],
-                                                                                                                   self.data_arrays[server]['scan']['scan'][0],
-                                                                                                                   row_proton_fluence.n,
-                                                                                                                   row_proton_fluence.s,
-                                                                                                                   row_proton_tid.n,
-                                                                                                                   row_proton_tid.s))
+            self.data_arrays[server]['scan']['row_primary_fluence'] = row_primary_fluence.n
+            self.data_arrays[server]['scan']['row_primary_fluence_error'] = row_primary_fluence.s
+            self.data_arrays[server]['scan']['row_tid'] = row_tid.n
+            self.data_arrays[server]['scan']['row_tid_error'] = row_tid.s
 
             # Add to overall fluence
-            self._row_fluence_hist[server][self.data_arrays[server]['scan']['row'][0]] += row_proton_fluence
+            self._row_fluence_hist[server][self.data_arrays[server]['scan']['row'][0]] += row_primary_fluence
 
             # Append data to table within this interpretation cycle
             self.data_flags[server]['scan'] = True
 
             # ETA time and n_scans
-            _mean_proton_fluence = np.mean(self._row_fluence_hist[server]).n
+            _mean_primary_fluence = np.mean(self._row_fluence_hist[server]).n
             row_scan_time = self.data_arrays[server]['scan']['row_stop_timestamp'][0] - self.data_arrays[server]['scan']['row_start_timestamp'][0]
 
             try:
                 # Check damage type
-                if self.data_arrays[server]['irrad']['aim_damage'][0] == bytes('NIEL', encoding='ascii'):
-                    # Get remaining proton fluence
-                    remainder_NIEL = self.data_arrays[server]['irrad']['aim_value'][0] / self._daq_params[server]['kappa'][0] - _mean_proton_fluence
-                    eta_n_scans = int(remainder_NIEL / row_proton_fluence.n)
+                if self.data_arrays[server]['irrad']['aim_damage'][0] == bytes('neq', encoding='ascii'):
+                    # Get remaining primary fluence
+                    aim_primary = self.data_arrays[server]['irrad']['aim_value'][0] / self._daq_params[server]['kappa'][0]
+                    remainder_primary = aim_primary - _mean_primary_fluence
+                    eta_n_scans = int(remainder_primary / row_primary_fluence.n)
+
+                elif self.data_arrays[server]['irrad']['aim_damage'][0] == bytes('tid', encoding='ascii'):
+                    remainder_tid = self.data_arrays[server]['irrad']['aim_value'][0]
+                    remainder_tid -= analysis.formulas.tid_per_scan(primary_fluence=_mean_primary_fluence,
+                                                                    stopping_power=self._daq_params[server]['stopping_power'])
+                    eta_n_scans = int(remainder_tid / analysis.formulas.tid_per_scan(primary_fluence=row_primary_fluence.n,
+                                                                                     stopping_power=self._daq_params[server]['stopping_power']))
+                # Remainder is primary fluence
                 else:
-                    remainder_TID = self.data_arrays[server]['irrad']['aim_value'][0] - analysis.formulas.tid_scan(proton_fluence=_mean_proton_fluence,
-                                                                                                                   stopping_power=self._daq_params[server]['stopping_power'])
-                    eta_n_scans = int(remainder_TID / analysis.formulas.tid_scan(proton_fluence=row_proton_fluence.n,
-                                                                                 stopping_power=self._daq_params[server]['stopping_power']))
+                    remainder_primary = self.data_arrays[server]['irrad']['aim_value'][0] - _mean_primary_fluence
+                    eta_n_scans = int(remainder_primary / row_primary_fluence.n)
+                    
 
                 eta_seconds = eta_n_scans * row_scan_time * self.data_arrays[server]['irrad']['n_rows'][0]
 
@@ -627,84 +632,85 @@ def _interpret_scan_data(self, server, data, meta):
             scan_data = {'meta': {'timestamp': meta['timestamp'], 'name': server, 'type': 'scan'},
                          'data': {'fluence_hist': unumpy.nominal_values(self._row_fluence_hist[server]).tolist(),
                                   'fluence_hist_err': unumpy.std_devs(self._row_fluence_hist[server]).tolist(),
-                                  'row_mean_proton_fluence': (row_proton_fluence.n, row_proton_fluence.s),
-                                  'row_mean_tid': (row_proton_tid.n, row_proton_tid.s),
+                                  'row_primary_fluence': (row_primary_fluence.n, row_primary_fluence.s),
+                                  'row_tid': (row_tid.n, row_tid.s),
                                   'row': int(self.data_arrays[server]['scan']['row'][0]),
                                   'eta_seconds': eta_seconds, 'eta_n_scans': eta_n_scans,
                                   'status': 'interpreted'}}
 
         elif data['status'] == 'scan_complete':
 
-            # Get scan proton fluence in each row
-            row_proton_fluences_last_scan = self.data_tables[server]['scan'].col('row_proton_fluence')[
+            # Get scan primary fluence in each row
+            row_primary_fluences_last_scan = self.data_tables[server]['scan'].col('row_primary_fluence')[
                                             -self.data_arrays[server]['irrad']['n_rows'][0]:]
 
-            # Get scan proton fluence error in each row
-            row_proton_fluences_last_scan_err = self.data_tables[server]['scan'].col('row_proton_fluence_error')[
+            # Get scan primary fluence error in each row
+            row_primary_fluences_last_scan_err = self.data_tables[server]['scan'].col('row_primary_fluence_error')[
                                                 -self.data_arrays[server]['irrad']['n_rows'][0]:]
 
-            # Calculate mean proton fluence of last scan
-            mean_scan_proton_fluence, mean_scan_proton_fluence_err = self._calc_mean_and_error(
-                data=unumpy.uarray(row_proton_fluences_last_scan,
-                                   row_proton_fluences_last_scan_err))
+            # Calculate mean primary fluence of last scan
+            mean_scan_primary_fluence, mean_scan_primary_fluence_err = self._calc_mean_and_error(
+                data=unumpy.uarray(row_primary_fluences_last_scan,
+                                   row_primary_fluences_last_scan_err))
 
             # Calculate absolute delivered fluence with this scan
-            abs_proton_fluence = ufloat(mean_scan_proton_fluence,
-                                        mean_scan_proton_fluence_err) + ufloat(
-                self.data_arrays[server]['damage']['scan_proton_fluence'][0],
-                self.data_arrays[server]['damage']['scan_proton_fluence_error'][0])
+            abs_primary_fluence = ufloat(mean_scan_primary_fluence,
+                                        mean_scan_primary_fluence_err) + ufloat(
+                self.data_arrays[server]['damage']['scan_primary_fluence'][0],
+                self.data_arrays[server]['damage']['scan_primary_fluence_error'][0])
 
             # Calculate absolute delivered TID with this scan
-            abs_tid = analysis.formulas.tid_scan(proton_fluence=abs_proton_fluence,
+            abs_tid = analysis.formulas.tid_per_scan(primary_fluence=abs_primary_fluence,
                                                  stopping_power=self._daq_params[server]['stopping_power'])
 
             # Completed scan number and timestamp of completion
             self.data_arrays[server]['damage']['timestamp'] = meta['timestamp']
             self.data_arrays[server]['damage']['scan'] = data['scan']
-            self.data_arrays[server]['damage']['scan_proton_fluence'] = abs_proton_fluence.n
-            self.data_arrays[server]['damage']['scan_proton_fluence_error'] = abs_proton_fluence.s
+            self.data_arrays[server]['damage']['scan_primary_fluence'] = abs_primary_fluence.n
+            self.data_arrays[server]['damage']['scan_primary_fluence_error'] = abs_primary_fluence.s
             self.data_arrays[server]['damage']['scan_tid'] = abs_tid.n
             self.data_arrays[server]['damage']['scan_tid_error'] = abs_tid.s
 
             # Log
             logging.info(
-                "Scan {}: ({:.2E} +- {:.2E}) protons / cm^2 and ({:.2E} +- {:.2E}) Mrad".format(data['scan'],
-                                                                                                abs_proton_fluence.n,
-                                                                                                abs_proton_fluence.s,
-                                                                                                abs_tid.n,
-                                                                                                abs_tid.s))
+                "Scan {}: ({:.2E} +- {:.2E}) {}s / cm^2 and ({:.2E} +- {:.2E}) Mrad".format(data['scan'],
+                                                                                            abs_primary_fluence.n,
+                                                                                            abs_primary_fluence.s,
+                                                                                            self._daq_params[server]['ion'].name,
+                                                                                            abs_tid.n,
+                                                                                            abs_tid.s))
 
             # Append data to table within this interpretation cycle
             self.data_flags[server]['damage'] = True
 
             scan_data = {'meta': {'timestamp': meta['timestamp'], 'name': server, 'type': 'damage'},
                          'data': {'scan': data['scan'],
-                                  'scan_proton_fluence': (abs_proton_fluence.n, abs_proton_fluence.s),
+                                  'scan_primary_fluence': (abs_primary_fluence.n, abs_primary_fluence.s),
                                   'scan_tid': (abs_tid.n, abs_tid.s)}}
 
         elif data['status'] == 'scan_finished':
 
             self.data_arrays[server]['result']['timestamp'] = meta['timestamp']
 
-            mean_result_proton_fluence = np.mean(self._row_fluence_hist[server])
-            mean_result_tid = analysis.formulas.tid_scan(proton_fluence=mean_result_proton_fluence, stopping_power=self._daq_params[server]['stopping_power'])
-            mean_result_neq_fluence = mean_result_proton_fluence * ufloat(*self._daq_params[server]['kappa'])
+            mean_result_primary_fluence = np.mean(self._row_fluence_hist[server])
+            mean_result_tid = analysis.formulas.tid_per_scan(primary_fluence=mean_result_primary_fluence, stopping_power=self._daq_params[server]['stopping_power'])
+            mean_result_neq_fluence = mean_result_primary_fluence * ufloat(*self._daq_params[server]['kappa'])
 
-            self.data_arrays[server]['result']['proton_fluence'] = mean_result_proton_fluence.n
-            self.data_arrays[server]['result']['proton_fluence_error'] = mean_result_proton_fluence.s
+            self.data_arrays[server]['result']['primary_fluence'] = mean_result_primary_fluence.n
+            self.data_arrays[server]['result']['primary_fluence_error'] = mean_result_primary_fluence.s
             self.data_arrays[server]['result']['tid'] = mean_result_tid.n
             self.data_arrays[server]['result']['tid_error'] = mean_result_tid.s
-            self.data_arrays[server]['result']['neutron_fluence'] = mean_result_neq_fluence.n
-            self.data_arrays[server]['result']['neutron_fluence_error'] = mean_result_neq_fluence.s
+            self.data_arrays[server]['result']['neq_fluence'] = mean_result_neq_fluence.n
+            self.data_arrays[server]['result']['neq_fluence_error'] = mean_result_neq_fluence.s
 
             # Append data to table within this interpretation cycle
             self.data_flags[server]['result'] = True
 
             scan_data = {'meta': {'timestamp': meta['timestamp'], 'name': server, 'type': 'result'},
                          'data': {'scan': int(self.data_arrays[server]['damage']['scan'][0]),
-                                  'proton_fluence': (mean_result_proton_fluence.n, mean_result_proton_fluence.s),
+                                  'primary_fluence': (mean_result_primary_fluence.n, mean_result_primary_fluence.s),
                                   'tid': (mean_result_tid.n, mean_result_tid.s),
-                                  'neutron_fluence': (mean_result_neq_fluence.n, mean_result_neq_fluence.s)}}
+                                  'neq_fluence': (mean_result_neq_fluence.n, mean_result_neq_fluence.s)}}
 
         return scan_data
 
diff --git a/irrad_control/processes/gui.py b/irrad_control/processes/gui.py
index 03eeb25a..4994bc81 100644
--- a/irrad_control/processes/gui.py
+++ b/irrad_control/processes/gui.py
@@ -483,7 +483,7 @@ def handle_data(self, data):
 
         elif data['meta']['type'] == 'damage':
 
-            #update_info(scan=data['data']['scan_proton_fluence'][0], unit='p/cm^2')
+            #update_info(scan=data['data']['scan_primary_fluence'][0], unit='p/cm^2')
             pass
 
         elif data['meta']['type'] == 'scan':
@@ -515,7 +515,7 @@ def handle_data(self, data):
                 # Check whether data is interpreted
             elif data['data']['status'] == 'interpreted':
                 self.monitor_tab.plots[server]['fluence_plot'].set_data(data)
-                #self.control_tab.update_info(row=data['data']['row_mean_proton_fluence'][0], unit='p/cm^2')
+                #self.control_tab.update_info(row=data['data']['row_primary_fluence'][0], unit='p/cm^2')
                 #self.control_tab.update_info(nscan=data['data']['eta_n_scans'])
 
                 if data['data']['eta_n_scans'] >= 0:
@@ -523,9 +523,8 @@ def handle_data(self, data):
                     # FIXME: more precise result would be helpful
                     pass
 
-                # Finish the scan programatically
-                if data['data']['eta_n_scans'] == 0:
-                    self.send_cmd(server, 'stage', 'finish')
+                # Finish the scan programatically, if wanted
+                self.control_tab.check_finish(server=server, eta_n_scans=data['data']['eta_n_scans'])
 
         elif data['meta']['type'] == 'temp_arduino':