diff --git a/synchrad/utils.py b/synchrad/utils.py
index 1e71202..5bbe435 100644
--- a/synchrad/utils.py
+++ b/synchrad/utils.py
@@ -5,6 +5,8 @@
 from scipy.interpolate import griddata
 from scipy.ndimage import gaussian_filter
 
+from .converters import tracksFromOPMD
+
 try:
     from tvtk.api import tvtk, write_data
     tvtk_installed = True
@@ -61,13 +63,6 @@ def get_full_spectrum(self, spect_filter=None, \
         if normalize_to_weights:
             val /= self.total_weight
 
-        if phot_num:
-            ax = self.Args['omega']
-            val /= ax[:,None,None]
-        else:
-            if lambda0_um is not None:
-                val *= J_in_um / lambda0_um
-
         return val
 
     def get_energy_spectrum(self, spect_filter=None, \
@@ -79,13 +74,13 @@ def get_energy_spectrum(self, spect_filter=None, \
         if self.Args['mode'] == 'far':
             theta_loc = 0.5 * (self.Args['theta'][1:] + self.Args['theta'][:-1])
             val_loc = 0.5 * (val[:,1:,:] + val[:,:-1,:])
-            int_theta = np.trapz( val_loc * np.sin(theta_loc)[None,:,None],
+            int_theta = np.trapezoid( val_loc * np.sin(theta_loc)[None,:,None],
                                   theta_loc, axis=1)
             val = self.Args['dph'] * int_theta.sum(-1)
 
         elif self.Args['mode'] == 'near':
             r_loc = self.Args['radius']
-            int_r = np.trapz( val * r_loc[None,:,None], r_loc, axis=1)
+            int_r = np.trapezoid( val * r_loc[None,:,None], r_loc, axis=1)
             val = self.Args['dph'] * int_r.sum(-1)
 
         return val
@@ -96,7 +91,13 @@ def get_energy(self, spect_filter=None, \
         val = self.get_energy_spectrum(spect_filter=spect_filter, \
           phot_num=phot_num, lambda0_um=lambda0_um, **kw_args)
 
-        val = np.trapz(val, self.Args['omega'])
+        if phot_num:
+            val /= self.Args['omega']
+        else:
+            if lambda0_um is not None:
+                val *= J_in_um / lambda0_um
+
+        val = np.trapezoid(val, self.Args['omega'])
         return val
 
     def get_spot(self, k0=None, spect_filter=None, \
@@ -106,13 +107,17 @@ def get_spot(self, k0=None, spect_filter=None, \
           phot_num=phot_num, lambda0_um=lambda0_um,  **kw_args)
 
         if k0 is None:
+            if phot_num:
+                ax = self.Args['omega']
+                val /= ax[:,None,None]
+            else:
+                if lambda0_um is not None:
+                    val *= J_in_um / lambda0_um
+
             if val.shape[0]>1:
-                val = np.trapz(val, self.Args['omega'], axis=0)
+                val = np.trapezoid(val, self.Args['omega'], axis=0)
             else:
                 val = val[0] * self.Args['dw']
-            #if val.shape[0]>1:
-            #    val = 0.5*(val[1:] + val[:-1])
-            #val = (val*self.Args['dw'][:, None, None]).sum(0)
         else:
             ax = self.Args['omega']
             indx = (ax<k0).sum()
diff --git a/tutorials/PIC/plot_spectrum.py b/tutorials/PIC/plot_spectrum.py
index 5fc1170..115dd41 100644
--- a/tutorials/PIC/plot_spectrum.py
+++ b/tutorials/PIC/plot_spectrum.py
@@ -53,11 +53,6 @@
 
     calc = SynchRad(file_spectrum="spectrum.h5")
 
-#    with h5py.File("spectrum.h5", "r") as f:
-#        calc.Data["radiation"] = {}
-#        for key in f["radiation"].keys():
-#            calc.Data["radiation"][key] = f[f"radiation/{key}"][...] / total_particle_weight
-
     # compute total emitted energy
     E_cutoff = 1.0  # keV
     eph_keV_m = 1.24e-9
@@ -68,21 +63,21 @@
     spect_filter = (energy_axis_full > E_cutoff)[:, np.newaxis, np.newaxis]
 
     energy_tot = calc.get_energy(
-        lambda0_um=1e6, phot_num=False, spect_filter=spect_filter
+        lambda0_um=1e6, spect_filter=spect_filter
     )
     energy_per_C = energy_tot/(e*calc.total_weight)
     print(f"Total energy emitted in >{E_cutoff:g} keV: {energy_per_C:g} J/C")
 
     # plot energy spectrum
     pyplot.figure()
-    energy_spectrum1D = calc.get_energy_spectrum(lambda0_um=1e6)
-    pyplot.semilogx(energy_axis, energy_spectrum1D)
+    energy_spectrum1D = calc.get_energy_spectrum()
+    pyplot.semilogx(energy_axis, energy_spectrum1D * 1e-3) # 1e-3 is for 0.1% units
     pyplot.xlabel("Photon energy (keV)")
     pyplot.ylabel("Brightness (ph./0.1%b.w./e$^-$)")
     pyplot.savefig("energy_spectrum_1D.png")
 
     # plot real-space emission
-    spotXY_far, ext_far = calc.get_spot_cartesian(bins=(512, 512), lambda0_um=1.0)
+    spotXY_far, ext_far = calc.get_spot_cartesian(bins=(512, 512), lambda0_um=1e6)
 
     fig = pyplot.figure(figsize=(8, 8))
     Plot2D(