<jupyter_start><jupyter_text>Basic components to understand alignment:* Linear Algebra, rotations* Linearization* Residual distribution* (Linear) Least Squares Method* Analyitcal solution for linear least squares (normal equations)* Block matrix structure (for the fast solution)* Lagrange multiplier* Track model (can be simple straight line, or more advanced e.g. GBL)* Degrees of freedom of the detector* Weak modes  Define Functions and Constants Used in the Notebook<jupyter_code>sigma_x = 0.8160  # Tobias (confirmed) number # 0.570 # TDR Number
sigma_y = 0.0163  # Tobias (confirmed) number # 0.020 # TDR Number
sigma_x = sigma_y

def sigma_phi(x, y):
    r2 = x**2 + y**2
    return ((-y / r2) ** 2 * sigma_x ** 2 + (x / r2) ** 2 * sigma_y ** 2) ** 0.5

def phi(x, y):
    return np.arctan2(y, x)

def w_phi(x, y):
    return 1 / sigma_phi(x, y) ** 2

w_x = 1 / sigma_x ** 2
w_y = 1 / sigma_y ** 2

C = np.diag([sigma_x ** 2, sigma_y ** 2])
Cinv = np.linalg.inv(C)
def space_point_prediction(z, cx, mx, cy, my):
    return np.array([
        [mx*z + cx],
        [my*z + cy],
    ])

def event_evaluation(cx, mx, cy, my, x, y, z):
    predictions = [space_point_prediction(_z, cx, mx, cy, my) for _z in z]
    delta = [np.array([[_x], [_y]]) - p for _x, _y, p in zip(x, y, predictions)]
    return sum(d.transpose() @ Cinv @ d for d in delta)[0][0]

def fit_and_calculate_residuals(df):
    df[f"fit"] = df[["x0", "x1", "x2", "y0", "y1", "y2", "z0", "z1", "z2"]].apply(
        lambda x: scipy.optimize.minimize(lambda m: event_evaluation(*m, x[0:3], x[3:6], x[6:9]), x0=[0, 0, 0 ,0]), axis=1)
    df[f"residual_x0"] = df[[f"z0", f"fit"]].apply(lambda x: space_point_prediction(x[0], *x[1].x)[0][0], axis=1) - df[f"x0"]
    df[f"residual_x1"] = df[[f"z1", f"fit"]].apply(lambda x: space_point_prediction(x[0], *x[1].x)[0][0], axis=1) - df[f"x1"]
    df[f"residual_x2"] = df[[f"z2", f"fit"]].apply(lambda x: space_point_prediction(x[0], *x[1].x)[0][0], axis=1) - df[f"x2"]
    df[f"residual_y0"] = df[[f"z0", f"fit"]].apply(lambda x: space_point_prediction(x[0], *x[1].x)[1][0], axis=1) - df[f"y0"]
    df[f"residual_y1"] = df[[f"z1", f"fit"]].apply(lambda x: space_point_prediction(x[0], *x[1].x)[1][0], axis=1) - df[f"y1"]
    df[f"residual_y2"] = df[[f"z2", f"fit"]].apply(lambda x: space_point_prediction(x[0], *x[1].x)[1][0], axis=1) - df[f"y2"]
    df[f"chi2"] = df["fit"].apply(lambda x: x.fun)
    
    df["r_x0"] = df["residual_x0"] / sigma_x
    df["r_x1"] = df["residual_x1"] / sigma_x
    df["r_x2"] = df["residual_x2"] / sigma_x
    df["r_y0"] = df["residual_y0"] / sigma_y
    df["r_y1"] = df["residual_y1"] / sigma_y
    df["r_y2"] = df["residual_y2"] / sigma_y
    
    return df


def plot_metrics(df):
    print(f'Layer 0 x Residuals: ({df[f"residual_x0"].mean() * 1000:.3f} +- {df[f"residual_x0"].std() * 1000:.3f})' + r' $\mu$m')
    print(f'Layer 1 x Residuals: ({df[f"residual_x1"].mean() * 1000:.3f} +- {df[f"residual_x1"].std() * 1000:.3f})' + r' $\mu$m')
    print(f'Layer 2 x Residuals: ({df[f"residual_x2"].mean() * 1000:.3f} +- {df[f"residual_x2"].std() * 1000:.3f})' + r' $\mu$m')
    print(f'Layer 0 y Residuals: ({df[f"residual_y0"].mean() * 1000:.3f} +- {df[f"residual_y0"].std() * 1000:.3f})' + r' $\mu$m')
    print(f'Layer 1 y Residuals: ({df[f"residual_y1"].mean() * 1000:.3f} +- {df[f"residual_y1"].std() * 1000:.3f})' + r' $\mu$m')
    print(f'Layer 2 y Residuals: ({df[f"residual_y2"].mean() * 1000:.3f} +- {df[f"residual_y2"].std() * 1000:.3f})' + r' $\mu$m')
    
    plot_residual(df[f"residual_x0"], "Layer 0", "x")
    plot_residual(df[f"residual_x1"], "Layer 1", "x")
    plot_residual(df[f"residual_x2"], "Layer 2", "x")
    plot_residual(df[f"residual_y0"], "Layer 0", "y")
    plot_residual(df[f"residual_y1"], "Layer 1", "y")
    plot_residual(df[f"residual_y2"], "Layer 2", "y")
    plot_chi2(df["chi2"])
    plt.show()
    plt.close()
def plot_residual(series, title, axis):
    if axis == "x":
        xlim = (-4*sigma_x, 4*sigma_x)
    if axis == "y":
        xlim = (-4*sigma_y, 4*sigma_y)
    _, edges, _ = plt.hist(series, bins=20, range=xlim, density=True)
    tmp = np.linspace(*xlim)
    plt.plot(tmp, scipy.stats.norm.pdf(tmp, loc=series.mean(), scale=series.std()), 
             label=r"$\mu$" + f" = {series.mean() * 1000:.3f}" + r" $\mu$m" + "\n" + r"$\sigma$" + f" = {series.std() * 1000:.3f}" + r" $\mu$m")
    plt.title(title)
    plt.xlabel(f" {axis} Residual")
    plt.ylabel("arb. units")
    plt.legend()
    plt.show()
    plt.close()
def plot_chi2(series):
    tmp = np.linspace(0, 20, num=100)
    plt.hist(series, bins=50, density=True, range=(0, 50))
    plt.plot(tmp, scipy.stats.chi2.pdf(tmp, df=2))
    plt.xlabel(r"$\chi^2$")
    plt.ylabel(r"arb. units")
    plt.show()
    plt.close()
