diff --git a/mbtrack2/tracking/FeedBack.py b/mbtrack2/tracking/FeedBack.py
new file mode 100644
index 0000000000000000000000000000000000000000..14a41e0ed5fbb215ab492e80c9754d2f6e72e4e8
--- /dev/null
+++ b/mbtrack2/tracking/FeedBack.py
@@ -0,0 +1,170 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue May  3 09:36:24 2022
+
+@author: sauret
+"""
+from mbtrack2.tracking.element import Element
+import numpy as np
+
+class FeedBackTransverse(Element) : 
+    """ 
+    Transverse feedback system to counter transverse instabilities.
+    
+    Parameters
+    ----------
+    ring : Synchrotron object
+    turn_delay : int, optional
+        Number of turn delay before applying kick.
+    tap_number : int, optional
+        Number of tap for the FIR filters.
+    ampl : float, optional
+        Gain of the FIR filters.
+    phase : float, optional
+        Phase of the FIR filters in deg.
+            
+    Attributes
+    ----------
+     x_pos : array
+         Retrieves and stores the position of the bunch, measured by the BPM, 
+         in horizontal plane.
+     y_pos : array
+         Retrieves and stores the position of the bunch, measured by the BPM, 
+         in vertical plane.
+     kick_x : array
+         Stores the numerical values of the kicks for the horizontal plane.
+     kick_y : array
+         Stores the numerical values of the kicks for the vertical plane.
+     coef_x : array
+         Stores the coefficients of the FIR filter for the horizontal plane.
+     coef_y : array
+         Stores the coefficients of the FIR filter for the vertical plane. 
+    """
+    
+    def __init__(self, ring, turn_delay=1, tap_number=5, ampl =1e-1, phase = -90):
+        
+        self.ring = ring
+        self.turn_delay = turn_delay
+        self.tap_number = tap_number
+        self.ampl = ampl
+        self.phase = phase
+        self.x_pos = np.zeros((self.tap_number,))
+        self.y_pos = np.zeros((self.tap_number,))
+        self.kick_x = np.zeros((self.turn_delay+1,))
+        self.kick_y = np.zeros((self.turn_delay+1,))
+        self.coef_x = self.get_fir_x()
+        self.coef_y = self.get_fir_y()
+        
+    def get_fir_x(self):        
+        """
+        Compute the FIR coefficients for the horizontal plane.
+        
+        Returns
+        -------
+        coef_x : array
+            Array containing the FIR coefficients for the horizontal plane.
+        """
+        phi_x = 2*np.pi*self.ring.tune[0]
+        zeta = (self.phase*2*np.pi)/360
+        coef_x = np.zeros((self.tap_number,))
+        
+        it = np.zeros((self.tap_number,))
+        for k in range(self.tap_number):
+            it[k] = (-k - self.turn_delay)
+            
+        cs_x = np.zeros((self.tap_number,))
+        sn_x = np.zeros((self.tap_number,))
+        for i in range(self.tap_number):
+            cs_x[i] = (np.cos(phi_x*it[i]))
+            sn_x[i] = (np.sin(phi_x*it[i]))
+            
+        CC_x = np.zeros((5, self.tap_number))
+        for j in range(self.tap_number):
+            CC_x[0][j] = 1
+            CC_x[1][j] = cs_x[j]
+            CC_x[2][j] = sn_x[j]
+            CC_x[3][j] = it[j]*sn_x[j]
+            CC_x[4][j] = it[j]*cs_x[j]
+            
+        TCC_x = np.transpose(CC_x)        
+        W_x = np.linalg.inv(CC_x.dot(TCC_x))
+        D_x = W_x.dot(CC_x)
+       
+        for n in range(self.tap_number):
+            a_x = self.ampl*(D_x[1][n]*np.cos(zeta) + D_x[2][n])
+            coef_x[n] = (a_x)
+            
+        return coef_x
+    
+    def get_fir_y(self):
+        """
+        Compute the FIR coefficients for the vertical plane.
+        
+        Returns
+        -------
+        coef_y : array
+            Array containing the FIR coefficients for the vertical plane.
+        """
+        phi_y = 2*np.pi*self.ring.tune[1]
+        zeta = (self.phase*2*np.pi)/360
+        coef_y = np.zeros((self.tap_number,))
+        
+        it = np.zeros((self.tap_number,))
+        for k in range(0,self.tap_number):
+            it[k] = (-k - self.turn_delay)
+                
+        cs_y = np.zeros((self.tap_number,))
+        sn_y = np.zeros((self.tap_number,))
+        for i in range(self.tap_number):
+            cs_y[i] = (np.cos(phi_y*it[i]))
+            sn_y[i] = (np.sin(phi_y*it[i]))
+            
+        CC_y = np.zeros((5, self.tap_number))
+        for j in range(self.tap_number):
+            CC_y[0][j] = 1
+            CC_y[1][j] = cs_y[j]
+            CC_y[2][j] = sn_y[j]
+            CC_y[3][j] = it[j]*sn_y[j]
+            CC_y[4][j] = it[j]*cs_y[j]
+            
+        TCC_y = np.transpose(CC_y)
+        W_y = np.linalg.inv(CC_y.dot(TCC_y))        
+        D_y = W_y.dot(CC_y)
+
+        for n in range(self.tap_number):
+            a_y = self.ampl*(D_y[1][n]*np.cos(zeta) + D_y[2][n]*np.sin(zeta))
+            coef_y[n] = (a_y)    
+            
+        return coef_y
+    
+    @Element.parallel
+    def track(self, bunch):
+        """
+        Tracking method for the feedback system
+        No bunch to bunch interaction, so written for Bunch object and
+        @Element.parallel is used to handle Beam object.
+        
+        Parameters
+        ----------
+        bunch : Bunch or Beam object
+        """
+        self.x_pos[0] = bunch.mean[0] 
+        self.y_pos[0] = bunch.mean[2]
+        
+        x_kick = 0
+        y_kick = 0     
+            
+        for k in range(self.tap_number):
+            
+            x_kick += self.coef_x[k]*self.x_pos[k]
+            y_kick += self.coef_y[k]*self.y_pos[k]
+            
+        self.kick_x[-1] = x_kick
+        self.kick_y[-1] = y_kick
+       
+        bunch["xp"] -= self.kick_x[0]
+        bunch["yp"] -= self.kick_y[0]
+        self.x_pos = np.roll(self.x_pos,1)
+        self.y_pos = np.roll(self.y_pos,1)
+        self.kick_x = np.roll(self.kick_x, -1)
+        self.kick_y = np.roll(self.kick_y, -1)
\ No newline at end of file
diff --git a/mbtrack2/tracking/monitors/plotting.py b/mbtrack2/tracking/monitors/plotting.py
index e41904e462e8ccaa4b2ec5c8f7e1b5402c694e86..d7ca1cee3bdb8f64ccc001df59276aeacebcf4f7 100644
--- a/mbtrack2/tracking/monitors/plotting.py
+++ b/mbtrack2/tracking/monitors/plotting.py
@@ -358,7 +358,7 @@ def plot_phasespacedata(filename, bunch_number, x_var, y_var, turn,
              "$\\delta$"]
     
     # find the index of "turn" in time array
-    turn_index = np.where(file["PhaseSpaceData_0"]["time"][:]==turn) 
+    turn_index = np.where(file[group]["time"][:]==turn) 
     
     if len(turn_index[0]) == 0:
         raise ValueError("Turn {0} is not found. Enter turn from {1}.".