Adding BPMs errors in measuring the optics

[elafmusa]

I want to include the effect of BPMs errors in the measured optics, for that i defined a parameter Systematic_error to the BPM definition that I can access and change later :

BPM = elements.Monitor('BPM', K = 0.0, Systematic_error = [0.01, 0.01])

I measured the beam optics using the following, my question is how to add the BPM errors effect in the case of dispersion, and beta, as i did in closed orbit by adding the errors on top to the measured orbit:

bpm_indexes = get_refpts(ring, elements.Monitor)
lindata0, tune, chrom, lindata = ring.linopt(get_chrom=True, refpts=bpm_indexes)

Eta_xx1 = lindata['dispersion'][:, 0]
Eta_yy1 = lindata['dispersion'][:, 2]
closed_orbitx1 = lindata['closed_orbit'][:, 0]
closed_orbity1 = lindata['closed_orbit'][:, 2]
for i in range(len(bpm_indexes)):
Sys_err = (ring[bpm_indexes[i]].Systematic_error)

    closed_orbitx2 = closed_orbitx1[i] +  Sys_err[0]*(1+ rand * np.random.randn())
    closed_orbity2 = closed_orbity1[i] + Sys_err[1]*(1+ rand * np.random.randn())

[simoneliuzzo]

Dear @elafmusa,

in matlab I used to have a wrapper function for findorbit* called findorbit*Err here how it looked like similar to yours.
I would always use findorbit*Err.m in the code, rather than findorbit*.m

The LOCO-like corrections are all based on orbit reading (ORM, dispersion). Thus modifying this single function had the desired result.

Concerning a modification of the beta and dispersino output of at.linopt (or matlab atlinopt) I think it is not needed. The beta functions and dispersion are not affected by measurement errors. What is affected is only the measurement itself, so orbit vs RF shift and orbit vs steerers kicks.

Also in a real machine, setting random or systematic BPM errors will have no impact on the optics. The impact will be visible only when the BPM measurement are used for some correction.

What you will observe is that in presence of BPM errors the performance of the optics correction will be limited by the BPM errors, due to the worst measurement. When looking at the residual beta beating and dispersion deviation in your model after correction, in presence of BPM errors, you will likely converge to a larger beta-beating and dispersion residual.

hope it is of some help

best regards
Simone

[simoneliuzzo]

I would also like to mention that we may look for a solution to include BPM errors in pyAT in a more solid way.

I think @swhite2401 and @lfarv may have already thought of this and may be have good ideas on how to do this.

It would be a great addition to AT.

I remember discussions concerning a BPM passmethod, but it is not as simple as it seems, since all coordinate transformations are then passed to the next element.

We may think for example to update findorbit.m and at.find_orbit to always consider BPM errors if they are available

A fully loaded BPM would then look like:

         FamName: 'bpm'
      PassMethod: 'IdentityPass'
          Length: 0
           Class: 'Monitor'
     OffsetError: [0 0]          # offset the closed orbit
    ReadingError: [0 0]        # random value added at each `findorbit* `
      ScaleError: [0 0]          # scale factor to bpm reading
        Rotation: [2×2 double]    # BPM rotation
              T1: [6×1 double]        # translation of the whole BPM, the golden orbit at this BPM. 
              T2: [6×1 double]
              R1: [6×6 double]       # rotation of the whole BPM, also to be removed from the real bpm reading for optimal corrections
              R2: [6×6 double]

[lnadolski]

hi @simoneliuzzo

Just a few questions, what is the
1/ Difference between Rotation and R1(R2)?
2/ is Offfset for BBA for instance and T1/T2 for alignment errors?
3/ BPM errors (Scale, Offset) depends of the mode (slow. acquisition like orbit) et Turn by Turn data for instance

[lfarv]

@lnadolski:
For you point 1, you are right, it' similar. But if you apply the translation on the BPM, you'll never see the result ! To get it, you should set a translation at the exit of the previous element, and the opposite one at the entrance of the next element. Not really convenient…

