Better plots.

pes_to_spec/model.py

@@ -336,21 +336,25 @@ class SelectRelevantLowResolution(TransformerMixin, BaseEstimator):
         sum_low_res = - np.mean(sum(list(X.values())), axis=0)
         peak_idx = self.estimate_prompt_peak(X)
         import matplotlib.pyplot as plt
-        fig = plt.figure(figsize=(8, 16))
+        fig = plt.figure(figsize=(8, 8))
         ax = plt.gca()
-        ax.plot(np.arange(peak_idx-300, peak_idx+300),
-                sum_low_res[peak_idx-300:peak_idx+300],
-                c="b",
+        first = peak_idx - 300
+        last = peak_idx + 300
+        first = max(0, first)
+        last = min(sum_low_res.shape[0]-1, last)
+        ax.plot(np.arange(first, last),
+                sum_low_res[first:last],
+                c="k",
                 label="Data")
         ax.set(title="",
-               xlabel="Photon Spectrometer channel",
-               ylabel="Sum of all Photon Spectrometer channels")
-        plt.axvline(peak_idx,
-                linewidth=3,
-                ls="--",
-                color='r',
-                label="Peak position")
-        ax.legend()
+               xlabel="Photon Spectrometer bin",
+               ylabel=r"$\sum$ PES channels")
+        #plt.axvline(peak_idx,
+        #        linewidth=3,
+        #        ls="--",
+        #        color='r',
+        #        label="Peak position")
+        #ax.legend(frameon=False)
         plt.savefig(filename)
         plt.close(fig)
 

pes_to_spec/test/offline_analysis.py