def plot_mc_residuals(df):
    return 
    df[f"residual_cx"] = df[f"fit"].apply(lambda x: x.x[0]) - df[f"cx"]
    df[f"residual_mx"] = df[f"fit"].apply(lambda x: x.x[1]) - df[f"mx"]
    df[f"residual_cy"] = df[f"fit"].apply(lambda x: x.x[2]) - df[f"cy"]
    df[f"residual_my"] = df[f"fit"].apply(lambda x: x.x[3]) - df[f"my"]
    
    print(f'cx Residuals: ({df[f"residual_cx"].mean():.3f} +- {df[f"residual_cx"].std():.3f})')
    print(f'my Residuals: ({df[f"residual_mx"].mean():.3f} +- {df[f"residual_mx"].std():.3f})')
    print(f'cy Residuals: ({df[f"residual_cy"].mean():.3f} +- {df[f"residual_cy"].std():.3f})')
    print(f'ny Residuals: ({df[f"residual_my"].mean():.3f} +- {df[f"residual_my"].std():.3f})')
    
    plot_residual(df[f"residual_cx"], r"$c_x$ Residual", "cx")
    plot_residual(df[f"residual_mx"], r"$m_x$ Residual", "mx")
    plot_residual(df[f"residual_cy"], r"$c_y$ Residual", "mx")
    plot_residual(df[f"residual_my"], r"$m_y$ Residual", "my")
    plt.show()
    plt.close()
def module_number(eta, phi):
    return (3 - phi) + (1 - eta) * 2


def unique_plane_id(eta, phi, plane):
    return 10 * module_number(eta, phi) + plane
def read_data_and_create_df_of_spacepoints(df):
    
    df["unique_plane_id"] = df[["eta_module", "phi_module", "plane"]].apply(lambda x: unique_plane_id(*x), axis=1)
    # remove events with multiple space points in any layer
    df = df.drop_duplicates(['event', 'unique_plane_id'], keep=False)
    # TODO: Remove that, perform track fits with all combinations instead and cut by chi2

    # reshape df (sp -> event)
    mod_sp = df.pivot(index='event', columns='unique_plane_id')[['unique_plane_id', 'station', 'x','y','z']]
    mod_sp.columns = [f'{x}{y}' for x,y in mod_sp.columns]
    
    mod_sp["station"] = mod_sp["station0"]
    for i in range(3):
        del mod_sp[f"station{i}"]
        
    # remove events with a missing space points
    mod_sp = mod_sp[mod_sp[['x0', 'x1', 'x2']].isna().sum(axis=1) == 0]

    return mod_sp<jupyter_output><empty_output><jupyter_text>Read in Data Data:* r600 EHN1, station3 only* all others Ti12, all stations, lots of statistics in station 0 and 3, little in 2 Decide which kind of data to use:* ehn1, ti12: real data* simulation: simulated cosmics* toy: toy sample<jupyter_code>DATA = "toy" # cosmics, simulation, toy
mc_truth = None

STATION = 3  # 1 ... 3
MODULE = [1, 1, 1] # 0 ... 7


if DATA == "ehn1":
    input_data = "faser-data/ehn1/*.root"
    df = []
    for f in glob(input_data):
        print(f)
        try:
            df.append(read_root(f))
        except OSError:
            pass
    
    df = pd.concat(df)
    mod_sp = read_data_and_create_df_of_spacepoints(df)
    mod_sp = mod_sp.query(f"station == {STATION}")#.query(f"module0 == {MODULE[0]}").query(f"module0 == {MODULE[1]}").query(f"module0 == {MODULE[2]}")
    
if DATA == "ti12":
    input_data = "faser-data/ti12/*.root"
    df = []
    for f in glob(input_data):
        print(f)
        try:
            df.append(read_root(f))
        except OSError:
            pass
    
    df = pd.concat(df)
    mod_sp = read_data_and_create_df_of_spacepoints(df)
    
    
if DATA == "simulation":
    input_data = "faser-data/simulation/*.root"
    df = []
    for f in glob(input_data):
        try:
            df.append(read_root(f))
        except OSError:
            pass
    for i in range(len(glob(input_data))):
        df[i]["event"] += i*100000
    df = pd.concat(df)
    mod_sp = read_data_and_create_df_of_spacepoints(df)

    
if DATA == "toy":
    N = 500
    z = np.array([1, 2, 3])
    cx = scipy.stats.norm.rvs(loc=0, scale=60, size=N)
    mx = scipy.stats.norm.rvs(loc=0, scale=0.05, size=N)
    cy = scipy.stats.norm.rvs(loc=0, scale=60, size=N)
    my = scipy.stats.norm.rvs(loc=0, scale=0.05, size=N)
    
    mc_truth = pd.DataFrame({
        "cx": cx,
        "mx": mx,
        "cy": cy,
        "my": my,
    })
    
    # Nominal points on straight line
    sp0 = [space_point_prediction(z[0], _cx, _mx, _cy, _my) for _cx, _mx, _cy, _my in zip(cx, mx, cy, my)]
    sp1 = [space_point_prediction(z[1], _cx, _mx, _cy, _my) for _cx, _mx, _cy, _my in zip(cx, mx, cy, my)]
    sp2 = [space_point_prediction(z[2], _cx, _mx, _cy, _my) for _cx, _mx, _cy, _my in zip(cx, mx, cy, my)]
    
    # Resolution effect
    n0 = scipy.stats.multivariate_normal.rvs(mean=np.array([0, 0]), cov=C, size=(N))
    n1 = scipy.stats.multivariate_normal.rvs(mean=np.array([0, 0]), cov=C, size=(N))
    n2 = scipy.stats.multivariate_normal.rvs(mean=np.array([0, 0]), cov=C, size=(N))
    
    # Only add noise to x and y
    position0 = np.array([sp[:, 0] + n.transpose() for sp, n in zip(sp0, n0)])
    position1 = np.array([sp[:, 0] + n.transpose() for sp, n in zip(sp1, n1)])
    position2 = np.array([sp[:, 0] + n.transpose() for sp, n in zip(sp2, n2)])
    
    mod_sp = pd.DataFrame({
        "x0": position0[:, 0],
        "x1": position1[:, 0],
        "x2": position2[:, 0],
        "y0": position0[:, 1],
        "y1": position1[:, 1],
        "y2": position2[:, 1],
        "z0": z[0] * np.ones(N),
        "z1": z[1] * np.ones(N),
        "z2": z[2] * np.ones(N),        
    })
    
    mod_sp = pd.concat([mod_sp, mc_truth], axis=1)<jupyter_output><empty_output><jupyter_text>Uncomment this if you want to have a statistic about how many tracks go across multiple modules and which.<jupyter_code>if DATA != "toy":
    with pd.option_context('display.max_rows', None, 'display.max_columns', None):  # more options can be specified also
        print(mod_sp.groupby(["module0", "module1", "module2"]).size())<jupyter_output><empty_output><jupyter_text>Intentionally Mis-AlignRealistic scenario: Precision of installation of modules is about 0.050 millimeter.<jupyter_code>mod_sp_orig = mod_sp.copy()