So for clarity, I am in favour of @simoneliuzzo's solution: postprocessing the results of find_orbit in an enclosing function. Whether a general method can cover all error cases is still to be seen.

[simoneliuzzo]

Dear @lnadolski,

I have been running simulations with rather large machine misalignment following a smooth curve.

An example in the figure below:

The blue crosses represent the set misalignment
The red dots are the closed orbit
The black circles (close to zero) the orbit that we really want to correct.

In this case the BPMs are also misaligned, and their misalignment is naturally defined in AT by T1, T2, R1 and R2.
T1, T2, R1 and R2 have no impact on the BPM (IdentityPass), nor on the orbit provided by AT (the red dots in the figure above).

The target orbit for a machine with such large errors is not zero, but the orbit that follows in the best way such large misalignments.
So, if correctly set, T1,T2, R1 and R2 may be used to represent for the BPMs the best knowledge in a machine with errors of the wished final position of the beam. While additional fields such as Offset, scale, Rotation, Reading represent instead the errors of the BPM itself.

Once the fields Offset/Rotation/Scale/Reading are set, they may take values corresponding to turn by turn or to slow acquisition.
I have no objection to consider duplicate all the additional fields to consider these two cases independently.

[simoneliuzzo]

Dear @elafmusa @lfarv and @lnadolski ,

I have these functions that I use in the pyAT ring simulator Tango Device server. They are copies of the matlab equivalents written originally by @lfarv. Rather low quality code, sorry.

https://gitlab.esrf.fr/accelerators/Simulators/EbsSimulator/optics_for_pyringsimulator/-/blob/master/pyringsimulator/bpm_modifications.py

USAGE:

 rel, tel, trand = bpm_modifications.bpm_matrices( ring, ind_bpms)

 h_err, v_err = bpm_modifications.bpm_process(
                    lb['closed_orbit'][:, 0],
                    lb['closed_orbit'][:, 2],
                    rel, tel, trand)

Fucntion above in case the link is not accessible

from builtins import isinstance, hasattr

import at
import numpy as np


def bpm_matrices(ring,ind_bpms):
    # read BPM matrices from lattice
    rel = []
    tel = []
    trand = []
    t1 = [0.0, 0.0]
    r1 = [[0.0, 0.0], [0.0, 0.0]]
    tb = [0.0, 0.0]
    rb = [[0.0, 0.0], [0.0, 0.0]]
    t_rand = [0.0, 0.0]

    for ib in ind_bpms:

        if hasattr(ring[ib], 'T1'):
            ii = [0, 2]
            t1 = [ring[ib].T1[iii] for iii in ii]
        else:
            t1 = [0.0, 0.0]

        if hasattr(ring[ib], 'R1'):
            r1[0][0] = ring[ib].R1[0][0]
            r1[0][1] = ring[ib].R1[0][2]
            r1[1][0] = ring[ib].R1[2][0]
            r1[1][1] = ring[ib].R1[2][2]

        else:
            r1 = [[1.0, 0.0], [0.0, 1.0]]

        if hasattr(ring[ib], 'Rotation'):
            if isinstance(ring[ib].Rotation, np.ndarray) or \
                isinstance(ring[ib].Rotation, np.float64 or \
                isinstance(ring[ib].Rotation, np.uint8)):

                dimr = ring[ib].Rotation.size
            else:
                dimr = len(ring[ib].Rotation)

            if dimr > 1:
                rb[0][0] = np.cos(ring[ib].Rotation[0])
                rb[0][1] = np.sin(ring[ib].Rotation[0])
                rb[1][0] = -np.sin(ring[ib].Rotation[0])
                rb[1][1] = np.cos(ring[ib].Rotation[0])
            else:
                rb[0][0] = np.cos(ring[ib].Rotation)
                rb[0][1] = np.sin(ring[ib].Rotation)
                rb[1][0] = -np.sin(ring[ib].Rotation)
                rb[1][1] = np.cos(ring[ib].Rotation)
        else:
            rb = [[1.0, 0.0], [0.0, 1.0]]

        if hasattr(ring[ib], 'Offset'):
            tb = ring[ib].Offset
        else:
            tb = [0.0, 0.0]

        if hasattr(ring[ib], 'Scale'):
            scale = ring[ib].Scale
        else:
            scale = [1.0, 1.0]

        if hasattr(ring[ib], 'Reading'):
            t_rand = ring[ib].Reading
        else:
            t_rand = [0.0, 0.0]

        t_rand = np.array(t_rand)
        t1 = np.array(t1)
        tb = np.array(tb)
        r1 = np.array(r1)
        rb = np.array(rb)
        scale = np.array(scale)

        trand.append(t_rand)

        r_mat_prod = np.matmul(r1, rb)

        """
        print('r_mat_prod')
        print(r_mat_prod)
        print('r_mat * t1')
        print(t1)
        print(np.matmul(r_mat_prod, t1))
        print('tb')
        print(tb)
        print('r_mat @ t1 + tb')
        print(((r_mat_prod @ t1) + tb))
        print('scale')
        print(scale)
        print('scale * ((r_mat_prod @ t1) + tb) ')
        print(scale * ((r_mat_prod @ t1) + tb) )
        """

        tt = (scale * ( (r_mat_prod @ t1) + tb)).transpose()
        # print('tt')
        # print(tt)

        tel.append(tt)

        # print('(np.tile(scale, 2, 1))')
        # print((np.tile(scale, (2,1))))
        # print('np.tile(scale, 2, 1) * r_mat_prod)')
        # print((np.tile(scale, (2,1)) * r_mat_prod))
        rel.append((np.tile(scale, (2,1)) * r_mat_prod))

    return rel, tel, trand


