Pythonic approach simple
Introduction
In this notebook, we are aiming to use python as scripting language for the MAD-X mask. The goal is to bring together the scripting flexibility of python and the optics capability from MAD-X.
The rationale is to use cpymad to interface python to MAD-X.
Beyond the standard python setup, to run this notebook you need the following packages
In addition you need AFS mounted and access to
/afs/cern.ch/eng/lhc/optics/
/afs/cern.ch/eng/tracking_tools
/afs/cern.ch/user/s/sterbini/public/tracking_tools # to be replaced once consolidated
Warning
We assume that the user knows MAD-X and python.
Danger
Values chosen are arbitrary!
Importing the packages
Here you are the "special" packages
from madxp import luminosity as lumi
from cpymad.madx import Madx
from madxp import cpymadTool as mt
import fillingpatterns as fp
This are standard packages.
from collections import OrderedDict
import numpy as np
import pandas as pd
from matplotlib import pylab as plt
from IPython.display import display
import time
import os
import warnings
import shutil
import urllib.request, json
warnings.filterwarnings('always')
To import the filling pattern you can use the approach proposed here.
Warning
Please note that the chose pattern is arbitratry.
with urllib.request.urlopen('https://raw.githubusercontent.com/PyCOMPLETE/FillingPatterns/master/examples/25ns_2744b_2736_2246_2370_240bpi_13inj_800ns_bs200ns_BCMS_5x48b.json') as url:
data = json.loads(url.read().decode())
bb_pattern = fp.FillingPattern(data['beam1'], data['beam2'])
Mask parameters
In the parameter_dict we will store the parameters needed for simulation.
Warning
Remember that MAD-X does not accept "strings" as variables.
# %% Definition of the parameters that are not knobs of the beam sequence (no strings please!)
parameter_dict={
# =============================================================================
# Beam parameters
# =============================================================================
## LHC beam 1 (clockwise), LHC beam 2 (clockwise), LHC beam 2 (counterclockwise)
'par_mylhcbeam': 1,
## beam normalized emittance [m rad]
'par_beam_norm_emit': 2.5e-6,
## [m]
'par_beam_sigt': 0.075,
## [-]
'par_beam_sige': 1.1e-4,
## [-]
'par_beam_npart': 1.16e11,
## [GeV]
'par_beam_energy_tot': 7000,
## [A]
'par_oct_current': 350,
## [-]
'par_chromaticity': 15,
## [MV]
'par_vrf_total': 16.,
## Tunes with fractional part
'par_qx0': 62.31, 'par_qy0': 60.32,
# =============================================================================
# Beam-Beam configuration
# =============================================================================
## install the BB elements [0,1]
'par_on_bb_switch': 1,
## if 1 lumi leveling in ip8 is applied and q/q' match is done with bb off [0,1]
'par_on_collision': 1,
## bunch separation [ns]
'par_b_t_dist': 25.,
## default value for the number of additionnal parasitic encounters inside D1
'par_n_inside_D1': 5,
## number of slices for head-on in IR1 [between 0 and 201]
'par_nho_IR1': 11, 'par_nho_IR2': 11, 'par_nho_IR5': 11, 'par_nho_IR8': 11,
## flag to install the Crab Cavities [0, 1]
'par_install_crabcavities': 0,
# can be negative positive or zero to switch of spectr in lhcb
'par_lhcb_polarity': 1.,
# =============================================================================
# Leveling in IP8
# =============================================================================
# leveled luminosity in IP8 (considered if par_on_collision=1) [Hz/cm2]
'par_lumi_ip8': 2e33,
# These variables define the number of Head-On collisions in the 4 IPs
'par_nco_IP1': bb_pattern.n_coll_ATLAS,
'par_nco_IP2': bb_pattern.n_coll_ALICE,
'par_nco_IP5': bb_pattern.n_coll_CMS,
'par_nco_IP8': bb_pattern.n_coll_LHCb,
# =============================================================================
# Errors and corrections
# =============================================================================
# Select seed for errors
'par_myseed': 0,
# Set this flag to correct the errors of D2 in the NLC
# (warning: for now only correcting b3 of D2, still in development)
'par_correct_for_D2': 0,
# Set this flag to correct the errors of MCBXF in the NLC
# (warning: this might be less reproducable in reality, use with care)
'par_correct_for_MCBX': 0,
'par_off_all_errors': 0,
'par_on_errors_LHC': 0,
'par_on_errors_MBH': 0,
'par_on_errors_Q5': 0,
'par_on_errors_Q4': 0,
'par_on_errors_D2': 0,
'par_on_errors_D1': 0,
'par_on_errors_IT': 0,
'par_on_errors_MCBRD': 0,
'par_on_errors_MCBXF': 0,
# =============================================================================
# Additional parameters
# =============================================================================
# parameter for having verbose output [0,1]
'par_verbose': 1,
# definition of the slicefactor used in the makethin
'par_slicefactor': 4,
# number of optics to use
'par_optics_number':30,
# Specify machine version
'ver_lhc_run' : 3, 'ver_hllhc_optics' : 0,
}
The make_sequence function
The user has to define a "make_sequence" function. The function will take the MAD-X handle, the beam number (1, 2 or 4) and, in this case, the slice factor (to make the sequence thin).
Info
It is important to know that we clearly define an interface that isolate the function from MAD-X.
Hint
Take the time to consider the different madx.call.
Info
This function refers to links (e.g., macro.madx) that will be defined later.
def make_sequence(mad, mylhcbeam, slicefactor):
'''
User-defined function to make the Run3 optics.
'''
start_time = time.time()
mad.input('option, -echo,warn, -info;')
# optics dependent macros
mad.call('macro.madx')
# optics dependent macros
mad.call('optics_indep_macros.madx')
assert mylhcbeam in [1, 2, 4], "Invalid mylhcbeam (it should be in [1, 2, 4])"
if mylhcbeam in [1, 2]:
mad.call('optics_runII/2018/lhc_as-built.seq')
else:
mad.call('optics_runII/2018/lhcb4_as-built.seq')
# New IR7 MQW layout and cabling
mad.call('optics_runIII/RunIII_dev/IR7-Run3seqedit.madx')
# Makethin part
if slicefactor > 0:
mad.input(f'slicefactor={slicefactor};') # the variable in the macro is slicefactor
mad.call('optics_runII/2018/toolkit/myslice.madx')
mad.beam()
for my_sequence in ['lhcb1','lhcb2']:
if my_sequence in list(mad.sequence):
mad.input(f'use, sequence={my_sequence}; makethin, sequence={my_sequence}, style=teapot, makedipedge=false;')
else:
warnings.warn('The sequences are not thin!')
# Cycling w.r.t. to IP3 (mandatory to find closed orbit in collision in the presence of errors)
for my_sequence in ['lhcb1','lhcb2']:
if my_sequence in list(mad.sequence):
mad.input(f'seqedit, sequence={my_sequence}; flatten; cycle, start=IP3; flatten; endedit;')
my_output_dict = get_status(mad)
elapsed_time = time.time() - start_time
my_output_dict['elapsed_time'] = elapsed_time
return my_output_dict
The load_optics function
This is the (very simple) function to load the the optics file.
Info
This function refers to links (e.g., optics.madx) that will be defined later.
Hint
We prefer to use links more than variables so that, in the folder of the simulation you have (at least some of) the links used from your mask.
def load_optics(mad):
'''
User-defined function load the optics file.
'''
start_time = time.time()
# nothing very special
mad.call(f'optics.madx')
my_output_dict = get_status(mad)
elapsed_time = time.time() - start_time
my_output_dict['elapsed_time'] = elapsed_time
return my_output_dict
Auxiliary functions
To improve the analysis and debugging of the mask is useful to have some auxiliary functions.
Info
We plan to move them to a package.
Warning
Please check the files that will be removed by the clean_folder function.
def get_status(mad):
'''
Return the status of the variables, sequences, beams and tables of a MAD-X object (mad).
'''
start_time = time.time()
variables=mt.variables_dict(mad)
my_output_dict= {'constant_df': variables['constant_df'],
'independent_variable_df': variables['independent_variable_df'],
'dependent_variable_df': variables['dependent_variable_df'],
'sequences_df': mt.sequences_df(mad),
'beams_df': mt.beams_df(mad),
'tables_list': list(mad.table)}
elapsed_time = time.time() - start_time
my_output_dict['elapsed_time'] = elapsed_time
return my_output_dict
def run_module(mad, module_name):
'''
Run the module_name in the MAD-X object (mad).
'''
start_time = time.time()
mad.call(f'modules/{module_name}')
my_output_dict = get_status(mad)
elapsed_time = time.time() - start_time
my_output_dict['elapsed_time']= elapsed_time
return my_output_dict
def read_parameters(mad, parameter_dict):
'''
Assign the parameter_dict to the MAD-X object (mad).
'''
start_time = time.time()
parameter_dict['par_qx00']=int(parameter_dict['par_qx0'])
parameter_dict['par_qy00']=int(parameter_dict['par_qy0'])
parameter_dict['par_tsplit']=parameter_dict['par_qx00']-parameter_dict['par_qy00']
assert parameter_dict['par_nco_IP5']==parameter_dict['par_nco_IP1']
assert parameter_dict['par_qx00']-parameter_dict['par_qy00']==parameter_dict['par_tsplit']
assert 'par_mylhcbeam' in parameter_dict
assert 'par_beam_norm_emit' in parameter_dict
assert 'par_optics_number' in parameter_dict, 'Optics file not defined.'
for i in parameter_dict:
if isinstance(parameter_dict[i], (float,int)):
mad.input(f'{i}={parameter_dict[i]};')
my_output_dict = get_status(mad)
elapsed_time = time.time() - start_time
my_output_dict['elapsed_time']= elapsed_time
return my_output_dict
def clean_folder(file_string='fc.* *parquet twiss* log.madx stdout.madx bb_lenses.dat last_twiss.0.gz temp', rm_links=True):
'''
Remove the folder from the MAD-X output.
'''
if rm_links:
os.system('find -type l -delete')
os.system('rm -rf '+ file_string)
def check_links(my_path='.'):
'''
Checks the validity of the links.
'''
symlinks = [i for i in os.listdir(my_path)]
for i in symlinks:
if os.path.islink(i):
assert os.path.exists(i), f'Link to {i} is broken.'
Make links
Warning
Establish the correct links is very important since the links are directly used in the python functions and in the MAD-X modules.
Warning
The link to the /afs/cern.ch/user/s/sterbini/public/tracking_tools/modules will be removed once the folder will be merged with te official repository /afs/cern.ch/eng/tracking-tools/modules.
#%% Make links for setting the enviroments
clean_folder()
# Main path
os.symlink('/afs/cern.ch/eng/tracking-tools', 'tracking_tools')
# Mask code folder
os.symlink('/afs/cern.ch/user/s/sterbini/public/tracking_tools/modules', 'modules')
# Machine folder
os.symlink('tracking_tools/machines', 'machines')
# Toolkit folder
os.symlink('tracking_tools/tools', 'tools')
# Beam-beam macros folder
os.symlink('tracking_tools/beambeam_macros', 'beambeam_macros')
# Errors folder
os.symlink('tracking_tools/errors', 'errors')
# RunII optics
os.symlink('/afs/cern.ch/eng/lhc/optics/runII', 'optics_runII')
# RunIII optics
os.symlink('/afs/cern.ch/eng/lhc/optics/runIII', 'optics_runIII')
# Load optics (magnet strengths)
op = int(parameter_dict['par_optics_number'])
os.symlink(f'optics_runIII/RunIII_dev/2022_V1/PROTON/opticsfile.{op}', 'optics.madx')
# General macros
# optics dependent macros
#os.symlink('optics_runII/2018/toolkit/macro.madx', 'macro.madx')
os.symlink('/afs/cern.ch/work/s/sterbini/tracking_tools/tools/macro.madx', 'macro.madx')
# optics independent macros
os.symlink('tools/optics_indep_macros.madx', 'optics_indep_macros.madx')
check_links()
Launch MAD-X from python
command_log_file='log.madx'
stdout_file='stdout.madx'
with open(stdout_file, 'w') as myFile:
mad = Madx(stdout=myFile,command_log=command_log_file)
Preliminary optics checks
In the following cell we execute the functions
- read_parameters
- make_sequence
- load_optics
then we run the module module_00_check_optics.madx.
In a similar way we could ran all the remaing module of the mask. Before executing them we will analyze and do sanity checks on the optics.
Info
While executing the code we build up a work-flow dictionary (my_workflow_dict). Indeed, all the fuction return a dictionary with metadata used for debugging or automatically asserting the work-flow. We will do plenty of example in the following.
# Start making MAD-X operation
my_workflow_dict = OrderedDict()
my_workflow_dict['read_parameters'] = read_parameters(mad, parameter_dict)
my_workflow_dict['make_sequence'] = make_sequence(mad, parameter_dict['par_mylhcbeam'], parameter_dict['par_slicefactor'])
my_workflow_dict['load_optics'] = load_optics(mad)
my_workflow_dict['check_optics'] = run_module(mad,'module_00_check_optics.madx')
Setting the crossing angles and separations
# This is the relation from DA studies from Nikos and Stéphane.
def from_beta_to_xing_angle_urad(beta_m):
return 0.5*(132.47 + 58.3959 * np.sqrt(beta_m) + 30.0211 * beta_m)/np.sqrt(beta_m)
knob_dict={
'on_sep1': 0,
'on_sep5': 0,
'on_sep2h': 2,
'on_sep2v': 0,
'on_x2h': 0,
'on_x2v': 200,
'on_sep8h': 0,
'on_sep8v': 1,
'on_x8h': 0,
'on_x8v': 135,
'on_disp': 1,
'on_alice': 7000/parameter_dict['par_beam_energy_tot'],
'on_lhcb': 7000/parameter_dict['par_beam_energy_tot'],
'on_sol_atlas': 7000/parameter_dict['par_beam_energy_tot'],
'on_sol_cms': 7000/parameter_dict['par_beam_energy_tot'],
'on_sol_alice': 7000/parameter_dict['par_beam_energy_tot'],
}
betx_ip1 = mad.globals['betx_ip1']
knob_dict['on_x1'] = from_beta_to_xing_angle_urad(betx_ip1)
knob_dict['on_x5'] = from_beta_to_xing_angle_urad(betx_ip1)
for i in knob_dict:
mad.input(f'{i} = {knob_dict[i]};')
mad.input('on_sep8=on_sep8v;')
mad.input('on_sep2=on_sep2h;')
True
Saving the reference CO
my_workflow_dict['save_crossing'] = run_module(mad,'module_01_save_crossing.madx')
execution_df=pd.DataFrame(my_workflow_dict).transpose()
About Luminosity
from madxp import luminosity as lumi
B1=mad.sequence.lhcb1.beam
B2=mad.sequence.lhcb2.beam
#check the frequency
assert B1.freq0==B2.freq0
mad.twiss(sequence='lhcb1'); B1_DF=mt.twiss_df(mad.table.twiss)
mad.twiss(sequence='lhcb2'); B2_DF=mt.twiss_df(mad.table.twiss)
def check_luminosity(B1,B2,B1_DF,B2_DF):
for myIP in ['IP1','IP2', 'IP5', 'IP8']:
B1_IP=B1_DF.loc[myIP.lower()+':1']
B2_IP=B2_DF.loc[myIP.lower()+':1']
aux=lumi.L(f=B1.freq0*1e6, nb=parameter_dict['par_nco_'+myIP],
N1=B1.npart, N2=B2.npart,
energy_tot1=B1.energy, energy_tot2=B2.energy,
deltap_p0_1=B1.sige, deltap_p0_2=B2.sige,
epsilon_x1=B1.exn, epsilon_x2=B2.exn,
epsilon_y1=B1.eyn, epsilon_y2=B2.eyn,
sigma_z1=B1.sigt, sigma_z2=B2.sigt,
beta_x1=B1_IP.betx, beta_x2=B2_IP.betx,
beta_y1=B1_IP.bety, beta_y2=B2_IP.bety,
alpha_x1=B1_IP.alfx, alpha_x2=B2_IP.alfx,
alpha_y1=B1_IP.alfy, alpha_y2=B2_IP.alfy,
dx_1=B1_IP.dx, dx_2=B2_IP.dx,
dpx_1=B1_IP.dpx, dpx_2=B2_IP.dpx,
dy_1=B1_IP.dy, dy_2=B2_IP.dy,
dpy_1=B1_IP.dpy, dpy_2=B2_IP.dpy,
x_1=B1_IP.x, x_2=B2_IP.x,
px_1=B1_IP.px, px_2=B2_IP.px,
y_1=B1_IP.y, y_2=B2_IP.y,
py_1=B1_IP.py, py_2=B2_IP.py, verbose=False)
print(f'Luminosity at {myIP}: {aux} Hz/cm^2')
check_luminosity(B1,B2,B1_DF,B2_DF)
Luminosity at IP1: 2.1084853766379376e+34 Hz/cm^2
Luminosity at IP2: 0.0 Hz/cm^2
Luminosity at IP5: 2.1065742028275344e+34 Hz/cm^2
Luminosity at IP8: 0.0 Hz/cm^2
Info
The IP2/IP8 are separated.
# we are inverting the beams
check_luminosity(B1,B2,B1_DF,B2_DF)
print('After B1/2 inversion')
check_luminosity(B2,B1,B2_DF,B1_DF)
Luminosity at IP1: 2.1084853766379376e+34 Hz/cm^2
Luminosity at IP2: 0.0 Hz/cm^2
Luminosity at IP5: 2.1065742028275344e+34 Hz/cm^2
Luminosity at IP8: 0.0 Hz/cm^2
After B1/2 inversion
Luminosity at IP1: 2.1084853766379376e+34 Hz/cm^2
Luminosity at IP2: 0.0 Hz/cm^2
Luminosity at IP5: 2.1065742028275344e+34 Hz/cm^2
Luminosity at IP8: 0.0 Hz/cm^2
Success
If one invert the B½ you get the same luminosity.
Leveling in intensity
We can "level" the beam intensity to the average luminosity in IP1 and IP5.
from scipy.optimize import least_squares
print('\n==== Offset Levelling ====')
L_target=2e+34
starting_guess=B1.npart
def function_to_minimize(N):
B1_IP=B1_DF.loc['ip1:1']
B2_IP=B2_DF.loc['ip1:1']
L_IP1=lumi.L(f=B1.freq0*1e6, nb=parameter_dict['par_nco_IP1'],
N1=N, N2=N,
energy_tot1=B1.energy, energy_tot2=B2.energy,
deltap_p0_1=B1.sige, deltap_p0_2=B2.sige,
epsilon_x1=B1.exn, epsilon_x2=B2.exn,
epsilon_y1=B1.eyn, epsilon_y2=B2.eyn,
sigma_z1=B1.sigt, sigma_z2=B2.sigt,
beta_x1=B1_IP.betx, beta_x2=B2_IP.betx,
beta_y1=B1_IP.bety, beta_y2=B2_IP.bety,
alpha_x1=B1_IP.alfx, alpha_x2=B2_IP.alfx,
alpha_y1=B1_IP.alfy, alpha_y2=B2_IP.alfy,
dx_1=B1_IP.dx, dx_2=B2_IP.dx,
dpx_1=B1_IP.dpx, dpx_2=B2_IP.dpx,
dy_1=B1_IP.dy, dy_2=B2_IP.dy,
dpy_1=B1_IP.dpy, dpy_2=B2_IP.dpy,
x_1=B1_IP.x, x_2=B2_IP.x,
px_1=B1_IP.px, px_2=B2_IP.px,
y_1=B1_IP.y, y_2=B2_IP.y,
py_1=B1_IP.py, py_2=B2_IP.py, verbose=False)
B1_IP=B1_DF.loc['ip5:1']
B2_IP=B2_DF.loc['ip5:1']
L_IP5=lumi.L(f=B1.freq0*1e6, nb=parameter_dict['par_nco_IP5'],
N1=N, N2=N,
energy_tot1=B1.energy, energy_tot2=B2.energy,
deltap_p0_1=B1.sige, deltap_p0_2=B2.sige,
epsilon_x1=B1.exn, epsilon_x2=B2.exn,
epsilon_y1=B1.eyn, epsilon_y2=B2.eyn,
sigma_z1=B1.sigt, sigma_z2=B2.sigt,
beta_x1=B1_IP.betx, beta_x2=B2_IP.betx,
beta_y1=B1_IP.bety, beta_y2=B2_IP.bety,
alpha_x1=B1_IP.alfx, alpha_x2=B2_IP.alfx,
alpha_y1=B1_IP.alfy, alpha_y2=B2_IP.alfy,
dx_1=B1_IP.dx, dx_2=B2_IP.dx,
dpx_1=B1_IP.dpx, dpx_2=B2_IP.dpx,
dy_1=B1_IP.dy, dy_2=B2_IP.dy,
dpy_1=B1_IP.dpy, dpy_2=B2_IP.dpy,
x_1=B1_IP.x, x_2=B2_IP.x,
px_1=B1_IP.px, px_2=B2_IP.px,
y_1=B1_IP.y, y_2=B2_IP.y,
py_1=B1_IP.py, py_2=B2_IP.py, verbose=False)
return 0.5*(L_IP1+L_IP5)-L_target
aux=least_squares(function_to_minimize, starting_guess)
print(aux)
print(f"\nLuminosity after levelling: {function_to_minimize(aux['x'][0])+L_target} Hz/cm^2")
==== Offset Levelling ====
active_mask: array([0.])
cost: 2.658455991569832e+36
fun: array([2.30584301e+18])
grad: array([8.16213177e+41])
jac: array([[3.5397604e+23]])
message: '`xtol` termination condition is satisfied.'
nfev: 5
njev: 5
optimality: 8.16213177232973e+41
status: 3
success: True
x: array([1.13001999e+11])
Luminosity after levelling: 2.0000000000000001e+34 Hz/cm^2
Then we move from python to MAD-X and we recompute the luminosity.
mad.sequence.lhcb1.beam.npart=aux.x[0]
mad.sequence.lhcb2.beam.npart=aux.x[0]
B1=mad.sequence.lhcb1.beam
B2=mad.sequence.lhcb2.beam
check_luminosity(B1,B2,B1_DF,B2_DF)
Luminosity at IP1: 2.000906831219997e+34 Hz/cm^2
Luminosity at IP2: 0.0 Hz/cm^2
Luminosity at IP5: 1.9990931687800033e+34 Hz/cm^2
Luminosity at IP8: 0.0 Hz/cm^2
(2.0011449627161415+1.9988550372838583)/2
2.0
Success
The "intensity leveling" is working.
Offset in IP8
We can now level the IP8 by separation.
Info
In the original mask before leveling the on_disp is forced to 0 before the leveling (and re-established after the leveling macros). We are NOT forcing it to 0 in the following.
print('\n==== Offset Levelling ====')
L_target=2e+32
starting_guess=1e-5
B1_IP=B1_DF.loc['ip8:1']
B2_IP=B2_DF.loc['ip8:1']
def function_to_minimize(delta):
aux=lumi.L(f=B1.freq0*1e6, nb=parameter_dict['par_nco_IP8'],
N1=B1.npart, N2=B2.npart,
energy_tot1=B1.energy, energy_tot2=B2.energy,
deltap_p0_1=B1.sige, deltap_p0_2=B2.sige,
epsilon_x1=B1.exn, epsilon_x2=B2.exn,
epsilon_y1=B1.eyn, epsilon_y2=B2.eyn,
sigma_z1=B1.sigt, sigma_z2=B2.sigt,
beta_x1=B1_IP.betx, beta_x2=B2_IP.betx,
beta_y1=B1_IP.bety, beta_y2=B2_IP.bety,
alpha_x1=B1_IP.alfx, alpha_x2=B2_IP.alfx,
alpha_y1=B1_IP.alfy, alpha_y2=B2_IP.alfy,
dx_1=B1_IP.dx, dx_2=B2_IP.dx,
dpx_1=B1_IP.dpx, dpx_2=B2_IP.dpx,
dy_1=B1_IP.dy, dy_2=B2_IP.dy,
dpy_1=B1_IP.dpy, dpy_2=B2_IP.dpy,
x_1=B1_IP.x, x_2=B2_IP.x,
px_1=B1_IP.px, px_2=B2_IP.px,
y_1=delta, y_2=-delta,
py_1=B1_IP.py, py_2=B2_IP.py, verbose=False)
return aux-L_target
aux=least_squares(function_to_minimize, starting_guess, bounds=(0, 3e-4))
print(aux)
print(f"\nLuminosity after levelling: {function_to_minimize(aux['x'][0])+L_target} Hz/cm^2")
==== Offset Levelling ====
active_mask: array([0])
cost: 2.9074086151568487e+41
fun: array([7.62549489e+20])
grad: array([-2.21281711e+58])
jac: array([[-2.90186688e+37]])
message: '`xtol` termination condition is satisfied.'
nfev: 9
njev: 9
optimality: 5.688322531645036e+54
status: 3
success: True
x: array([4.29375201e-05])
Luminosity after levelling: 2.0000000000076256e+32 Hz/cm^2
Setting the python result in MAD-X
print(f"BEFORE: on_sep8v={mad.globals['on_sep8v']}")
mad.globals['on_sep8v']=aux.x[0]*1e3
print(f"AFTER: on_sep8v={mad.globals['on_sep8v']}")
BEFORE: on_sep8v=1.0
AFTER: on_sep8v=0.04293752005282778
Sanity check
B1=mad.sequence.lhcb1.beam
B2=mad.sequence.lhcb2.beam
#check the frequency
assert B1.freq0==B2.freq0
mad.twiss(sequence='lhcb1'); B1_DF=mt.twiss_df(mad.table.twiss)
mad.twiss(sequence='lhcb2'); B2_DF=mt.twiss_df(mad.table.twiss)
check_luminosity(B1,B2,B1_DF,B2_DF)
Luminosity at IP1: 2.00108737607358e+34 Hz/cm^2
Luminosity at IP2: 0.0 Hz/cm^2
Luminosity at IP5: 1.9988077575957086e+34 Hz/cm^2
Luminosity at IP8: 1.9751712055946312e+32 Hz/cm^2
B1_IP=B1_DF.loc['ip8:1']
B2_IP=B2_DF.loc['ip8:1']
print(B1_IP['y'])
print(B2_IP['y'])
4.308203834994383e-05
-4.295580603092084e-05
Info
The finite precision of the knob can impact on the final luminosity (and, very marginally, the second order effect on dispersion, optics...) that we get.