x_shift = np.zeros(3)
x_shift = np.array([0.1, 0.1, 0.1])
y_shift = np.zeros(3)
y_shift = np.array([-0.06, 0.03, 0.03])
gamma = np.zeros(3)
gamma = np.array([0, 0, 0]) / 1e3 # mrad

# Rotate
mod_sp["x0"] = mod_sp[["x0", "y0"]].apply(lambda x: np.cos(gamma[0]) * x[0] - np.sin(gamma[0]) * x[1], axis=1)
mod_sp["y0"] = mod_sp[["x0", "y0"]].apply(lambda x: np.sin(gamma[0]) * x[0] + np.cos(gamma[0]) * x[1], axis=1)
mod_sp["x1"] = mod_sp[["x1", "y1"]].apply(lambda x: np.cos(gamma[1]) * x[0] - np.sin(gamma[1]) * x[1], axis=1)
mod_sp["y1"] = mod_sp[["x1", "y1"]].apply(lambda x: np.sin(gamma[1]) * x[0] + np.cos(gamma[1]) * x[1], axis=1)
mod_sp["x2"] = mod_sp[["x2", "y2"]].apply(lambda x: np.cos(gamma[2]) * x[0] - np.sin(gamma[2]) * x[1], axis=1)
mod_sp["y2"] = mod_sp[["x2", "y2"]].apply(lambda x: np.sin(gamma[2]) * x[0] + np.cos(gamma[2]) * x[1], axis=1)

# Shift
mod_sp["x0"] += x_shift[0]
mod_sp["x1"] += x_shift[1]
mod_sp["x2"] += x_shift[2]
mod_sp["y0"] += y_shift[0]
mod_sp["y1"] += y_shift[1]
mod_sp["y2"] += y_shift[2]
plt.scatter(mod_sp_orig.x1, mod_sp_orig.y1, marker='.', label="Original")
plt.scatter(mod_sp.x1, mod_sp.y1, marker='.', label="Mis-Aligned")
#plt.xlim(-20, 20)
#plt.ylim(-20, 20)
plt.legend()
plt.show()
plt.close()
z = np.array([mod_sp.z0.unique()[0], mod_sp.z1.unique()[0], mod_sp.z2.unique()[0]])
# Initial Global Shift
x_shift_global = np.mean(x_shift)
y_shift_global = np.mean(y_shift)

# Initial Global Shear
shearing = scipy.optimize.minimize(lambda m: event_evaluation(
    *m, 
    x_shift, 
    y_shift, 
    [(z[0] - z[1]), 0, (z[2]-z[1])]
), x0=[0,0,0,0])
x_shear = space_point_prediction(np.array([(z[0] - z[1]), 0, (z[2]-z[1])]), 0, shearing.x[1], 0, shearing.x[3])[0][0]
y_shear = space_point_prediction(np.array([(z[0] - z[1]), 0, (z[2]-z[1])]), 0, shearing.x[1], 0, shearing.x[3])[1][0] 
gamma_global = np.mean(gamma)
curl = scipy.optimize.minimize(lambda m: event_evaluation(
    *m, 
    [1+gamma[0], 1+gamma[1], 1+gamma[2]],
    [1-gamma[0], 1-gamma[1], 1-gamma[2]],
    [(z[0] - z[1]), 0, (z[2]-z[1])]
), x0=[0,0,0,0])
curl
gamma_curl = space_point_prediction(np.array([(z[0] - z[1]), 0, (z[2]-z[1])]), 0, curl.x[1], 0, curl.x[3])[0][0] ** 2 \
    + space_point_prediction(np.array([(z[0] - z[1]), 0, (z[2]-z[1])]), 0, curl.x[1], 0, curl.x[3])[1][0] ** 2
gamma_curl
print("Chosen mis-alignment:")
print(f"Mis-Align in x-Axis (mm)   : {x_shift[0]: .4f} \t {x_shift[1]: .4f} \t {x_shift[2]: .4f}")
print(f"Mis-Align in y-Axis (mm)   : {y_shift[0]: .4f} \t {y_shift[1]: .4f} \t {y_shift[2]: .4f}")
print(f"Mis-Align in gamma (mrad)  : {gamma[0]*1000: .4f} \t {gamma[1]*1000: .4f} \t {gamma[2]*1000: .4f}")

print("The chosen mis-alignment has the following global shifts and shears implicit:")
print(f"                              Layer 0 \t Layer 1 \t Layer 2")
print(f"Shifts in x-Axis (mm)       : {x_shift_global: .4f} \t {x_shift_global: .4f} \t {x_shift_global: .4f}")
print(f"Shifts in y-Axis (mm)       : {y_shift_global: .4f} \t {y_shift_global: .4f} \t {y_shift_global: .4f}")
print(f"Shear along x-Axis (mm)     : {x_shear[0]: .4f} \t {x_shear[1]: .4f} \t {x_shear[2]: .4f}")
print(f"Shear along y-Axis (mm)     : {y_shear[0]: .4f} \t {y_shear[1]: .4f} \t {y_shear[2]: .4f}")
print(f"Rotation along z-Axis (mrad): {gamma_global*1000: .4f} \t {gamma_global*1000: .4f} \t {gamma_global*1000: .4f}")<jupyter_output>Chosen mis-alignment:
Mis-Align in x-Axis (mm)   :  0.1000 	  0.1000 	  0.1000
Mis-Align in y-Axis (mm)   : -0.0600 	  0.0300 	  0.0300
Mis-Align in gamma (mrad)  :  0.0000 	  0.0000 	  0.0000
The chosen mis-alignment has the following global shifts and shears implicit:
                              Layer 0 	 Layer 1 	 Layer 2
Shifts in x-Axis (mm)       :  0.1000 	  0.1000 	  0.1000
Shifts in y-Axis (mm)       :  0.0000 	  0.0000 	  0.0000
Shear along x-Axis (mm)     :  0.0000 	  0.0000 	 -0.0000
Shear along y-Axis (mm)     : -0.0450 	  0.0000 	  0.0450
Rotation along z-Axis (mrad):  0.0000 	  0.0000 	  0.0000<jupyter_text>Data Exploration<jupyter_code>plot_df = mod_sp[:30]