def bpm_process(hor_pos, ver_pos, rel, tel, trand):
    # modify bpm readings according to offset, rotation and scale factors in ring
    new_hor_pos = []
    new_ver_pos = []

    pos = np.column_stack((hor_pos, ver_pos))
    # print(len(pos))

    for p, trasl, scal, ran in zip(pos, tel, rel, trand):
        # print('pos: {} , scale: {} , trasl {}, random {}'.format(p, scal, trasl, ran))
        # print(scal @ p)
        # print(scal @ p + trasl)
        new_pos =  scal @ p + trasl + ran*np.random.normal((2,1))
        # print('new pos: {} '.format(new_pos))
        new_hor_pos.append(new_pos[0])
        new_ver_pos.append(new_pos[1])

    return new_hor_pos, new_ver_pos


def bpm_matrices_direct(ring, ind_bpms, hor, ver):
    # read BPM matrices from lattice
    rel = []
    tel = []
    trand = []
    t1 = []
    r1 = []
    tb = []
    rb = []
    t_rand = []
    new_h = []
    new_v = []

    for ib, h, v in zip(ind_bpms, hor, ver):

        if hasattr(ring[ib], 'T1'):
            t1 = ring[ib].T1
        else:
            t1 = [0.0, 0.0]

        if hasattr(ring[ib], 'R1'):
            r1[0][0] = ring[ib].R1[0][0]
            r1[0][1] = ring[ib].R1[0][2]
            r1[1][0] = ring[ib].R1[2][0]
            r1[1][1] = ring[ib].R1[2][2]

        else:
            r1 = [[1.0, 0.0], [0.0, 1.0]]

        if hasattr(ring[ib], 'Rotation'):
            rb[0][0] = np.cos(ring[ib].Rotation[0])
            rb[0][1] = np.sin(ring[ib].Rotation[0])
            rb[1][0] = -np.sin(ring[ib].Rotation[0])
            rb[1][1] = np.cos(ring[ib].Rotation[0])

        else:
            rb = [[1.0, 0.0], [0.0, 1.0]]

        if hasattr(ring[ib], 'Offset'):
            tb = ring[ib].Offset
        else:
            tb = [0.0, 0.0]

        if hasattr(ring[ib], 'Scale'):
            scale = ring[ib].Scale
        else:
            scale = [1.0, 1.0]

        if hasattr(ring[ib], 'Reading'):
            t_rand = ring[ib].Reading
        else:
            t_rand = [0.0, 0.0]

        t_rand = np.array(t_rand)
        t1 = np.array(t1)
        tb = np.array(tb)
        r1 = np.array(r1)
        rb = np.array(rb)
        scale = np.array(scale)

        trand.append(t_rand)

        r_mat_prod = np.matmul(r1, rb)

        tel.append(np.dot(scale, np.matmul(r_mat_prod, t1) + tb))

        rel.append(np.dot(np.array([scale, scale]), r_mat_prod))

        new_h.append(h)
        new_v.append(v)

    return new_h, new_v


