diff --git a/mbtrack2/tracking/monitors/plotting.py b/mbtrack2/tracking/monitors/plotting.py
index e407151aa17b7450205afcf80f6909c56ff5f970..5cfeeba9f0386e99009cc42a8cbdef521d2af05d 100644
--- a/mbtrack2/tracking/monitors/plotting.py
+++ b/mbtrack2/tracking/monitors/plotting.py
@@ -1220,6 +1220,7 @@ def plot_cavitydata(filename,
                     bunch_number=0,
                     turn=None,
                     cm_lim=None,
+                    max = False,
                     show_objective=False):
     """
     Plot data recorded by CavityMonitor.
@@ -1273,6 +1274,10 @@ def plot_cavitydata(filename,
     cavity_data = file[cavity_name]
 
     time = np.array(cavity_data["time"])
+    if max :
+        maxT = np.max(np.nonzero(time))-1
+    else:
+        maxT = len(time)
 
     ph = {"cavity": 0, "beam": 1, "generator": 2, "ig": 3}
     labels = ["Cavity", "Beam", "Generator", "Generator"]
@@ -1307,21 +1312,21 @@ def plot_cavitydata(filename,
 
         fig, ax = plt.subplots()
         twin = ax.twinx()
-        p1, = ax.plot(time,
-                      np.abs(data[ph[phasor]]) * units_val[ph[phasor]],
+        p1, = ax.plot(time[:maxT],
+                      np.abs(data[ph[phasor]][:maxT]) * units_val[ph[phasor]],
                       color="r",
                       label=ylabel1)
-        p2, = twin.plot(time,
-                        np.angle(data[ph[phasor]]),
+        p2, = twin.plot(time[:maxT],
+                        np.angle(data[ph[phasor]])[:maxT],
                         color="b",
                         label=ylabel2)
         if show_objective:
-            o1, = ax.plot(time,
-                          np.array(cavity_data["Vc"]) * 1e-6,
+            o1, = ax.plot(time[:maxT],
+                          np.array(cavity_data["Vc"][:maxT]) * 1e-6,
                           "r--",
                           label="Cavity voltage objective value [MV]")
-            o2, = twin.plot(time,
-                            np.array(cavity_data["theta"]),
+            o2, = twin.plot(time[:maxT],
+                            np.array(cavity_data["theta"][:maxT]),
                             "b--",
                             label="Cavity phase objective value [rad]")
         ax.set_xlabel("Turn number")
diff --git a/mbtrack2/tracking/rf.py b/mbtrack2/tracking/rf.py
index 46d452f1bb610b3a95bd9147fd9763eb88bdd307..10711aa582e9260e6eff70e98e2d4c995d79ed2b 100644
--- a/mbtrack2/tracking/rf.py
+++ b/mbtrack2/tracking/rf.py
@@ -584,7 +584,7 @@ class CavityResonator():
     def DFB_ig_phasor(self):
         """Last current generator phasor of each bunch in [A]"""
         for FB in self.feedback:
-            if isinstance(FB, (ProportionalIntegralLoop, DirectFeedback)):
+            if isinstance(FB, DirectFeedback):
                 return FB.DFB_ig_phasor
         return np.zeros(self.ring.h)
 
@@ -1349,6 +1349,7 @@ class ProportionalIntegralLoop():
         self.ring = ring
         self.cav_res = cav_res
         self.Ig2Vg_mat = np.zeros((self.ring.h, self.ring.h), dtype=complex)
+        self.ig_modulation_signal = np.zeros(self.ring.h, dtype=complex)
         self.gain = gain
         self.FF = FF
 
@@ -1365,14 +1366,6 @@ class ProportionalIntegralLoop():
             raise ValueError("Bad parameter set : delay or every")
         self.sample_num = int(sample_num)
 
-        # IIR filter
-        ## IIR_cutoff ; cutoff frequecny in Hz. If 0, cutoff frequency is infinity. 
-        if IIR_cutoff == 0:
-            self.IIRcoef = 1.0
-        else:
-            self.IIRcoef = self.every * self.ring.T1 * IIR_cutoff 
-        self.IIRout = self.cav_res.Vc
-
         # init lists for FB process
         self.ig_phasor = np.ones(self.ring.h, dtype=complex) * self.Vg2Ig(
             self.cav_res.generator_phasor)
@@ -1424,6 +1417,7 @@ class ProportionalIntegralLoop():
         # update vc_previous for next turn
         self.vc_previous = self.cav_res.cavity_phasor_record[-self.sample_num:]
 
+        self.ig_phasor = self.ig_phasor + self.ig_modulation_signal
         self.ig_phasor_record = self.ig_phasor
 
         if apply_changes:
@@ -1498,7 +1492,7 @@ class ProportionalIntegralLoop():
             T = self.ring.T1 * self.every
             alpha = np.cos(omega*T)-1
             tmp = alpha * alpha - 2 * alpha
-            if alpha > 0:
+            if tmp > 0:
                 self.IIRcoef = alpha + np.sqrt(tmp)
             else:
                 self.IIRcoef = self.every * self.ring.T1 * cutoff * 2 * np.pi
@@ -1507,6 +1501,12 @@ class ProportionalIntegralLoop():
     def IIR(self,input):
         self.IIRout = (1-self.IIRcoef)*self.IIRout+self.IIRcoef*input
         return self.IIRout
+    
+    @property
+    def IIRcutoff(self):
+        """Return IIR cutoff frequency."""
+        T = self.ring.T1 * self.every
+        return 1.0/2.0/np.pi/T * np.arccos((2-2*self.IIRcoef-self.IIRcoef*self.IIRcoef)/2/(1-self.IIRcoef))
 
 class DirectFeedback(ProportionalIntegralLoop):
     """