ax = plt.axes(projection='3d')

ax.scatter3D(plot_df.z0, plot_df.x0, plot_df.y0, label="Layer 0")#, c=zdata, cmap='Greens');
ax.scatter3D(plot_df.z1, plot_df.x1, plot_df.y1, label="Layer 1")#, c=zdata, cmap='Greens');
ax.scatter3D(plot_df.z2, plot_df.x2, plot_df.y2, label="Layer 2")#, c=zdata, cmap='Greens');

ax.set_xlabel("z")
ax.set_ylabel("x")
ax.set_zlabel("y")

#ax.set_ylim(-130, 130)
#ax.set_zlim(-130, 130)

plt.legend()
plt.tight_layout()
plt.show()
plt.close()
plt.hist(mod_sp.x0, bins=50, range=(-130, 130), label="Layer 0", density=True, ls="solid", histtype="step")
plt.hist(mod_sp.x1, bins=50, range=(-130, 130), label="Layer 1", density=True, ls="dashed", histtype="step")
plt.hist(mod_sp.x2, bins=50, range=(-130, 130), label="Layer 2", density=True, ls="dotted", histtype="step")
plt.xlabel("x Position of Spacepoint in mm")
plt.ylabel("arb. units")
plt.xlim(-130, 130)
plt.legend()
plt.tight_layout()
plt.show()
plt.close()


plt.hist(mod_sp.y0, bins=50, range=(-130, 130), label="Layer 0", density=True, ls="solid", histtype="step")
plt.hist(mod_sp.y1, bins=50, range=(-130, 130), label="Layer 1", density=True, ls="dashed", histtype="step")
plt.hist(mod_sp.y2, bins=50, range=(-130, 130), label="Layer 2", density=True, ls="dotted", histtype="step")
plt.xlabel("y Position of Spacepoint in mm ")
plt.ylabel("arb. units")
plt.xlim(-130, 130)
plt.legend()
plt.tight_layout()
plt.show()
plt.close()


plt.hist(mod_sp.z0, bins=9, label="Layer 0", density=True, ls="solid", histtype="step")
plt.hist(mod_sp.z1, bins=9, label="Layer 1", density=True, ls="dashed", histtype="step")
plt.hist(mod_sp.z2, bins=9, label="Layer 2", density=True, ls="dotted", histtype="step")
plt.xlabel("z Position of Spacepoint in mm ")
plt.ylabel("arb. units")
# plt.xlim(-130, 130)
plt.legend()
plt.tight_layout()
plt.show()
plt.close()<jupyter_output><empty_output><jupyter_text>Uncomment this if you want to have a look at the initial distributions first.<jupyter_code>%%time
print(f"Number of Tracks: {len(mod_sp)}")
mod_sp = fit_and_calculate_residuals(mod_sp)
plot_metrics(mod_sp)<jupyter_output><empty_output><jupyter_text>Uncomment this if you want to select only resonable tracks.Tracks with a chi2>20 have a probability of less than 0.001%.This is only relevant when using simulation/data to remove bad tracks, and not relevant for the toy example because the toy does not have bad tracks.<jupyter_code>mod_sp = mod_sp.query("chi2 < 20")
print(f"Number of Tracks: {len(mod_sp)}")
mod_sp = fit_and_calculate_residuals(mod_sp)
plot_metrics(mod_sp)<jupyter_output><empty_output><jupyter_text>Least Square with Local and Global Parameters  Calculate to Relative Values<jupyter_code>#mod_sp = fit_and_calculate_residuals(mod_sp)
N_GLOBAL = 6
N_LOCAL = 4
N_LAGRANGE = 4
N_LAYERS = 3<jupyter_output><empty_output><jupyter_text>Local Parameters<jupyter_code>def dr_dalpha(z):
    drdalpha_x = np.array([-1, -z, 0, 0]) / sigma_x
    drdalpha_y = np.array([0, 0, -1, -z]) / sigma_y
    return np.array([
        drdalpha_x,
        drdalpha_y,
    ])
def Gamma(z):
    return sum(dr_dalpha(_z).transpose() @ dr_dalpha(_z) for _z in z)

def beta(x, y, z):
    return sum(dr_dalpha(_z).transpose() @ np.array([_x, _y]) for _x, _y, _z in zip(x, y, z))<jupyter_output><empty_output><jupyter_text>Solver. Run this if you want to have just local track fits without alignment parameters.<jupyter_code>def solve_local_parameters(df, w_y=w_y, verbose=True):
    
    beta_i = df[["r_x0", "r_x1", "r_x2", "r_y0", "r_y1", "r_y2", "z0", "z1", "z2"]].apply(lambda x: beta(x[0:3], x[3:6], x[6:]), axis=1).values
    Gamma_i = df[["z0", "z1", "z2"]].apply(lambda x: Gamma(x) , axis=1).values
    return np.array([np.linalg.inv(Gamma) @ beta for Gamma, beta in zip(Gamma_i, beta_i)])
local_solutions = solve_local_parameters(mod_sp)
print("Local Solutions (Compare numbers to toy generator settings).")
print(f"cx: {local_solutions[:, 0].mean():.4f} +- {local_solutions[:, 0].std():.4f}")
print(f"mx: {local_solutions[:, 1].mean():.4f} +- {local_solutions[:, 1].std():.4f}")
print(f"cy: {local_solutions[:, 2].mean():.4f} +- {local_solutions[:, 2].std():.4f}")
print(f"my: {local_solutions[:, 3].mean():.4f} +- {local_solutions[:, 3].std():.4f}")<jupyter_output><empty_output><jupyter_text>Global ParametersModify dr_da with the proper derivatives if you want to align something other than x and y<jupyter_code>def dr_da(x, y, z):
    drda_x = np.array([-1, 0]) / sigma_x
    drda_y = np.array([0, -1]) / sigma_y
    return np.array([
        drda_x,
        drda_y,
    ])
def constrain_x_shift():
    tmp = np.zeros(N_GLOBAL)
    i = N_GLOBAL // N_LAYERS
    tmp[0::i] = np.ones(N_LAYERS)
    return tmp

def constrain_y_shift():
    tmp = np.zeros(N_GLOBAL)
    i = N_GLOBAL // N_LAYERS
    tmp[1::i] = np.ones(N_LAYERS)
    return tmp

def constrain_x_shear(z):
    tmp = np.zeros(N_GLOBAL)
    i = N_GLOBAL // N_LAYERS
    tmp[0::i] = z
    return tmp