def example_bpm_err():

    import at
    import time

    # path_to_lattice_model = "/machfs/liuzzo/EBS/beamdyn/matlab/optics/sr/theory/betamodel.mat"
    # ring = at.load_lattice(path_to_lattice_model, key='betamodel')

    path_to_lattice_model = "/machfs/liuzzo/EBS/minibeta_ID31_x7m_y1m/ANALYSIS/DA_TLT_IE/" \
                            "minibeta_ID31_x7m_y1m_53_30_betamodel/" \
                            "RefLatModelminibeta_ID31_x7m_y1m_53_30_betamodel_ns10/" \
                            "Errors_0003.mat"

    # path_to_lattice_model = "/machfs/liuzzo/EBS/S28F_all_BM_22Aug2021/ANALYSIS/" \
    #                        "DA_TLT_IE_ErrorsCorrection/S28F_all_BM_22Aug2021_53_30/" \
    #                        "RefLatModelS28F_all_BM_22Aug2021_53_30_ns10/" \
    #                        "Errors_0002.mat"

    ring = at.load_lattice(path_to_lattice_model, key='rcor')

    ring = at.Lattice(ring)

    # add some orbit
    ring[15].PolynomB[0] = 20e-6
    ring[15].PolynomA[0] = 20e-6
    print(ring[15])
    time.sleep(1)

    elems = ring.get_elements('BPM*')
    ind_bpms = [ring.index(elem) for elem in elems]

    # introduce known errors at given BPMs
    print(ind_bpms[45])
    print(ring[ind_bpms[45]])

    ring[ind_bpms[45]].Offset = [-234e-6, 334e-6]

    ring[ind_bpms[318]].Offset = [-24e-6, -334e-6]

    ring[ind_bpms[300]].Rotation = [3.14]

    ring[ind_bpms[200]].Scale = [1.1, 0.6]

    ring[ind_bpms[100]].Reading = [100e-6, 60e-6]

    # ring.radiation_off()
    # _orb_, orb = at.find_sync_orbit(ring, refpts=ind_bpms)

    ring.radiation_on()
    _orb_, orb = at.find_orbit6(ring, refpts=ind_bpms)

    rel, tel, trand = bpm_matrices(ring, ind_bpms)

    h = orb[:, 0]
    v = orb[:, 2]

    print(len(h))
    print(len(v))
    print(len(tel))
    print(len(rel))
    print(len(trand))

    he, ve = bpm_process(h, v, rel, tel, trand)

    # read again to see effect of Reading error

    # _orb_, orb = at.find_sync_orbit(ring, refpts=ind_bpms)

    ring.radiation_on()
    _orb_, orb = at.find_orbit6(ring, refpts=ind_bpms)

    h2 = orb[:, 0]
    v2 = orb[:, 2]
    he2, ve2 = bpm_process(h2, v2, rel, tel, trand)

    # plot effect of errors
    import matplotlib.pyplot as plt

    fig, (ax1, ax2) = plt.subplots(2, 1)
    ax1.plot(h, label='Horizontal initial')
    ax1.plot(he, label='Error')
    ax1.plot(he2, label='Error2')
    ax1.legend()

    ax2.plot(v, label='vertical initial')
    ax2.plot(ve, label='Error')
    ax2.plot(ve2, label='Error2')
    ax2.legend()

    plt.show()


if __name__ == '__main__':

    # to test:
    # >>> python
    # >>> from pyringsimulator import bpm_modifications as bpm_mod
    # >>> bpm_mod.example_bpm_err()

    # example_bpm_err()
    pass