Apply rotation matrix determined by separate fixed point - python

I'm applying a rotation matrix to a group of points with the aim to align the points along the horizontal axis. Using below, the xy points I want to adjust are recorded in x and y.

This works fine but I'm hoping to transform all the points using the angle between X_Ref and Y_Ref and X_Fixed and Y_Fixed. I'm also hoping to transform the points so X_Ref and Y_Ref is at 0,0 once the rotation is completed.

The rotated points currently don't adjust for this. I'm not sure if I should account for the reference point prior to rotating or afterwards.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
import pandas as pd

df = pd.DataFrame({  
    'Period' : ['1','1','1','1'],        
    'Label' : ['A','B','C','D'],                             
    'x' : [2.0,3.0,3.0,2.0],
    'y' : [2.0,3.0,-1.0,0.0],     
    'X_Ref' : [1.0,1.0,1.0,1.0],
    'Y_Ref' : [1.0,1.0,1.0,1.0],        
    'X_Fixed' : [0.0,0.0,0.0,0.0],
    'Y_Fixed' : [0.0,0.0,0.0,0.0],                
    })

np.random.seed(1)

xy = df[['x','y']].values
Ref = df[['X_Ref','Y_Ref']].values
Fix = df[['X_Fixed','Y_Fixed']].values

fig, ax = plt.subplots()

plot_kws = {'alpha': 0.75,
        'edgecolor': 'white',
        'linewidths': 0.75}

ax.scatter(xy[:, 0], xy[:, 1], **plot_kws)
ax.scatter(Ref[:, 0], Ref[:, 1], marker = 'x')
ax.scatter(Fix[:, 0], Fix[:, 1], marker = '+')

pca = PCA(2)

# Fit the PCA object, but do not transform the data
pca.fit(xy)

# pca.components_ : array, shape (n_components, n_features)
# cos theta
ct = pca.components_[0, 0]
# sin theta
st = pca.components_[0, 1]

# One possible value of theta that lies in [0, pi]
t = np.arccos(ct)

# If t is in quadrant 1, rotate CLOCKwise by t
if ct > 0 and st > 0:
    t *= -1
# If t is in Q2, rotate COUNTERclockwise by the complement of theta
elif ct < 0 and st > 0:
    t = np.pi - t
# If t is in Q3, rotate CLOCKwise by the complement of theta
elif ct < 0 and st < 0:
    t = -(np.pi - t)
# If t is in Q4, rotate COUNTERclockwise by theta, i.e., do nothing
elif ct > 0 and st < 0:
    pass

# Manually build the ccw rotation matrix
rotmat = np.array([[np.cos(t), -np.sin(t)], 
                   [np.sin(t),  np.cos(t)]])

# Apply rotation to each row of 'm'. The output (m2)
# will be the rotated FIFA input coordinates.
m2 = (rotmat @ xy.T).T

# Center the rotated point cloud at (0, 0)
m2 -= m2.mean(axis=0)

Initial distribution:

After rotation (I have manually inserted scatter at 0,0):

Intended output:

from Apply rotation matrix determined by separate fixed point - python