def constrain_y_shear(z):
    tmp = np.zeros(N_GLOBAL)
    i = N_GLOBAL // N_LAYERS
    tmp[1::i] = z
    return tmp

def constrain_gamma_rotation(x, y):
    tmp = np.zeros(N_GLOBAL)
    i = N_GLOBAL // N_LAYERS
    # tmp[2::i] = ((x*sigma_x)**2 + (y*sigma_y)**2) ** 0.5 * np.ones(N_LAYERS)
    # tmp[2::i] = np.ones(N_LAYERS)
    tmp[0::i] = np.ones(N_LAYERS)
    return tmp
    
def constrain_gamma_curl(x, y, z):
    tmp = np.zeros(N_GLOBAL)
    i = N_GLOBAL // N_LAYERS
    #tmp[2::i] = ((x*sigma_x)**2 + (y*sigma_y)**2) ** 0.5 * z
    #tmp[2::i] = z
    tmp[0::i] = ((x*sigma_x)**2 + (y*sigma_y)**2) ** 0.5 * z
    return tmp<jupyter_output><empty_output><jupyter_text>Modify the lagrange multipliers if you change what alignment correction should be applied.<jupyter_code>def C(x, y, z, N_global=N_GLOBAL, N_lagrange=N_LAGRANGE, N_layers=N_LAYERS):
    
    N = N_global + N_lagrange
    C = np.zeros([N, N])
    sub_matrix_size = N_global // N_layers
    
    for i in range(N_layers):
        C[
            i*sub_matrix_size:(i+1)*sub_matrix_size,
            i*sub_matrix_size:(i+1)*sub_matrix_size
        ] = dr_da(x[i], y[i], z[i]).transpose() @ dr_da(x[i], y[i], z[i])
        
    
    L = np.zeros([N_global, N_lagrange])
    L[:, 0] = constrain_x_shift()
    L[:, 1] = constrain_x_shear(z)
    L[:, 2] = constrain_y_shift()
    L[:, 3] = constrain_y_shear(z)
    #L[:, 0] = constrain_gamma_rotation(x, y)
    #L[:, 1] = constrain_gamma_curl(x, y, z)    
    
    C[N_global:N, :N_global] = L.transpose()
    C[:N_global, N_global:N] = L
    
    return C
def G(x, y, z, N_global=N_GLOBAL, N_lagrange=N_LAGRANGE, N_local=N_LOCAL, N_layers=N_LAYERS):
    
    G = np.zeros([N_global + N_lagrange, N_local])
    sub_matrix_size = N_global // N_layers
    
    for i in range(N_layers):
        G[i*sub_matrix_size:(i+1)*sub_matrix_size, :] = dr_da(x[i], y[i], z[i]).transpose() @ dr_dalpha(z[i])
    
    return G
def b(x, y, N_global=N_GLOBAL, N_lagrange=N_LAGRANGE, N_layers=N_LAYERS):
    b = np.zeros([N_global + N_lagrange])
    sub_matrix_size = N_global // N_layers
    
    for i in range(N_layers):
        b[i*sub_matrix_size:(i+1)*sub_matrix_size] = dr_da(x[i], y[i], z[i]).transpose() @ np.array([x[i], y[i]])
    return b<jupyter_output><empty_output><jupyter_text>Solver<jupyter_code>def solve_global_parameters(df, n_layers, w_y=w_y, verbose=True):
    
    print(f"Processing {len(df)} tracks for the alignment.")
    
    b_i = df[[
        *[f"r_x{i}" for i in range(n_layers)], 
        *[f"r_y{i}" for i in range(n_layers)]]
    ].apply(lambda x: b(x[:n_layers].values, x[n_layers:].values), axis=1).values
    
    beta_i = df[[
        *[f"r_x{i}" for i in range(n_layers)], 
        *[f"r_y{i}" for i in range(n_layers)],
        *[f"z{i}" for i in range(n_layers)], 
    ]].apply(lambda x: beta(x[:n_layers].values, x[n_layers:2*n_layers].values, x[2*n_layers:].values), axis=1).values
    
    C_i = df[[
        *[f"r_x{i}" for i in range(n_layers)], 
        *[f"r_y{i}" for i in range(n_layers)],
        *[f"z{i}" for i in range(n_layers)],
    ]].apply(lambda x: C(x[:n_layers].values, x[n_layers:2*n_layers].values, x[2*n_layers:].values), axis=1).values
        
    Gamma_i = df[[
        f"z{i}" for i in range(n_layers)
    ]].apply(lambda x: Gamma(x) , axis=1).values
        
    G_i = df[[
        *[f"r_x{i}" for i in range(n_layers)], 
        *[f"r_y{i}" for i in range(n_layers)],
        *[f"z{i}" for i in range(n_layers)],  
    ]].apply(lambda x: G(x[:n_layers].values, x[n_layers:2*n_layers].values, x[2*n_layers:].values), axis=1).values
    
    if verbose:
        for i in range(len(b_i)):
            tmp = np.linalg.inv(Gamma_i[i]) @ beta_i[i]
            print(f"Pure Local Solution [{i}]: " + f"c = {tmp[0]:.3f}, m = {tmp[1]:.3f}")
    Cprime = sum(C_i) - sum(G @ np.linalg.inv(Gamma) @ G.transpose() for G, Gamma in zip(G_i, Gamma_i))
    bprime = sum(b_i) - sum(G @ (np.linalg.inv(Gamma) @ beta) for G, Gamma, beta in zip(G_i, Gamma_i, beta_i))
    if verbose:
        print("Global Solution:")
        print(np.linalg.inv(Cprime) @ bprime)
    
    Cprime_inv = np.linalg.inv(Cprime)
    return Cprime_inv @ bprime<jupyter_output><empty_output><jupyter_text>Align all layers Pre-Alignment<jupyter_code>%%time
mod_sp = fit_and_calculate_residuals(mod_sp)
plot_metrics(mod_sp)
if mc_truth is not None:
    plot_mc_residuals(mod_sp)
%%time
result = solve_global_parameters(mod_sp, n_layers=3, verbose=False)<jupyter_output>Processing 500 tracks for the alignment.
CPU times: user 1.54 s, sys: 94.3 ms, total: 1.63 s
Wall time: 1.49 s<jupyter_text>Alignment result.The x resolution is a factor of 50 worse than the y resolution and therefore needs a factor of 250 more statistics to get the same precision. Values differing from the initial setting for the x-Axis shifts are expected when the sample is small.<jupyter_code># Use this if only xy alignment
z = np.array([mod_sp.z0.unique()[0], mod_sp.z1.unique()[0], mod_sp.z2.unique()[0]])