@@ -104,19 +104,27 @@ def plot_result(filename: str,
     #ax.plot(spec_raw_pe, spec_pred["deconvolved"], c='g', ls='-.', lw=3, label="Wiener filter result")
     Y = np.amax(spec_smooth)
     ax.legend(frameon=False, borderaxespad=0, loc='upper left')
-    ax.set(title=f"Beam intensity: {intensity*1e-3:.1f} mJ", #avg(stat unc) = {unc_stat}, avg(pca unc) = {unc_pca}",
+    ax.set_title(f"Beam intensity: {intensity*1e-3:.1f} mJ", loc="left")
+    ax.set(
            xlabel="Photon energy [eV]",
            ylabel="Intensity",
-           ylim=(0, 2*Y))
+           ylim=(0, 1.5*Y))
     if pes is not None:
         ax2 = plt.axes([0,0,1,1])
         # Manually set the position and relative size of the inset axes within ax1
-        ip = InsetPosition(ax, [0.65,0.6,0.35,0.4])
+        #ip = InsetPosition(ax, [0.65,0.6,0.35,0.4])
+        ip = InsetPosition(ax, [0.7,0.7,0.35,0.4])
         ax2.set_axes_locator(ip)
-        ax2.plot(pes_bin, pes[pes_to_show][pes_bin], c='black', lw=3)
+        if pes_to_show == "sum":
+            pes_plot = sum([pes[k][pes_bin] for k in pes.keys()])
+            pes_label = r"$\sum$ PES channels"
+        else:
+            pes_plot = pes[pes_to_show][pes_bin]
+            pes_label = pes_to_show
+        ax2.plot(pes_bin, pes_plot, c='black', lw=3)
         ax2.set(title=f"Low-resolution example data",
                 xlabel="Bin",
-                ylabel=f"{pes_to_show}",
+                ylabel=pes_label,
                 ylim=(0, None),
                 #labelsize=SMALL_SIZE,
                 #xticklabels=dict(fontdict=dict(fontsize=SMALL_SIZE)),
@@ -431,65 +439,67 @@ def main():
         fig.savefig(os.path.join(args.directory, "rmse.png"))
         plt.close(fig)
 
-        # SPEC integral w.r.t XGM intensity
-        fig = plt.figure(figsize=(12, 8))
-        gs = GridSpec(1, 1)
-        ax = fig.add_subplot(gs[0, 0])
-        sns.regplot(x=np.sum(spec_raw_int_t, axis=1)*de, y=xgm_flux_t[:,0], color='r', robust=True, ax=ax)
-        ax.set(title=f"",
-               xlabel="SPEC (raw) integral",
-               ylabel="XGM Intensity [uJ]",
-               )
-        fig.savefig(os.path.join(args.directory, "xgm_vs_intensity.png"))
-        plt.close(fig)
-
-        # SPEC integral w.r.t XGM intensity
-        fig = plt.figure(figsize=(12, 8))
-        gs = GridSpec(1, 1)
-        ax = fig.add_subplot(gs[0, 0])
-        sns.regplot(x=np.sum(spec_raw_int_t, axis=-1)*de, y=np.sum(spec_pred["expected"], axis=(-1, -2))*de, color='r', robust=True, ax=ax)
-        ax.set(title=f"",
-               xlabel="SPEC (raw) integral",
-               ylabel="Predicted integral",
-               )
-        fig.savefig(os.path.join(args.directory, "expected_vs_intensity.png"))
-        plt.close(fig)
-
-        fig = plt.figure(figsize=(12, 8))
-        gs = GridSpec(1, 1)
-        ax = fig.add_subplot(gs[0, 0])
-        sns.regplot(x=np.sum(spec_pred["expected"], axis=(-1, -2))*de, y=xgm_flux_t[:,0], color='r', robust=True, ax=ax)
-        ax.set(title=f"",
-               xlabel="Predicted integral",
-               ylabel="XGM intensity [uJ]",
-               )
-        fig.savefig(os.path.join(args.directory, "xgm_vs_expected.png"))
-        plt.close(fig)
+        ## SPEC integral w.r.t XGM intensity
+        #fig = plt.figure(figsize=(12, 8))
+        #gs = GridSpec(1, 1)
+        #ax = fig.add_subplot(gs[0, 0])
+        #sns.regplot(x=np.sum(spec_raw_int_t, axis=1)*de, y=xgm_flux_t[:,0], color='r', robust=True, ax=ax)
+        #ax.set(title=f"",
+        #       xlabel="SPEC (raw) integral",
+        #       ylabel="XGM Intensity [uJ]",
+        #       )
+        #fig.savefig(os.path.join(args.directory, "xgm_vs_intensity.png"))
+        #plt.close(fig)
+
+        ## SPEC integral w.r.t XGM intensity
+        #fig = plt.figure(figsize=(12, 8))
+        #gs = GridSpec(1, 1)
+        #ax = fig.add_subplot(gs[0, 0])
+        #sns.regplot(x=np.sum(spec_raw_int_t, axis=-1)*de, y=np.sum(spec_pred["expected"], axis=(-1, -2))*de, color='r', robust=True, ax=ax)
+        #ax.set(title=f"",
+        #       xlabel="SPEC (raw) integral",
+        #       ylabel="Predicted integral",
+        #       )
+        #fig.savefig(os.path.join(args.directory, "expected_vs_intensity.png"))
+        #plt.close(fig)
+
+        #fig = plt.figure(figsize=(12, 8))
+        #gs = GridSpec(1, 1)
+        #ax = fig.add_subplot(gs[0, 0])
+        #sns.regplot(x=np.sum(spec_pred["expected"], axis=(-1, -2))*de, y=xgm_flux_t[:,0], color='r', robust=True, ax=ax)
+        #ax.set(title=f"",
+        #       xlabel="Predicted integral",
+        #       ylabel="XGM intensity [uJ]",
+        #       )
+        #fig.savefig(os.path.join(args.directory, "xgm_vs_expected.png"))
+        #plt.close(fig)
 
     first, last = model.get_low_resolution_range()
-    first = max(0, first+300)
-    last = min(last-100, pes_raw_t["channel_1_D"].shape[1]-1)
-    pes_to_show = 'channel_1_D'
+    first = max(0, first+250)
+    last = min(last, pes_raw_t["channel_1_D"].shape[1]-1)
+    pes_to_show = 'sum'
     # plot
     high_int_idx = np.argsort(xgm_flux_t[:,0])
-    for idx in high_int_idx[-10:]:
-        tid = test_tids[idx]
-        plot_result(os.path.join(args.directory, f"test_{tid}.png"),
-                   {k: item[idx, 0, ...] if k != "pca"
-                       else item[0, ...]
-                       for k, item in spec_pred.items()},
-                    spec_smooth[idx, :] if showSpec else None,
-                    spec_raw_pe_t[idx, :] if showSpec else None,
-                    #spec_raw_int_t[idx, :] if showSpec else None,
-                    intensity=xgm_flux_t[idx,0],
-                    pes={k: -item[idx, :]
-                         for k, item in pes_raw_t.items()},
-                    pes_to_show=pes_to_show,
-                    pes_bin=np.arange(first, last),
-                    )
-        for ch in channels:
-            plot_pes(os.path.join(args.directory, f"test_pes_{tid}_{ch}.png"),
-                     pes_raw_t[ch][idx, first:last], first, last)
+    for q in [10, 25, 50, 75, 100]:
+        qi = int(len(high_int_idx)*(q/100.0))
+        for idx in high_int_idx[qi-10:qi]:
+            tid = test_tids[idx]
+            plot_result(os.path.join(args.directory, f"test_q{q}_{tid}.png"),
+                       {k: item[idx, 0, ...] if k != "pca"
+                           else item[0, ...]
+                           for k, item in spec_pred.items()},
+                        spec_smooth[idx, :] if showSpec else None,
+                        spec_raw_pe_t[idx, :] if showSpec else None,
+                        #spec_raw_int_t[idx, :] if showSpec else None,
+                        intensity=xgm_flux_t[idx,0],
+                        pes={k: -item[idx, :]
+                             for k, item in pes_raw_t.items()},
+                        pes_to_show=pes_to_show,
+                        pes_bin=np.arange(first, last),
+                        )
+            #for ch in channels:
+            #    plot_pes(os.path.join(args.directory, f"test_pes_{tid}_{ch}.png"),
+            #             pes_raw_t[ch][idx, first:last], first, last)
 
 if __name__ == '__main__':
     main()