diff --git a/mbtrack2/tracking/monitors/plotting.py b/mbtrack2/tracking/monitors/plotting.py
index 5f256b91974cc7b73bc0f2a5bbe2f1593af40116..30773006a323561d397a8d6834f7f8639174c0e1 100644
--- a/mbtrack2/tracking/monitors/plotting.py
+++ b/mbtrack2/tracking/monitors/plotting.py
@@ -1107,6 +1107,8 @@ def plot_cavitydata(filename, cavity_name, phasor="cavity",
             - "mean" plots the mean phasor voltage and angle versus time for 
             non empty bunches. Needs BeamMonitor data with same save_every as 
             CavityMonitor.
+            - "tuner_diff" plots the tuner difference from optimal tuning 
+            condition (for a given bunch) and the tuning angle versus time.
     bunch_number : int, optional
         Bunch number to select. The default is 0.
     turn : int, optional
@@ -1134,7 +1136,7 @@ def plot_cavitydata(filename, cavity_name, phasor="cavity",
     units = [" voltage [MV]", " voltage [MV]", " voltage [MV]", " current [A]"]
     units_val = [1e-6, 1e-6, 1e-6, 1]
     
-    if plot_type == "bunch" or plot_type == "mean":
+    if (plot_type == "bunch") or (plot_type == "mean"):
     
         if plot_type == "bunch":
             data = [cavity_data["cavity_phasor_record"][bunch_number,:], 
@@ -1267,5 +1269,24 @@ def plot_cavitydata(filename, cavity_name, phasor="cavity",
         ax.set_ylabel("Power [kW]")
         plt.legend()
         
+    if plot_type == "tuner_diff":
+        data0 = np.angle(cavity_data["cavity_phasor_record"][bunch_number,:])
+        data1 = np.array(cavity_data["psi"])
+        data2 = np.array(cavity_data["theta_g"])
+        
+        ylabel1 = "Tuner diff. from optimal [rad]"
+        ylabel2 = r"Tuning angle $\psi$ [rad]"
+        fig, ax = plt.subplots()
+        twin = ax.twinx()
+        p1, = ax.plot(time, data0-data2+data1, color="r",label=ylabel1)
+        p2, = twin.plot(time, data1, color="b", label=ylabel2)
+        ax.set_xlabel("Turn number")
+        ax.set_ylabel(ylabel1)
+        twin.set_ylabel(ylabel2)
+        plots = [p1, p2]
+        ax.legend(handles=plots, loc="best")
+        ax.yaxis.label.set_color("r")
+        twin.yaxis.label.set_color("b")
+        
     file.close()
     return fig
diff --git a/mbtrack2/tracking/rf.py b/mbtrack2/tracking/rf.py
index 36fa232d8ded8d8b50aac97379cfab35e430d0c5..7f38c0bb0d17c35bb18ffa649a77d3f6f0be2774 100644
--- a/mbtrack2/tracking/rf.py
+++ b/mbtrack2/tracking/rf.py
@@ -993,6 +993,7 @@ class ProportionalLoop():
         Phase gain of the feedback.
     delay : int
         Feedback delay in unit of turns.
+        Must be supperior or equal to 1.
 
     """
     def __init__(self, ring, cav_res, gain_A, gain_P, delay):
@@ -1001,6 +1002,8 @@ class ProportionalLoop():
         self.gain_A = gain_A
         self.gain_P = gain_P
         self.delay = int(delay)
+        if delay < 1:
+            raise ValueError("delay must be >= 1.")
         self.volt_delay = np.ones(self.delay)*self.cav_res.Vc
         self.phase_delay = np.ones(self.delay)*self.cav_res.theta
     
@@ -1108,7 +1111,6 @@ class ProportionalIntegralLoop():
     gain : float list
         Proportional gain (Pgain) and integral gain (Igain) of the feedback 
         system in the form of a list [Pgain, Igain].
-        For IQ feedback, same gain set is used for I and Q.
         In case of normal conducting cavity (QL~1e4), the Pgain of ~1.0 and 
         Igain of ~1e4(5) are usually used.
         In case of super conducting cavity (QL > 1e6), the Pgain of ~100 
@@ -1162,9 +1164,9 @@ class ProportionalIntegralLoop():
        By using ig instead of Vg, the cavity response can be taken account.
        3) ig changes are reflected to Vg after the specifed delay (delay) of the system
 
-    Vc,IQ-->(rot)-->(-)-->(V->I,fac)-->PI-->Ig,IQ {--> Vg,IQ}
-                     |
-                     Ref
+    Vc-->(rot)-->(-)-->(V->I,fac)-->PI-->Ig --> Vg
+                  |
+                 Ref
 
     Examples
     --------
@@ -1299,7 +1301,7 @@ class ProportionalIntegralLoop():
 
 class DirectFeedback(ProportionalIntegralLoop):
     """
-    Direct RF FB (DFB) via generator current (ig) using PI controller [1].
+    Direct RF FB (DFB) via generator current (ig) using PI controller [1,2].
     
     The DirectFeedback inherits from ProportionalIntegralLoop so all 
     mandatory parameters from ProportionalIntegralLoop should be passed as
@@ -1364,6 +1366,9 @@ class DirectFeedback(ProportionalIntegralLoop):
     [1] : Akai, K. (2022). Stability analysis of rf accelerating mode with 
     feedback loops under heavy beam loading in SuperKEKB. PRAB, 25(10), 
     102002.
+    [2] : N. Yamamoto et al. (2023). Stability survey of a double RF system 
+    with RF feedback loops for bunch lengthening in a low-emittance synchrotron 
+    ring. In Proc. IPAC'23. doi:10.18429/JACoW-IPAC2023-WEPL161
 
     """
     def __init__(self, DFB_gain, DFB_phase_shift, DFB_sample_num=None,