print("Alignment Constants")
print(f"                       Layer 0 \t Layer 1 \t Layer 2 \t Global \t Shear ")
print(f"Shifts in x-Axis (mm): {result[0]: .4f} \t {result[2]: .4f} \t {result[4]: .4f} \t {sum(result[0:6:2]): .4f} \t {sum(result[0:6:2] * z): .4f}")
print(f"Shifts in y-Axis (mm): {result[1]: .4f} \t {result[3]: .4f} \t {result[5]: .4f} \t {sum(result[1:6:2]): .4f} \t {sum(result[1:6:2] * z): .4f}")
print(f"Lagrange Multipliers:  {result[6]: .3f}, {result[7]: .3f}, {result[8]: .3f}, {result[9]: .3f}")

_result = copy(result)
_result[0:6:2] += x_shift_global + x_shear
_result[1:6:2] += y_shift_global + y_shear
print("Test for Closure, apply the implicit global shifts and shears, which the alignment centers around 0:")
print(f"                       Layer 0 \t Layer 1 \t Layer 2")
print(f"Shifts in x-Axis (mm): {_result[0]: .4f} \t {_result[2]: .4f} \t {_result[4]: .4f}")
print(f"Shifts in y-Axis (mm): {_result[1]: .4f} \t {_result[3]: .4f} \t {_result[5]: .4f}")

mod_sp_aligned = mod_sp.copy()
mod_sp_aligned["x0"] -= result[0]
mod_sp_aligned["x1"] -= result[2]
mod_sp_aligned["x2"] -= result[4]
mod_sp_aligned["y0"] -= result[1]
mod_sp_aligned["y1"] -= result[3]
mod_sp_aligned["y2"] -= result[5]<jupyter_output>Alignment Constants
                       Layer 0 	 Layer 1 	 Layer 2 	 Global 	 Shear 
Shifts in x-Axis (mm): -0.0004 	  0.0009 	 -0.0004 	 -0.0000 	 -0.0000
Shifts in y-Axis (mm): -0.0150 	  0.0300 	 -0.0150 	 -0.0000 	 -0.0000
Lagrange Multipliers:  -0.000,  0.000, -0.000, -0.000
Test for Closure, apply the implicit global shifts and shears, which the alignment centers around 0:
                       Layer 0 	 Layer 1 	 Layer 2
Shifts in x-Axis (mm):  0.0996 	  0.1009 	  0.0996
Shifts in y-Axis (mm): -0.0600 	  0.0300 	  0.0300<jupyter_text>Use this for only rotation<jupyter_code>z = np.array([mod_sp.z0.unique()[0], mod_sp.z1.unique()[0], mod_sp.z2.unique()[0]])

print("Alignment Constants")
print(f"                        Layer 0 \t Layer 1 \t Layer 2 \t Global \t Shear/Curl ")
print(f"Rotation z-Axis (mrad): {result[0]*1000: .4f} \t {result[1]*1000: .4f} \t {result[2]*1000: .4f} \t {sum(result[0:3]*1000): .4f} \t {sum(result[0:3]*1000 * z): .4f}")
print(f"Lagrange Multipliers: {result[0]: .3f}, {result[1]: .3f}")


mod_sp_aligned = mod_sp.copy()
# Apply rotation alignment proper
# mod_sp_aligned["x0"] = mod_sp_aligned[["x0", "y0"]].apply(lambda x: np.cos(-result[2]) * x[0] - np.sin(-result[2]) * x[1], axis=1)
# mod_sp_aligned["y0"] = mod_sp_aligned[["x0", "y0"]].apply(lambda x: np.sin(-result[2]) * x[0] + np.cos(-result[2]) * x[1], axis=1)
# mod_sp_aligned["x1"] = mod_sp_aligned[["x1", "y1"]].apply(lambda x: np.cos(-result[5]) * x[0] - np.sin(-result[5]) * x[1], axis=1)
# mod_sp_aligned["y1"] = mod_sp_aligned[["x1", "y1"]].apply(lambda x: np.sin(-result[5]) * x[0] + np.cos(-result[5]) * x[1], axis=1)
# mod_sp_aligned["x2"] = mod_sp_aligned[["x2", "y2"]].apply(lambda x: np.cos(-result[8]) * x[0] - np.sin(-result[8]) * x[1], axis=1)
# mod_sp_aligned["y2"] = mod_sp_aligned[["x2", "y2"]].apply(lambda x: np.sin(-result[8]) * x[0] + np.cos(-result[8]) * x[1], axis=1)

# Apply rotation alignment linearized
mod_sp_aligned["x0"] = mod_sp_aligned[["x0", "y0"]].apply(lambda x: x[0] - (-result[0]) * x[1], axis=1)
mod_sp_aligned["y0"] = mod_sp_aligned[["x0", "y0"]].apply(lambda x: -result[0] * x[0] + x[1], axis=1)
mod_sp_aligned["x1"] = mod_sp_aligned[["x1", "y1"]].apply(lambda x: x[0] - (-result[1]) * x[1], axis=1)
mod_sp_aligned["y1"] = mod_sp_aligned[["x1", "y1"]].apply(lambda x: -result[1] * x[0] + x[1], axis=1)
mod_sp_aligned["x2"] = mod_sp_aligned[["x2", "y2"]].apply(lambda x: x[0] - (-result[2]) * x[1], axis=1)
mod_sp_aligned["y2"] = mod_sp_aligned[["x2", "y2"]].apply(lambda x: -result[2] * x[0] + x[1], axis=1)<jupyter_output><empty_output><jupyter_text>Use this for xy + rotation<jupyter_code>z = np.array([mod_sp.z0.unique()[0], mod_sp.z1.unique()[0], mod_sp.z2.unique()[0]])

print("Alignment Constants")
print(f"                        Layer 0 \t Layer 1 \t Layer 2 \t Global \t Shear/Curl ")
print(f"Shifts in x-Axis (mm) : {result[0]: .4f} \t {result[3]: .4f} \t {result[6]: .4f} \t {sum(result[0:9:3]): .4f} \t {sum(result[0:9:3] * z): .4f}")
print(f"Shifts in y-Axis (mm) : {result[1]: .4f} \t {result[4]: .4f} \t {result[7]: .4f} \t {sum(result[1:9:3]): .4f} \t {sum(result[1:9:3] * z): .4f}")
print(f"Rotation z-Axis (mrad): {result[2]*1000: .4f} \t {result[5]*1000: .4f} \t {result[8]*1000: .4f} \t {sum(result[2:9:3]*1000): .4f} \t {sum(result[2:9:3]*1000 * z): .4f}")
print(f"Lagrange Multipliers: {result[9]: .3f}, {result[10]: .3f}, {result[11]: .3f}, {result[12]: .3f},")# {result[13]: .3f}, {result[14]: .3f}")

_result = copy(result)
_result[0:N_GLOBAL:3] += x_shift_global + x_shear
_result[1:N_GLOBAL:3] += y_shift_global + y_shear
_result[2:N_GLOBAL:3] += gamma_global
print()
print("Test for Closure, apply the implicit global shifts and shears, which the alignment centers around 0:")
print(f"                        Layer 0 \t Layer 1 \t Layer 2")
print(f"Shifts in x-Axis (mm) : {_result[0]: .4f} \t {_result[3]: .4f} \t {_result[6]: .4f}")
print(f"Shifts in y-Axis (mm) : {_result[1]: .4f} \t {_result[4]: .4f} \t {_result[7]: .4f}")
print(f"Rotation z-Axis (mrad): {_result[2]*1000: .4f} \t {_result[5]*1000: .4f} \t {_result[8]*1000: .4f}")

mod_sp_aligned = mod_sp.copy()
mod_sp_aligned["x0"] -= result[0]
mod_sp_aligned["x1"] -= result[3]
mod_sp_aligned["x2"] -= result[6]
mod_sp_aligned["y0"] -= result[1]
mod_sp_aligned["y1"] -= result[4]
mod_sp_aligned["y2"] -= result[7]

# Apply rotation alignment proper
# mod_sp_aligned["x0"] = mod_sp_aligned[["x0", "y0"]].apply(lambda x: np.cos(-result[2]) * x[0] - np.sin(-result[2]) * x[1], axis=1)
# mod_sp_aligned["y0"] = mod_sp_aligned[["x0", "y0"]].apply(lambda x: np.sin(-result[2]) * x[0] + np.cos(-result[2]) * x[1], axis=1)
# mod_sp_aligned["x1"] = mod_sp_aligned[["x1", "y1"]].apply(lambda x: np.cos(-result[5]) * x[0] - np.sin(-result[5]) * x[1], axis=1)
# mod_sp_aligned["y1"] = mod_sp_aligned[["x1", "y1"]].apply(lambda x: np.sin(-result[5]) * x[0] + np.cos(-result[5]) * x[1], axis=1)
# mod_sp_aligned["x2"] = mod_sp_aligned[["x2", "y2"]].apply(lambda x: np.cos(-result[8]) * x[0] - np.sin(-result[8]) * x[1], axis=1)
# mod_sp_aligned["y2"] = mod_sp_aligned[["x2", "y2"]].apply(lambda x: np.sin(-result[8]) * x[0] + np.cos(-result[8]) * x[1], axis=1)

# Apply rotation alignment linearized
mod_sp_aligned["x0"] = mod_sp_aligned[["x0", "y0"]].apply(lambda x: x[0] - (-result[2]) * x[1], axis=1)
mod_sp_aligned["y0"] = mod_sp_aligned[["x0", "y0"]].apply(lambda x: -result[2] * x[0] + x[1], axis=1)
mod_sp_aligned["x1"] = mod_sp_aligned[["x1", "y1"]].apply(lambda x: x[0] - (-result[5]) * x[1], axis=1)
mod_sp_aligned["y1"] = mod_sp_aligned[["x1", "y1"]].apply(lambda x: -result[5] * x[0] + x[1], axis=1)
mod_sp_aligned["x2"] = mod_sp_aligned[["x2", "y2"]].apply(lambda x: x[0] - (-result[8]) * x[1], axis=1)
mod_sp_aligned["y2"] = mod_sp_aligned[["x2", "y2"]].apply(lambda x: -result[8] * x[0] + x[1], axis=1)<jupyter_output><empty_output><jupyter_text>Post-Alignment<jupyter_code>%%time
mod_sp_aligned = fit_and_calculate_residuals(mod_sp_aligned)
plot_metrics(mod_sp_aligned)
if mc_truth is not None:
    plot_mc_residuals(mod_sp_aligned)
plt.scatter(mod_sp_orig.x1, mod_sp_orig.y1, marker='.', label="Original")
plt.scatter(mod_sp_aligned.x1, mod_sp_aligned.y1, marker='.', label="Aligned") 
plt.xlim(-150,150)
plt.ylim(-150,150)
plt.legend()
plt.show()
plt.close()

plt.scatter(mod_sp_orig.x1, mod_sp_orig.y1, marker='.', label="Original")
plt.scatter(mod_sp.x1, mod_sp.y1, marker='.', label="Mis-Aligned") 
plt.xlim(-150,150)
plt.ylim(-150,150)
plt.legend()
plt.show()
plt.close()

plt.scatter(mod_sp.x1, mod_sp.y1, marker='.', label="Mis-Aligned")
plt.scatter(mod_sp_aligned.x1, mod_sp_aligned.y1, marker='.', label="Aligned") 
plt.xlim(-150,150)
plt.ylim(-150,150)
plt.legend()
plt.show()
plt.close()
assert False, "Anything below is not really doing anything"<jupyter_output><empty_output><jupyter_text>IterationsDoes not do anything really.<jupyter_code>for i in range(1000):
    result = solve_global_parameters(mod_sp_aligned, n_layers=3, verbose=False)
    
    z = np.array([mod_sp.z0.unique()[0], mod_sp.z1.unique()[0], mod_sp.z2.unique()[0]])

    print("Alignment Constants")
    print(f"                        Layer 0 \t Layer 1 \t Layer 2 \t Global \t Shear/Curl ")
    print(f"Shifts in x-Axis (mm) : {result[0]: .4f} \t {result[3]: .4f} \t {result[6]: .4f} \t {sum(result[0:9:3]): .4f} \t {sum(result[0:9:3] * z): .4f}")
    print(f"Shifts in y-Axis (mm) : {result[1]: .4f} \t {result[4]: .4f} \t {result[7]: .4f} \t {sum(result[1:9:3]): .4f} \t {sum(result[1:9:3] * z): .4f}")
    print(f"Rotation z-Axis (mrad): {result[2]*1000: .4f} \t {result[5]*1000: .4f} \t {result[8]*1000: .4f} \t {sum(result[2:9:3]*1000): .4f} \t {sum(result[2:9:3]*1000 * z): .4f}")
    print(f"Lagrange Multipliers: {result[9]: .3f}, {result[10]: .3f}, {result[11]: .3f}, {result[12]: .3f}, {result[13]: .3f}, {result[14]: .3f}")

    _result = copy(result)
    _result[0:N_GLOBAL:3] += x_shift_global + x_shear
    _result[1:N_GLOBAL:3] += y_shift_global + y_shear
    _result[2:N_GLOBAL:3] += gamma_global
    print()
    print("Test for Closure, apply the implicit global shifts and shears, which the alignment centers around 0:")
    print(f"                        Layer 0 \t Layer 1 \t Layer 2")
    print(f"Shifts in x-Axis (mm) : {_result[0]: .4f} \t {_result[3]: .4f} \t {_result[6]: .4f}")
    print(f"Shifts in y-Axis (mm) : {_result[1]: .4f} \t {_result[4]: .4f} \t {_result[7]: .4f}")
    print(f"Rotation z-Axis (mrad): {_result[2]*1000: .4f} \t {_result[5]*1000: .4f} \t {_result[8]*1000: .4f}")
    
    mod_sp_aligned["x0"] -= result[0]
    mod_sp_aligned["x1"] -= result[3]
    mod_sp_aligned["x2"] -= result[6]
    mod_sp_aligned["y0"] -= result[1]
    mod_sp_aligned["y1"] -= result[4]
    mod_sp_aligned["y2"] -= result[7]

    mod_sp_aligned["x0"] = mod_sp_aligned[["x0", "y0"]].apply(lambda x: np.cos(-result[2]) * x[0] - np.sin(-result[2]) * x[1], axis=1)
    mod_sp_aligned["y0"] = mod_sp_aligned[["x0", "y0"]].apply(lambda x: np.sin(-result[2]) * x[0] + np.cos(-result[2]) * x[1], axis=1)
    mod_sp_aligned["x1"] = mod_sp_aligned[["x1", "y1"]].apply(lambda x: np.cos(-result[5]) * x[0] - np.sin(-result[5]) * x[1], axis=1)
    mod_sp_aligned["y1"] = mod_sp_aligned[["x1", "y1"]].apply(lambda x: np.sin(-result[5]) * x[0] + np.cos(-result[5]) * x[1], axis=1)
    mod_sp_aligned["x2"] = mod_sp_aligned[["x2", "y2"]].apply(lambda x: np.cos(-result[8]) * x[0] - np.sin(-result[8]) * x[1], axis=1)
    mod_sp_aligned["y2"] = mod_sp_aligned[["x2", "y2"]].apply(lambda x: np.sin(-result[8]) * x[0] + np.cos(-result[8]) * x[1], axis=1)

    mod_sp_aligned["r_x0"] = mod_sp_aligned["x0"] / sigma_x
    mod_sp_aligned["r_x1"] = mod_sp_aligned["x1"] / sigma_x
    mod_sp_aligned["r_x2"] = mod_sp_aligned["x2"] / sigma_x
    mod_sp_aligned["r_y0"] = mod_sp_aligned["y0"] / sigma_y
    mod_sp_aligned["r_y1"] = mod_sp_aligned["y1"] / sigma_y
    mod_sp_aligned["r_y2"] = mod_sp_aligned["y2"] / sigma_y
    
    

plt.scatter(mod_sp_orig.x1, mod_sp_orig.y1, marker='.', label="Original")
plt.scatter(mod_sp_aligned.x1, mod_sp_aligned.y1, marker='.', label="Aligned") 
plt.xlim(-150,150)
plt.ylim(-150,150)
plt.legend()
plt.show()
plt.close()
%%time
plot_metrics(mod_sp_aligned)
if mc_truth is not None:
    plot_mc_residuals(mod_sp_aligned)<jupyter_output><empty_output><jupyter_text>Rotations ...Just testing some things. We define 3 points along a straight line<jupyter_code>p0 = np.array((1, 0, 1))
p1 = np.array((1, 0, 2))
p2 = np.array((1, 0, 3))
def format_coordinates(p0, p1, p2):
    x = np.array([p0[0], p1[0], p2[0]])
    y = np.array([p0[1], p1[1], p2[1]])
    z = np.array([p0[2], p1[2], p2[2]])
    return x, y, z
def dr_dalpha(z):
    drdalpha_x = np.array([1, z, 0, 0])
    drdalpha_y = np.array([0, 0, 1, z])
    return np.array([
        drdalpha_x,
        drdalpha_y,
    ])

def Gamma(z):
    return sum(dr_dalpha(_z).transpose() @ dr_dalpha(_z) for _z in z)

def beta(x, y, z):
    return sum(dr_dalpha(_z).transpose() @ np.array([_x, _y]) for _x, _y, _z in zip(x, y, z))<jupyter_output><empty_output><jupyter_text>Nominal (No Rotation)Fit the straight line with 2 layers, check if the third point is still on the straight line.<jupyter_code>x, y, z = format_coordinates(p0, p1, p2)
x, y, z
track_parameters = np.linalg.inv(Gamma(z[:2])) @ beta(x[:2], y[:2], z[:2])
space_point_prediction(z[2], *track_parameters).transpose(), p2[:2]<jupyter_output><empty_output><jupyter_text>Real RotationApply real rotations, relative to the distance in z.<jupyter_code>U = lambda x: np.array([
    [np.cos(x), -np.sin(x), 0],
    [np.sin(x), np.cos(x), 0],
    [0, 0, 1],
])
rotation = 1e-3
x, y, z = format_coordinates(U(0*rotation) @ p0, U(1*rotation) @ p1, U(2*rotation) @ p2)
x, y, z
track_parameters = np.linalg.inv(Gamma(z[:2])) @ beta(x[:2], y[:2], z[:2])
space_point_prediction(z[2], *track_parameters).transpose(), x[2], y[2]<jupyter_output><empty_output><jupyter_text>Linearized RotationApply linearized rotations, relative to the distance in z.<jupyter_code>V = lambda x: np.array([
    [1, -x, 0],
    [x, 1, 0],
    [0, 0, 1]
])
rotation = 1e-3
x, y, z = format_coordinates(V(0*rotation) @ p0, V(1*rotation) @ p1, V(2*rotation) @ p2)
x, y, z
track_parameters = np.linalg.inv(Gamma(z[:2])) @ beta(x[:2], y[:2], z[:2])
space_point_prediction(z[2], *track_parameters).transpose(), x[2], y[2]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20751950568152067, 'fraction_numerical': 0.050007854636853956, 'mean_word_length': 4.06363515842503, 'pattern_counts': {'":': 20, '<': 91, '<?xml version=': 0, '>': 92, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '23', 'dir': 'jupyter-structured-clean-dedup', 'id': '553407', 'license': 'MIT', 'n_tokens_mistral': 17404, 'n_tokens_neox': 16478, 'n_words': 3346, 'path': 'Full-Alignment.ipynb', 'repo_name': 'MarkusPrim/faser-alignment', 'size': '38194'}