diff --git a/README.md b/README.md
index ae71497..f1c4499 100644
--- a/README.md
+++ b/README.md
@@ -86,8 +86,8 @@ By default the jobs will be submitted to HTCondor. If you want to use a differen
       ```python
       myWS.init_study('myStudy')
       ```
-   
-   1. edit the `config.py` file to define the study (e.g. define `.mask` file, add scan parameters, etc...);
+
+   1. edit the `config.py` file to define the parameters you want to scan for your study as well as their range. This also includes adapting the `.mask` file and parametrizing the variables that are to be scanned (use `%` sign to declare these parameters followed by the parameter name defined in `config.py`, e.g. `%SEEDRAN`). Note that PySixDesk will set up a scan making all the possible combinations between defined parameter values (Cartesian product);
    
    1. load definition of study in `config.py` and create/update database:
    
@@ -96,12 +96,14 @@ By default the jobs will be submitted to HTCondor. If you want to use a differen
       myStudy.update_db() # only needed for a new study or when parameters are changed
       ```
 
-   1. prepare and submit pre-processing jobs (e.g. MADX job and sixtrack one turn jobs), and collect results:
+   This initializes a database file. Depending on the type of database you choose in the `config.py` file, it is either an SQLite database stored locally, or, alternatively, a MySQL database hosted by the CERN DBonDemand Service.
+
+   1. prepare and submit pre-processing jobs (i.e. MADX job and sixtrack one turn jobs), and collect results:
    
       ```python
       myStudy.prepare_preprocess_input()
       myStudy.submit(0, 5) # 0 stands for preprocess job, 5 is trial number 
-      myStudy.collect_result(0) # collect results locally
+      myStudy.collect_result(0) # collect results locally once jobs are finished
       ```
 
    1. prepare and submit actual sixtrack jobs, and collect results:
@@ -110,8 +112,10 @@ By default the jobs will be submitted to HTCondor. If you want to use a differen
       myStudy.prepare_sixtrack_input()
       or
       myStudy.prepare_sixtrack_input(True) #True: submit jobs to Boinc
+
       myStudy.submit(1, 5) # 1 stands for sixtrack job, 5 is trial number 
-      myStudy.collect_result(1) # 1 stands for sixtrack job 
+
+      myStudy.collect_result(1) # 1 stands for sixtrack job; collect results once finished
       or
       myStudy.collect_result(1, True) # True: collect results from boinc spool directory
       ```
diff --git a/pysixdesk/lib/machineparams.py b/pysixdesk/lib/machineparams.py
index bd57eb9..1ca36d8 100644
--- a/pysixdesk/lib/machineparams.py
+++ b/pysixdesk/lib/machineparams.py
@@ -1,5 +1,3 @@
-from .utils import merge_dicts
-
 # LHC settings
 LHC = {}
 # TODO: fill these in.
@@ -36,8 +34,8 @@
                 tune_x=62.31,
                 tune_y=60.32,
                 )
-LHC['inj'] = merge_dicts(_LHC_DEF, _LHC_INJ)
-LHC['col'] = merge_dicts(_LHC_DEF, _LHC_COL)
+LHC['inj'] = {**_LHC_DEF, **_LHC_INJ}
+LHC['col'] = {**_LHC_DEF, **_LHC_COL}
 
 # HLLHC settings
 HLLHC = {}
@@ -74,5 +72,5 @@
                   tune_x=62.31,
                   tune_y=60.32,
                   )
-HLLHC['inj'] = merge_dicts(_HLLHC_DEF, _HLLHC_INJ)
-HLLHC['col'] = merge_dicts(_HLLHC_DEF, _HLLHC_COL)
+HLLHC['inj'] = {**_HLLHC_DEF, **_HLLHC_INJ}
+HLLHC['col'] = {**_HLLHC_DEF, **_HLLHC_COL}
diff --git a/pysixdesk/lib/preprocess.py b/pysixdesk/lib/preprocess.py
index e887518..423318b 100755
--- a/pysixdesk/lib/preprocess.py
+++ b/pysixdesk/lib/preprocess.py
@@ -190,7 +190,13 @@ def sixtrack_prep_job(self, fort_cfg, source_prefix=None,
             madx_fc3 = source_prefix / madx_fc3
 
         # concatenate
-        utils.concatenate_files([dest, madx_fc3], output_file)
+        # utils.concatenate_files([dest, madx_fc3], output_file)
+        utils.sandwich(dest,
+                       output_file,
+                       path_prefix=source_prefix,
+                       logger=self._logger)
+        # with open(output_file, 'r') as fp:
+        #     print(''.join(fp.readlines()))
         # if not source.samefile(output_file):
         #     utils.diff(source, output_file, logger=self._logger)
 
@@ -338,6 +344,8 @@ def madx_prep(self, output_file='madx_in'):
         # show diff
         # utils.diff(self.madx_cfg["mask_file"], output_file,
         #            logger=self._logger)
+        # with open(output_file, 'r') as fp:
+        #     print(''.join(fp.readlines()))
 
     def madx_run(self, mask):
         """Runs madx.
@@ -361,6 +369,8 @@ def madx_run(self, mask):
             raise Exception(content)
         else:
             self._logger.info("MADX has completed properly!")
+        # with open('madx_stdout', 'r') as fp:
+        #     print(''.join(fp.readlines()))
 
     def madx_job(self):
         """Controls madx job execution.
@@ -495,5 +505,5 @@ def write_oneturnresult(self):
         where = f"task_id={job.task_id}"
         job_table['status'] = 'incomplete'
         job_table['mtime'] = int(time.time() * 1E7)
-        job.db.update('preprocess_wu'. job_table, where)
+        job.db.update('preprocess_wu', job_table, where)
         raise e
diff --git a/pysixdesk/lib/sixtrack.py b/pysixdesk/lib/sixtrack.py
index bd8eb99..6caa09e 100755
--- a/pysixdesk/lib/sixtrack.py
+++ b/pysixdesk/lib/sixtrack.py
@@ -284,8 +284,34 @@ def sixtrack_prep_job(self, fort_cfg, source_prefix=None,
         if source_prefix is not None:
             madx_fc3 = source_prefix / madx_fc3
 
+        utils.sandwich(dest,
+                       output_file,
+                       path_prefix=source_prefix,
+                       logger=self._logger)
+        # with open(output_file, 'r') as fp:
+        #     lines = fp.readlines()
+        #     print(''.join(lines))
+
+
+        # Ugly hack to generate the particle distribution for use with DIST block
+        if self.fort_cfg['dist_type'] == 'polar':
+            amps = json.loads(fort_cfg['phase_space_var1'])
+            angle = fort_cfg['phase_space_var2']
+            num_particle_pairs = int(fort_cfg['nss'] / 2.)
+            os.system(f'python init_polar_dist.py {amps[0]} {amps[1]} ' +
+                      f'{angle} {num_particle_pairs}')
+        elif self.fort_cfg['dist_type'] == 'cartesian':
+            amps_x = json.loads(fort_cfg['phase_space_var1'])
+            amps_y = json.loads(fort_cfg['phase_space_var2'])
+            num_particle_pairs = int(fort_cfg["nss"] / 2.)
+            os.system(
+                f'python init_cartesian_dist.py {amps_x[0]} {amps_x[1]} ' +
+                f'{amps_y[0]} {amps_y[1]} {num_particle_pairs}')
+        else:
+            raise ValueError('Unknown dist_type')
+
         # concatenate
-        utils.concatenate_files([dest, madx_fc3], output_file)
+        # utils.concatenate_files([dest, madx_fc3], output_file)
 
     def sixtrack_run(self, output_file):
         """Runs sixtrack.
@@ -438,7 +464,8 @@ def boinc_submit(self, job_name):
         """
         # zip all the input files, e.g. fort.3 fort.2 fort.8 fort.16
         input_zip = job_name + '.zip'
-        inputs = ['fort.2', 'fort.3', 'fort.8', 'fort.16'] + self.cr_inputs
+        inputs = ['fort.2', 'fort.3', 'fort.8', 'fort.16', 'fort.34'] + self.cr_inputs
+        inputs += ['input_dist.txt', 'fort.3.aper']
         with zipfile.ZipFile(input_zip, 'w', zipfile.ZIP_DEFLATED) as ziph:
             for infile in inputs:
                 if infile in os.listdir('.'):
diff --git a/pysixdesk/lib/study.py b/pysixdesk/lib/study.py
index 23828fa..3ade157 100644
--- a/pysixdesk/lib/study.py
+++ b/pysixdesk/lib/study.py
@@ -83,7 +83,7 @@ def _defaults(self):
             'fc.34': 'fort.34'}
 
         self.oneturn_sixtrack_input['input'] = dict(self.madx_output)
-        self.sixtrack_input['temp'] = 'fort.3'
+        self.sixtrack_input['fort_file'] = 'fort.3'
         self.sixtrack_output = ['fort.10']
 
         self.db_info['db_type'] = 'sql'
@@ -248,6 +248,7 @@ def _update_db_params(self):
         madx_keys = self.params.madx.keys()
         six_keys = (list(self.params.sixtrack.keys()) +
                     list(self.params.phasespace.keys()))
+
         # preprocess_wu already in the table
         prep_entries = self.db.select('preprocess_wu', madx_keys)
         prep_wu_done = set(prep_entries)
@@ -394,6 +395,9 @@ def _prep_sixtrack_cfg(self):
         self.sixtrack_config['sixtrack']['output_files'] = json.dumps(self.sixtrack_output)
         self.sixtrack_config['sixtrack']['test_turn'] = json.dumps(self.env['test_turn'])
         self.sixtrack_config['six_results'] = self.tables['six_results']
+        # TODO: UGLY HACK TO ADD THE aperture_losses FILE...
+        # self.sixtrack_config['aperture_losses'] = self.tables['aperture_losses']
+
         if self.collimation:
             self.sixtrack_config['aperture_losses'] = self.tables['aperture_losses']
             self.sixtrack_config['collimation_losses'] = self.tables['collimation_losses']
@@ -445,17 +449,18 @@ def _run_calcs(self):
                 self._logger.info('Queued update of sixtack_wu/wu_id:'
                                   f'{row["wu_id"]} with {update_cols}.')
 
-        # update everything at once
-        self._logger.info(f'Updating {len(calc_out_to_be_updated)} rows of '
-                          'sixtrack_wu.')
-        # turn list of dicts into dict of lists
-        calc_out_to_be_updated = {k: [dic[k] for dic in calc_out_to_be_updated]
-                                  for k in calc_out_to_be_updated[0]}
-        where_to_be_updated = {k: [dic[k] for dic in where_to_be_updated]
-                               for k in where_to_be_updated[0]}
-        self.db.updatem('sixtrack_wu',
-                        calc_out_to_be_updated,
-                        where=where_to_be_updated)
+        if calc_out_to_be_updated:
+            # update everything at once
+            self._logger.info(f'Updating {len(calc_out_to_be_updated)} rows of '
+                              'sixtrack_wu.')
+            # turn list of dicts into dict of lists
+            calc_out_to_be_updated = {k: [dic[k] for dic in calc_out_to_be_updated]
+                                      for k in calc_out_to_be_updated[0]}
+            where_to_be_updated = {k: [dic[k] for dic in where_to_be_updated]
+                                   for k in where_to_be_updated[0]}
+            self.db.updatem('sixtrack_wu',
+                            calc_out_to_be_updated,
+                            where=where_to_be_updated)
 
     def update_db(self, db_check=False):
         '''Update the database whith the user-defined parameters'''
@@ -545,11 +550,12 @@ def query(index):
                 print(query_list)
                 for i in wus:
                     print(i)
+            return results
 
         if job == 0 or job == 2:
-            query(0)
+            return query(0)
         if job == 1 or job == 2:
-            query(1)
+            return query(1)
 
     def submit(self, typ, trials=5, *args, **kwargs):
         '''Sumbit the preporcess or sixtrack jobs to htctondor.
@@ -627,6 +633,8 @@ def collect_result(self, typ, boinc=False):
             info_sec['boinc_results'] = self.env['boinc_results']
             info_sec['boinc'] = boinc
             info_sec['st_pre'] = self.st_pre
+            # TODO: UGLY HACK TO PARSE THE APERTURE_LOSSES FILE...
+            # config['aperture_losses'] = self.tables['aperture_losses']
 
             info_sec['outs'] = Table.result_table(self.sixtrack_output)
             if self.collimation:
@@ -670,7 +678,9 @@ def prepare_sixtrack_input(self, resubmit=False, boinc=False, groupby=None,
             action = 'resubmit'
         else:
             constraints = "status='incomplete' and preprocess_id in %s" % str(
-                preprocess_outs[0])
+                preprocess_outs[0]).replace(',)', ')')  # if there is a single
+            # job, the comma of the tuple (1,) breaks the sql query because
+            # there is no proper sql sanitization...
             action = 'submit'
         names = self.tables['sixtrack_wu'].keys()
         outputs = self.db.select('sixtrack_wu',
diff --git a/pysixdesk/lib/study_params.py b/pysixdesk/lib/study_params.py
index 30cbcff..a82b88e 100644
--- a/pysixdesk/lib/study_params.py
+++ b/pysixdesk/lib/study_params.py
@@ -1,6 +1,7 @@
 import re
 import logging
 
+from math import pi
 from pathlib import Path
 from collections import OrderedDict
 from itertools import product
@@ -9,7 +10,7 @@
 
 from . import machineparams
 from .constants import PROTON_MASS
-from .utils import PYSIXDESK_ABSPATH, merge_dicts
+from .utils import PYSIXDESK_ABSPATH, linspace
 
 
 class StudyParams:
@@ -46,7 +47,7 @@ def __init__(self,
         # comment regexp
         self._reg_comment = re.compile(r'^(\s?!|\s?/).*', re.MULTILINE)
         # placeholder pattern regexp
-        self._reg = re.compile(r'%([a-zA-Z0-9_]+/?)')
+        self._reg = re.compile(r'%(?!FILE|%)([a-zA-Z0-9_]+/?)')
         self.fort_path = fort_path
         self.mask_path = mask_path
         # initialize empty calculation queue
@@ -85,18 +86,17 @@ def __init__(self,
                             "writebins": 1,
                             }
         self.machine_defaults = machine_defaults
-        self.defaults = merge_dicts(self.f3_defaults, self.machine_defaults)
+        self.defaults = {**self.f3_defaults, **self.machine_defaults}
         # phasespace params
         # TODO: find sensible defaults for the phasespace parameters.
-        amp = [8, 10, 12]  # The amplitude
-        self.phasespace = {"amp": list(zip(amp, amp[1:])),
-                           "kang": list(range(1, 1 + 1)),
-                           "kmax": 5,
+        self.phasespace = {"phase_space_var1": [],
+                           "phase_space_var2": [],
                            }
 
         self.madx = self.find_patterns(self.mask_path)
         self.sixtrack = self.find_patterns(self.fort_path,
                                            mandatory=['chrom_eps', 'CHROM'])
+        self.sixtrack['dist_type'] = 'None'
 
     @property
     def _sixtrack_only(self):
@@ -117,6 +117,10 @@ def oneturn(self):
         sixtrack['toggle_coll/'] = '/'
         return sixtrack
 
+    @staticmethod
+    def da_angles(start=0, end=pi/2, n=7):
+            return linspace(start, end, n + 2)[1: -1]  # exclusive
+
     def keys(self):
         """Gets the keys of 'self.madx', 'self.sixtrack' and 'self.phasespace'.
 
@@ -180,7 +184,7 @@ def find_patterns(self, file_path, keep_none=True, mandatory=None):
 
     @staticmethod
     def _combinations_prep(**kwargs):
-        '''Sanitizes the paramter values.
+        '''Sanitizes the parameter values.
 
         Args:
             **kwargs: Parameter name = parameter value.
diff --git a/pysixdesk/lib/utils.py b/pysixdesk/lib/utils.py
index 27afbd9..406ec76 100644
--- a/pysixdesk/lib/utils.py
+++ b/pysixdesk/lib/utils.py
@@ -59,6 +59,58 @@ def download_output(filenames, dest, zp=True):
                 shutil.copy(filename, dest)
 
 
+def sandwich(in_file, out_file, path_prefix='', logger=None):
+    '''
+    Looks for any patterns matching '^%FILE:.*' in in_file, then replaces the
+    match with the content of the file following the match.
+    A path prefix can be specified to look for matched file in another
+    directory. If the matched file is not found, comment out the pattern.
+    Example:
+        contents of in_file:
+            aaaaaaaaa
+            %FILE:insert.txt
+            aaaaaaaaa
+        contents of insert.txt:
+            bbbbbbbbb
+            bbbbbbbbb
+        writes to out_file:
+            aaaaaaaaa
+            bbbbbbbbb
+            bbbbbbbbb
+            aaaaaaaaa
+    TODO: maybe it's best to return the the sandwiched lines and not write the
+    file, so that more parsing can happen in memory without reopening and
+    writing files. Also might be best to take as input the contents itself for
+    the same reason.
+    '''
+
+    if logger is not None:
+        display = logger.warning
+    else:
+        display = print
+
+    with open(in_file, 'r') as f:
+        in_lines = f.read()
+
+    reg = re.compile(r'^%FILE:.*', re.MULTILINE)
+    for m in re.finditer(reg, in_lines):
+        m_str = m.group()
+        try:
+            if path_prefix is None:
+                path_prefix = ''
+            fname = os.path.join(path_prefix, m_str.split(':')[1].lstrip())
+            with open(fname, 'r') as f:
+                file_lines = f.read()
+            in_lines = re.sub(f'{m_str}', file_lines, in_lines)
+        except FileNotFoundError as e:
+            display(e)
+            display(f'Commenting out {m_str} for {out_file}')
+            in_lines = re.sub(f'{m_str}', f'/{m_str}', in_lines)
+
+    with open(out_file, 'w') as out:
+        out.write(in_lines)
+
+
 def check_fort3_block(fort3, block):
     '''Check the existence of the given block in fort.3'''
 
@@ -238,15 +290,13 @@ def condor_logger(name):
     return logger
 
 
-def merge_dicts(x, y):
-    """Merges two dicts.
-
-    Returns:
-        dict: Merged dict
-    """
-    z = x.copy()
-    z.update(y)
-    return z
+def linspace(a, b, n):
+    '''Numpyless linear spacing function.
+    '''
+    if n < 2:
+        return a
+    diff = (float(b) - a)/(n - 1)
+    return [diff * i + a for i in range(n)]
 
 
 class ProgressBar(object):
diff --git a/templates/config.py b/templates/config.py
index fc3134a..97f26a1 100644
--- a/templates/config.py
+++ b/templates/config.py
@@ -4,6 +4,8 @@
 import os
 import logging
 
+from collections import OrderedDict
+
 from pysixdesk.lib import submission
 from pysixdesk import Study
 from pysixdesk.lib.study_params import StudyParams
@@ -12,6 +14,7 @@
 from pysixdesk.lib.study_params import set_output_keys
 from math import sqrt, pi, sin, cos, tan
 from pysixdesk.lib import machineparams
+from pysixdesk.lib import utils
 
 # logger configuration
 logger = logging.getLogger('pysixdesk')
@@ -39,65 +42,120 @@ def __init__(self, name='study', location=os.getcwd()):
         self.paths['boinc_spool'] = '/afs/cern.ch/work/b/boinc/boinctest'
         self.boinc_vars['appName'] = 'sixtracktest'
 
+        # Need to use more recent release of Sixtrack than what is defined in PySixDesk by default.
+        self.paths['sixtrack_exe'] = ('/afs/cern.ch/project/sixtrack/build/latest/' + 
+            'SixTrack_50403-9e43f07_zlib_crlibm_rn_Linux_gfortran_static_avx2_x86_64_64bit_double')
+
         # Database type
         self.db_info['db_type'] = 'sql'
         # self.db_info['db_type'] = 'mysql'
 
-        self.oneturn = True  # Switch for oneturn sixtrack job
+        self.oneturn = False  # Switch for oneturn sixtrack job
         # self.collimation = True
 
         # All parameters are case-sensitive
         # the name of mask file
-        mask_file = 'hl10.mask'
+        mask_file = 'example.mask'
         self.madx_input["mask_file"] = mask_file
+        self.madx_output = {'fc.2': 'fort.2',
+                            'fc.3': 'fort.3.mad',
+                            'fc.3.aux': 'fort.3.aux',
+                            'fc.8': 'fort.8',
+                            'fc.16': 'fort.16',
+                            'fc.34': 'fort.34'}
+
         # All parameters are case-sensitive
         self.params = StudyParams(mask_path=os.path.join(self.study_path, mask_file),
                                   fort_path=os.path.join(self.study_path, 'fort.3'),
-                                  machine_defaults=machineparams.HLLHC['col'])
+                                  machine_defaults=machineparams.LHC['inj'])
         self.params['Runnam'] = name
-        amp = [8, 10, 12]
-        self.params['turnss'] = int(1e2)  # number of turns to track (must be int)
-        self.params['nss'] = 30
-        self.params['amp'] = list(zip(amp, amp[1:]))  # Take pairs
-        self.params['kang'] = list(range(1, 1 + 2))  # The angle
-        self.params['kmax'] = 5
-        # self.params['toggle_coll/'] = '/'
-        self.params['seed_ran'] = [1, 2]
-        self.params['i_mo'] = list(range(100, 200 + 1, 100))
-
-        @set_input_keys(['kang', 'kmax'])
-        @set_output_keys(['angle'])
-        def calc_angle(kang, kmax):
-            return kang / (kmax + 1)  # 1-->pi/2, 0.5-->pi/4, ...
-        self.params.calc_queue.append(calc_angle)
-
-        @set_input_keys(['angle'])
-        @set_output_keys(['ratio'])
-        def calc_ratio(angle):
-            ratio = abs(tan((pi / 2) * angle))
-            if ratio < 1e-15:
-                ratio = 0.
-            else:
-                ratio = ratio ** 2
-            return ratio
-        self.params.calc_queue.append(calc_ratio)
-
-        # should it not be betax and betax2 ?
-        @set_requirements({'oneturn_sixtrack_results': ['betax', 'betax2']})
-        @set_input_keys(['angle', 'ratio', 'emit_norm_x', 'e_0', 'pmass', 'amp'])
-        @set_output_keys(['ax0s', 'ax1s'])
-        def calc_amp(angle, ratio, emit, e_0, pmass, amp, betax=None, betax2=None):
-            gamma = e_0 / pmass
-            factor = sqrt(emit / gamma)
-            ax0t = factor * (sqrt(betax) + sqrt(betax2 * ratio) * cos((pi / 2) * angle))
-            return amp[0] * ax0t, amp[1] * ax0t
-        self.params.calc_queue.append(calc_amp)
+        self.params['turnss'] = int(500)  # number of turns to track (must be int)
+
+
+        # Particle initialization
+        self.params['dist_type'] = 'cartesian'  # type of particle distribution (cartesian or polar)
+        amp = [float(i) for i in range(5, 16, 1)]
+
+        if self.params['dist_type'] == 'polar':
+            self.params['phase_space_var1'] = list(zip(amp, amp[1:]))  # Take pairs # amp_x
+            self.params['phase_space_var2'] = utils.linspace(0, pi/2., 12)[1:-1] # angles
+            n_particles_per_amplitude_angle_pair = 10
+            self.params['nss'] = 2 * n_particles_per_amplitude_angle_pair  # account for twin particles
+        elif self.params['dist_type'] == 'cartesian':
+            self.params['phase_space_var1'] = list(zip(amp, amp[1:]))  # Take pairs # amp_x
+            self.params['phase_space_var2'] = list(zip(amp, amp[1:]))  # Take pairs # amp_y
+            n_particles_per_amp_pair = 10
+            self.params['nss'] = 2 * n_particles_per_amp_pair**2  # account for twin particles
+        else:
+            raise ValueError('Unknown dist_type')
+
+        self.params['seed_ran'] = 10
+        self.params['sixtrack_seed'] = 42
+
+        self.params['i_mo'] = -10.
+        self.params['tune_x'] = 62.28
+        self.params['tune_y'] = 60.31
+        self.params['int_tune_x'] = 62
+        self.params['int_tune_y'] = 60
+
+        q_prime = 20.
+        self.params['q_prime'] = q_prime
+        self.params['q_prime_x'] = q_prime
+        self.params['q_prime_y'] = q_prime
+
+        self.params['sig_e'] = 4.5e-4
+        self.params['sig_z'] = 0.11
+        self.params['bunch_charge'] = 1e11
+        self.params['emit_norm'] = 3.5
+        self.params['b_t_dist'] = 25.
+        self.params['e_0'] = 450000.
+        self.params['rf_vol'] = 8.
+
+        # @set_input_keys(['angle', 'amp'])
+        # @set_output_keys(['x0', 'x1', 'y0', 'y1'])
+        # def calc_dist_coords(angle, amp):
+        #     x0, x1 = cos(angle)*amp[0], cos(angle)*amp[1]
+        #     y0, y1 = sin(angle)*amp[0], sin(angle)*amp[1]
+        #     return x0, x1, y0, y1
+        # self.params.calc_queue.append(calc_dist_coords)
+
+        #     return angles
+        # self.params.calc_queue.append(set_angles)
+
+        # @set_input_keys(['kang', 'kmax'])
+        # @set_output_keys(['angle'])
+        # def calc_angle(kang, kmax):
+        #     return kang / (kmax + 1)  # 1-->pi/2, 0.5-->pi/4, ...
+        # self.params.calc_queue.append(calc_angle)
+
+        # @set_input_keys(['angle'])
+        # @set_output_keys(['ratio'])
+        # def calc_ratio(angle):
+        #     ratio = abs(tan((pi / 2) * angle))
+        #     if ratio < 1e-15:
+        #         ratio = 0.
+        #     else:
+        #         ratio = ratio ** 2
+        #     return ratio
+        # self.params.calc_queue.append(calc_ratio)
+
+        # # should it not be betax and betax2 ?
+        # @set_requirements({'oneturn_sixtrack_results': ['betax', 'betax2']})
+        # @set_input_keys(['angle', 'ratio', 'emit_norm_x', 'e_0', 'pmass', 'amp'])
+        # @set_output_keys(['ax0s', 'ax1s'])
+        # def calc_amp(angle, ratio, emit, e_0, pmass, amp, betax=None, betax2=None):
+        #     gamma = e_0 / pmass
+        #     factor = sqrt(emit / gamma)
+        #     ax0t = factor * (sqrt(betax) + sqrt(betax2 * ratio) * cos((pi / 2) * angle))
+        #     return amp[0] * ax0t, amp[1] * ax0t
+        # self.params.calc_queue.append(calc_amp)
 
         self.oneturn_sixtrack_input['fort_file'] = 'fort.3'
         self.oneturn_sixtrack_output = ['oneturnresult']
         self.sixtrack_input['fort_file'] = 'fort.3'
         self.preprocess_output = dict(self.madx_output)
         self.sixtrack_input['input'] = dict(self.preprocess_output)
+        self.sixtrack_input['additional_input'] = ['init_cartesian_dist.py', 'init_polar_dist.py']
 
         # # The parameters for collimation job
         # self.madx_output = {
diff --git a/templates/example.mask b/templates/example.mask
new file mode 100644
index 0000000..14de466
--- /dev/null
+++ b/templates/example.mask
@@ -0,0 +1,643 @@
+! Links definitions
+!option, warn,info;
+!###################################################################################################################################
+!############################## CREATE SOFT LINKS AND DIRECTORIES ##################################################################
+!###################################################################################################################################
+System,"rm -r temp";
+System,"rm -r output";
+System,"mkdir temp";
+System,"mkdir output";
+!######################################################################
+system,"ln -fns /afs/cern.ch/eng/lhc/optics/V6.503 db5";
+system,"ln -fns /afs/cern.ch/eng/lhc/optics/HLLHCV1.2 slhc";
+
+system,"ln -fns /afs/cern.ch/eng/lhc/optics/V6.503/FiDeL fidel";
+system,"ln -fns /afs/cern.ch/eng/lhc/optics/V6.503/WISE wise";
+system,"ln -fns /afs/cern.ch/eng/lhc/optics/runII/2016 optics2016";
+system,"ln -fns /afs/cern.ch/eng/lhc/optics/runII/2017 optics2017";
+system,"ln -fns /afs/cern.ch/eng/lhc/optics/runII/2018 optics2018";
+
+!######################################################################
+
+!######################################################################
+
+
+! --------------- TUNE AND CHROMA ----------------
+!********************************************
+
+NRJ= %e_0 * 1e-3;
+
+! A Landau octupole current 20A inj, -570A col
+I_MO= %i_mo;
+Imax_MO = 550.;
+Kmax_MO = 378600.;
+brho:=NRJ*1e9/clight;
+value, brho;
+value, clight;
+value, Imax_MO;
+MOpowering=Kmax_MO*I_MO/Imax_MO/brho;
+MOFpowering=MOpowering;
+MODpowering=MOpowering;
+value, MOpowering;
+
+b_t_dist := %b_t_dist; !25 bunch separation [ns]
+emit_norm_x := %emit_norm * 1e-6; ! [um] --> [m]
+emit_norm_y := %emit_norm * 1e-6; ! [um] --> [m]
+Nb_0 := %bunch_charge;
+sigt := %sig_z; ! bunch length [m] in collision
+sige := %sig_e;
+gamma_rel := NRJ/pmass;
+epsx := emit_norm_x /gamma_rel;
+epsy := emit_norm_y /gamma_rel;
+
+!********************************************
+
+qx00 = %int_tune_x;  !integer h tune
+qy00 = %int_tune_y;  !integer v tune
+tsplit = qx00 - qy00;   !new tune split
+qx0 = %tune_x;
+qy0 = %tune_y;
+qprimex = %q_prime;
+qprimey = %q_prime;
+
+qxinit = qx0;
+qyinit = qy0;
+ptcqxinit=qxinit-qx00;
+ptcqyinit=qyinit-qy00;
+qxinitint=qxinit-ptcqxinit;
+qyinitint=qyinit-ptcqyinit; ! WATCH OUT - ALL THESE VALUES ARE USED BY MATCHING ROUTINES AND THE DQMIN CALCULATION LATER ON IN THE SCIPT...
+dqxinit = qprimex;
+dqyinit = qprimey;
+
+!some values required by new correction macros
+tsplit= abs(qxinitint-qyinitint);   !new tune split
+qx00  = qxinitint;  !integer h tune
+qy00  = qyinitint;  !integer v tune
+
+value, qxinit,qyinit;
+value, qx00,qy00;
+value, tsplit;
+value, ptcqxinit,ptcqyinit;
+
+!######### CALL MACRO FACILITIES ###########################################
+
+call,file="slhc/toolkit/macro.madx";
+
+!###################################################################################################################################
+!!!!POINTS TO CHECK WHEN EDITING THE SCRIPT
+!###################################################################################################################################
+!! OPTIC: set optics below in the set up lattice section, just pick the .str file desired
+!!
+!! ERRORS: Magnetic errors are applied piecewise to keep orbit and coupling under control. Should ensure when changing seed/script
+!! that all relevant parts are varied. This is mostly redundancy to ensure all seeds and configs run. Can be reduced if confident
+!! the mask will run.
+!! LINEAR OPTICS: Have assumed correction of a2 and b2 errors in the IR magnets has already been applied at both injection and top
+!! energy. Block is there, but commented out.
+!! CORRS: a2 is compensated effectively using tracking data and the OMC codes. This is to allow setting of a well defined amp and
+!! phase of f1001 at the end of the mask, reflecting operational procedure for optics MDs.
+!! CORRS: b2 is corrected by application nominal corrections determined for the WISE zero uncertainty seed in the thin lattice.
+!! Probably underestimates the true beta-beating, but okayish for now.
+!! CORRS: b3. b3 at injection is compensated by application of nominal correction determined for the WISE zero uncertainty seed in
+!! the thin lattice, with the bad arc depowered and the correction redistributed to the other arcs as per operation.
+!! CORRS: b3. b3 at top energy is compensated by applying the MCS settings found in LSA.
+!! CORRS: b3. Note that at injection the b3 is complicated by the decay. Should be possible in the future to use timestamped wise
+!! predictions, alongside the real machine settings, but both have to changed consistently.
+!! CORRS: 2016 -> Spool piece settings are NOT included in the normal beam-processes! They have a separate BP. If you take spools
+!! from the regular beam-process you will be off by up to an order of magnitude.
+!! XING SCHEME: settings from 2016 will not be valid with 2012 sequences etc (ie definition is totally different, for some Runi
+!! sequences don't set values in urad, set fractions of nominal). Also with Runi sequences watch out for changing definitions of the
+!! xing scheme between strength files for different beams and optics. Always check explicitly with a twiss.
+!! XING SCHEME: Make sure internal xings from alice spectrometer are what is desired. This is set explicitly later on in the script.
+!! ORBIT: if want virtual correctors to match to zero or to target set explicitly towards end of script
+!! ORBIT: running with virtual correctors to match orbit will make this SLOW! for tracking. If don't care too much about orbit can
+!! match using old method, just check the RMS and P2P at the end to make sure it is sensible.
+!! MO: script assumes uniform powering (with exception of ROD.A56.B1) consistent with Landau damping knob. For Asymettric powering
+!! marco needs to be edited later on.
+!! BEAMBASED: Effective models for MCD feeddown based on injection measurements, and measured detuning at 40cm are implemented later
+!! on in the scipt.
+!!
+!! MATCHING: at end of script switch to PTC based tune matching. This is slow...
+
+
+!###################################################################################################################################
+!############################## OPTIONS ############################################################################################
+!###################################################################################################################################
+mylhcbeam=1;
+
+! Set this flag if the file is not used as a mask file (sets %SEEDRAN to 1)
+not_a_mask = 0;
+is_thin=1;
+
+!********************************************
+
+if(mylhcbeam<3){ call,file="optics2018/lhc_as-built.seq";   };
+if(mylhcbeam>3){ call,file="optics2018/lhcb4_as-built.seq"; };
+
+Option, -echo,-warn,-info;
+
+!###################################################################################################################################
+!############################## SET UP LATTICE 2017 ################################################################################
+!###################################################################################################################################
+
+!if(is_thin==1) { call,file="optics2018/toolkit/slice.madx"; };
+if(is_thin==1) { slicefactor=8;
+           beam;
+           call,file="optics2018/toolkit/myslice.madx";
+     beam;
+     use, sequence=lhcb1; makethin, sequence=lhcb1, style=teapot, makedipedge=true;
+     use, sequence=lhcb2; makethin, sequence=lhcb2, style=teapot, makedipedge=true;
+
+};
+
+!Cycling w.r.t. to IP3 (mandatory to find closed orbit in collision in the presence of errors)
+if(mylhcbeam==1) {seqedit,sequence=lhcb1;flatten;cycle,start=IP3;endedit;}
+if(mylhcbeam>1)  {seqedit,sequence=lhcb2;flatten;cycle,start=IP3;endedit;}
+
+
+if(is_thin==0){call, file="optics2018/PROTON/opticsfile.1";
+}
+else {         call, file="optics2018/PROTON/opticsfile.1";
+}
+
+NRJ = %e_0 * 1e-3;
+if (mylhcbeam<3){Beam, particle=proton, sequence=lhcb1, energy=NRJ, sigt=sigt, bv=1, sige=sige, NPART=Nb_0, kbunch=1, ex=epsx, ey=epsy;
+          Beam, particle=proton, sequence=lhcb2, energy=NRJ, bv=-1, sigt=sigt, NPART=Nb_0, sige=sige, kbunch=1, ex=epsx, ey=epsy;};
+if (mylhcbeam>3){Beam, particle=proton, sequence=lhcb2, energy=NRJ, sigt=sigt, bv=1, NPART=Nb_0, sige=sige, ex=epsx, ey=epsy;};
+
+if(mylhcbeam==1) {Use, sequence=lhcb1;}
+if(mylhcbeam>1) {Use, sequence=lhcb2;}
+
+
+ select,flag=twiss,clear;
+ select, flag=twiss, pattern="IP", column=name,s,betx,bety,x,y,px,py;
+ twiss, chrom, sequence=lhcb1, table, file=twiss.final.b1.dat;
+
+!###################################################################################################################################
+!################################ SETUP FOR ORBIT and ERROR CORRECTIONS ############################################################
+!###################################################################################################################################
+on_x1= 0    ; on_sep1= 0 ; on_o1= 0 ;
+on_x2= 0    ; on_sep2= 0 ; on_o2= 0 ; on_oe2=0; on_a2= 0 ;
+on_x5= 0    ; on_sep5= 0 ; on_o5= 0 ;
+on_x8= 0    ; on_sep8= 0 ; on_o8= 0 ; on_a8= 0 ;  on_sep8h= 0 ; on_x8v= 0 ;
+on_alice= 0 ;
+on_sol_alice=0;
+on_lhcb = 0 ;
+on_sol_atlas=0;
+on_sol_cms=0;
+
+  on_x1  :=  0.000  ;
+  on_sep1:=   0  ;
+  on_oh1 :=    0.000  ;
+  on_ov1 :=    0.000  ;
+  phi_IR1:=   90.000  ;
+  on_ssep1:=on_sep1;
+  on_xx1  := on_x1 ;
+  on_x5  :=  0.000  ;
+  on_sep5:=    0.0  ;
+  on_oh5 :=    0.000  ;
+  on_ov5 :=   0.000  ;
+  phi_IR5:=    0.000  ;
+  on_ssep5:=on_sep5;
+  on_xx5  := on_x5 ;
+  on_x2  :=  0.000  ;
+  on_sep2:=    0.000  ;
+  on_a2  :=    0.000  ;
+  on_o2  :=    0.000  ;
+  on_oh2 :=    0.000  ;
+  on_ov2 :=   0.000  ;
+  phi_IR2:=   90.000  ;
+  on_x8  := 0.000  ;
+  on_sep8:=   0.000  ;
+  on_a8  :=    0.000  ;
+  on_o8  :=    0.000  ;
+  on_oh8 :=    0.000  ;
+  on_ov8 :=    0.000  ;
+  phi_IR8:=  180.000  ;
+
+
+select,flag=twiss,clear;
+if(mylhcbeam==1) {twiss,sequence=lhcb1,table=nominal; };
+if(mylhcbeam>1) {twiss,sequence=lhcb2,table=nominal; };
+
+ beta.ip1=table(twiss,IP1,betx);
+
+! print nominal optics parameter at the MB, MQS and MSS for
+! b2, b3, b4, b5, a2 and a3 correction
+ select, flag=twiss, clear;
+ select, flag=twiss,pattern=MB\.   ,class=multipole,column=name,k0L,k1L,betx,bety,dx,mux,muy;
+ select, flag=twiss,pattern=MQT\.14,class=multipole,column=name,k0L,k1L,betx,bety,dx,mux,muy;
+ select, flag=twiss,pattern=MQT\.15,class=multipole,column=name,k0L,k1L,betx,bety,dx,mux,muy;
+ select, flag=twiss,pattern=MQT\.16,class=multipole,column=name,k0L,k1L,betx,bety,dx,mux,muy;
+ select, flag=twiss,pattern=MQT\.17,class=multipole,column=name,k0L,k1L,betx,bety,dx,mux,muy;
+ select, flag=twiss,pattern=MQT\.18,class=multipole,column=name,k0L,k1L,betx,bety,dx,mux,muy;
+ select, flag=twiss,pattern=MQT\.19,class=multipole,column=name,k0L,k1L,betx,bety,dx,mux,muy;
+ select, flag=twiss,pattern=MQT\.20,class=multipole,column=name,k0L,k1L,betx,bety,dx,mux,muy;
+ select, flag=twiss,pattern=MQT\.21,class=multipole,column=name,k0L,k1L,betx,bety,dx,mux,muy;
+ select, flag=twiss,class=MQS                      ,column=name,k0L,k1L,betx,bety,dx,mux,muy;
+ select, flag=twiss,class=MSS                      ,column=name,k0L,k1L,betx,bety,dx,mux,muy;
+ select, flag=twiss,class=MCO                      ,column=name,k0L,k1L,betx,bety,dx,mux,muy;
+ select, flag=twiss,class=MCD                      ,column=name,k0L,k1L,betx,bety,dx,mux,muy;
+ select, flag=twiss,class=MCS                      ,column=name,k0L,k1L,betx,bety,dx,mux,muy;
+ twiss,  file='temp/optics0_MB.mad';
+
+! print nominal optics parameter at the D1, MQX and triplet corrector
+! for triplet correction
+ select, flag=twiss, clear;
+ select, flag=twiss, pattern=MQX  , class=multipole, column=name,betx,bety,x,y;
+ select, flag=twiss, pattern=MBX  , class=multipole, column=name,betx,bety,x,y;
+ select, flag=twiss, pattern=MBRC , class=multipole, column=name,betx,bety,x,y;
+ select, flag=twiss, pattern=MQSX                  , column=name,betx,bety,x,y;
+ select, flag=twiss, pattern=MCSX                  , column=name,betx,bety,x,y;
+ select, flag=twiss, pattern=MCTX                  , column=name,betx,bety,x,y;
+ select, flag=twiss, pattern=MCOSX                 , column=name,betx,bety,x,y;
+ select, flag=twiss, pattern=MCOX                  , column=name,betx,bety,x,y;
+ select, flag=twiss, pattern=MCSSX                 , column=name,betx,bety,x,y;
+ twiss,  file='temp/optics0_inser.mad';
+
+
+Use, sequence=lhcb1;
+select,flag=twiss,clear;
+select, flag=twiss, pattern="IP1",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
+select, flag=twiss, pattern="IP2",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
+select, flag=twiss, pattern="IP5",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
+select, flag=twiss, pattern="IP8",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
+twiss,sequence=lhcb1,file=twiss.b1;system,"cat twiss.b1";
+
+Use, sequence=lhcb2;
+select,flag=twiss,clear;
+select, flag=twiss, pattern="IP1",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
+select, flag=twiss, pattern="IP2",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
+select, flag=twiss, pattern="IP5",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
+select, flag=twiss, pattern="IP8",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
+twiss, sequence=lhcb2,file=twiss.b2;system,"cat twiss.b2";
+
+Use, sequence=lhcb1;
+select,flag=twiss,clear;
+select, flag=twiss, pattern="IP1",   column=name,s,x,y,px,py;
+select, flag=twiss, pattern="IP2",   column=name,s,x,y,px,py;
+select, flag=twiss, pattern="IP5",   column=name,s,x,y,px,py;
+select, flag=twiss, pattern="IP8",   column=name,s,x,y,px,py;
+twiss, sequence=lhcb1, file=twiss.b1;system,"cat twiss.b1";
+
+Use, sequence=lhcb2;
+select,flag=twiss,clear;
+select, flag=twiss, pattern="IP1",   column=name,s,x,y,px,py;
+select, flag=twiss, pattern="IP2",   column=name,s,x,y,px,py;
+select, flag=twiss, pattern="IP5",   column=name,s,x,y,px,py;
+select, flag=twiss, pattern="IP8",   column=name,s,x,y,px,py;
+twiss, sequence=lhcb2, file=twiss.b2;system,"cat twiss.b2";
+
+if(mylhcbeam==1) {Use, sequence=lhcb1;}
+if(mylhcbeam>1) {Use, sequence=lhcb2;}
+
+
+if(mylhcbeam==1){
+   match;
+     global, q1=qxinit, q2=qyinit;
+     global, dq1=dqxinit, dq2=dqyinit;
+     vary,   name=dQx.b1, step=1.0E-7 ;
+     vary,   name=dQy.b1, step=1.0E-7 ;
+     vary,   name=dQpx.b1, step=1.0E-7 ;
+     vary,   name=dQpy.b1, step=1.0E-7 ;
+     lmdif,  calls=100, tolerance=1.0E-21;
+   endmatch;
+}
+if(mylhcbeam>1){
+   match;
+     global, q1=qxinit, q2=qyinit;
+     global, dq1=dqxinit, dq2=dqyinit;
+     vary,   name=dQx.b2, step=1.0E-7 ;
+     vary,   name=dQy.b2, step=1.0E-7 ;
+     vary,   name=dQpx.b2, step=1.0E-7 ;
+     vary,   name=dQpy.b2, step=1.0E-7 ;
+     lmdif,  calls=100, tolerance=1.0E-21;
+   endmatch;
+}
+
+
+!###################################################################################################################################
+!################################ CALL ERROR SUBROUTINES ###########################################################################
+!###################################################################################################################################
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+!         Call error subroutines (nominal machine and new IT/D1)
+!                         and error tables
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+print, text="call errror subroutines";
+call,file="optics2016/measured_errors/Msubroutines_new.madx"; ! think the new subroutines are only relevant for MSS - not used pre-2017 so shouldnt make a difference compared to old Msubroutines...
+call,file="optics2016/measured_errors/Msubroutines_MS_MSS_MO_new.madx";
+call,file="optics2016/toolkit/Orbit_Routines.madx";
+call,file="optics2016/toolkit/SelectLHCMonCor.madx";
+readtable, file="optics2016/measured_errors/rotations_Q2_integral.tab";
+call,file="optics2016/errors/macro_error.madx";  ! some macros for error generation
+print, text="error subroutines finished";
+
+
+
+
+!###################################################################################################################################
+!################################ APPLY THE MAGNETIC ERRORS ########################################################################
+!###################################################################################################################################
+!###################################################################################################################################
+!###################################################################################################################################
+
+
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+!                Switch ON/OFF some multipole
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+exec, ON_ALL_MULT;
+ON_LSF:=1;            ! to be set to 0 if run is performed on BOINC. ! double check if this is still true - and if wise actually includes any of these errors
+
+ ON_A1s  =  0 ; ON_A1r  =  0 ; ON_B1s  =  0 ; ON_B1r  =  0 ;
+ ON_A2s  =  0 ; ON_A2r  =  0 ; ON_B2s  =  0 ; ON_B2r  =  0 ;
+ ON_A3s  =  1 ; ON_A3r  =  1 ; ON_B3s  =  1 ; ON_B3r  =  1 ;
+ ON_A4s  =  1 ; ON_A4r  =  1 ; ON_B4s  =  1 ; ON_B4r  =  1 ;
+ ON_A5s  =  1 ; ON_A5r  =  1 ; ON_B5s  =  1 ; ON_B5r  =  1 ;
+ ON_A6s  =  1 ; ON_A6r  =  1 ; ON_B6s  =  1 ; ON_B6r  =  1 ;
+ ON_A7s  =  1 ; ON_A7r  =  1 ; ON_B7s  =  1 ; ON_B7r  =  1 ;
+ ON_A8s  =  1 ; ON_A8r  =  1 ; ON_B8s  =  1 ; ON_B8r  =  1 ;
+ ON_A9s  =  1 ; ON_A9r  =  1 ; ON_B9s  =  1 ; ON_B9r  =  1 ;
+ ON_A10s =  1 ; ON_A10r =  1 ; ON_B10s =  1 ; ON_B10r =  1 ;
+ ON_A11s =  1 ; ON_A11r =  1 ; ON_B11s =  1 ; ON_B11r =  1 ;
+ ON_A12s =  ON_LSF; ON_A12r = ON_LSF; ON_B12s = ON_LSF; ON_B12r =  ON_LSF;
+ ON_A13s =  ON_LSF; ON_A13r = ON_LSF; ON_B13s = ON_LSF; ON_B13r =  ON_LSF;
+ ON_A14s =  ON_LSF; ON_A14r = ON_LSF; ON_B14s = ON_LSF; ON_B14r =  ON_LSF;
+ ON_A15s =  ON_LSF; ON_A15r = ON_LSF; ON_B15s = ON_LSF; ON_B15r =  ON_LSF;
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+! Define the seed number (for statistical error assignment in the new IT/D1)
+! and call the corresponding measured error table for nominal LHC magnets
+!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+!For the madx file
+! eoption,seed=1;
+! if (NRJ==450)  {readtable, file="wise/After_sector_3-4_repair/injection/injection_errors-emfqcs-1.tfs";}
+!!!!for the mask
+eoption,seed=1;
+! if (NRJ==450) {readtable, file="wise/After_sector_3-4_repair/injection/injection_errors-emfqcs-%seed_ran.tfs" ;}
+print, text="read errors from table";
+if (NRJ<4999.9999) {readtable, file="wise/After_sector_3-4_repair/injection/injection_errors-emfqcs-%seed_ran.tfs" ;}
+
+
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+!             apply field errors to MB and MQ magnets
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+print, text="apply field errors to MB and MQ magnets";
+ call,   file="optics2016/measured_errors/Efcomp_MB.madx"  ;
+ call,   file="optics2016/measured_errors/Efcomp_MQ.madx"  ;
+print, text="done applying field errors to MB and MQ magnets";
+
+!!!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+!!!           apply field errors to all other quads and dipoles
+!!!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+!!! Nominal LHC magnets
+print, text="apply field errors to all other quads and dips";
+ call,   file="optics2016/measured_errors/Efcomp_MBRB.madx"; ! for IP4
+ call,   file="optics2016/measured_errors/Efcomp_MBRS.madx"; ! for IP4
+ call,   file="optics2016/measured_errors/Efcomp_MBRC.madx"; ! this is the D2
+ call,   file="optics2016/measured_errors/Efcomp_MBX.madx" ; ! D in IP2 and 8
+ call,   file="optics2016/measured_errors/Efcomp_MBW.madx" ; ! this is for IP7 and IP3
+ call,   file="optics2016/measured_errors/Efcomp_MBXW.madx"; ! this is the D1 in IP1 and IP5
+!
+ call,   file="optics2016/measured_errors/Efcomp_MQM.madx" ;
+ call,   file="optics2016/measured_errors/Efcomp_MQMC.madx";
+ call,   file="optics2016/measured_errors/Efcomp_MQML.madx";
+ call,   file="optics2016/measured_errors/Efcomp_MQTL.madx";
+ call,   file="optics2016/measured_errors/Efcomp_MQW.madx" ;
+ call,   file="optics2016/measured_errors/Efcomp_MQY.madx" ;
+ call,   file="optics2016/measured_errors/Efcomp_MQX.madx" ; !!! Assume local corrections for coupling and optics errors in the IRs have been applied
+print, text="done applying field errors to all other quads and dips";
+
+
+!##############################################################################
+!########## GET NOMINAL CORRECTIONS FOR DIPOLES & TRIPLETS ETC ################
+!##############################################################################
+if(is_thin==1) {
+  value,NRJ;
+
+  select, flag=error, clear;
+  select, flag=error, pattern=MB\.,class=multipole;
+  esave,  file="temp/MB.errors";
+  system, "optics2016/errors/corr_MB_v5";
+  call, file="temp/MB_corr_setting.mad";
+  value, kcs.a12b1;
+  value, kcs.a78b1;
+  value, kcs.a12b2;
+  value, kcs.a78b2;
+
+  KSS.a12B1  := 0 ;
+  KSS.a23B1  := 0 ;
+  KSS.a34B1  := 0 ;
+  KSS.a45B1  := 0 ;
+  KSS.a56B1  := 0 ;
+  KSS.a67B1  := 0 ;
+  KSS.a78B1  := 0 ;
+  KSS.a81B1  := 0 ;
+
+  KCO.a12B1  := 0 ;
+  KCO.a23B1  := 0 ;
+  KCO.a34B1  := 0 ;
+  KCO.a45B1  := 0 ;
+  KCO.a56B1  := 0 ;
+  KCO.a67B1  := 0 ;
+  KCO.a78B1  := 0 ;
+  KCO.a81B1  := 0 ;
+};
+
+!MQS BEAM1
+CMRS.b1      :=   0.000000000000E+00 ;
+CMIS.b1      :=   0.000000000000E+00 ;
+CMRS.b1_sq   :=   0.000000000000E+00 ;
+CMIS.b1_sq   :=   0.000000000000E+00 ;
+ona2_b1      :=   0.000000000000E+00 ;
+ KQS.R1B1   :=  (   0.000000000000E+00) * ona2_b1 + ( 0.266278479040E-01) * CMRS.b1 + (-0.899733288016E-02) * CMIS.b1 + ( 0.000000000000E+00) * CMRS.b1_sq + ( 0.000000000000E+00) * CMIS.b1_sq  ;
+ KQS.L2B1   :=  (   0.000000000000E+00) * ona2_b1 + ( 0.266278479040E-01) * CMRS.b1 + (-0.899733288016E-02) * CMIS.b1 + ( 0.000000000000E+00) * CMRS.b1_sq + ( 0.000000000000E+00) * CMIS.b1_sq  ;
+KQS.A23B1   :=  (   0.000000000000E+00) * ona2_b1 + ( 0.142516736842E-01) * CMRS.b1 + ( 0.848602983914E-02) * CMIS.b1 + ( 0.302173156154E-01) * CMRS.b1_sq + ( 0.109843179604E-01) * CMIS.b1_sq  ;
+ KQS.R3B1   :=  (   0.000000000000E+00) * ona2_b1 + (-0.171205193764E-01) * CMRS.b1 + (-0.807870546221E-02) * CMIS.b1 + (-0.402300957758E-01) * CMRS.b1_sq + (-0.822964594698E-02) * CMIS.b1_sq  ;
+ KQS.L4B1   :=  (   0.000000000000E+00) * ona2_b1 + (-0.171205193764E-01) * CMRS.b1 + (-0.807870546221E-02) * CMIS.b1 + (-0.402300957758E-01) * CMRS.b1_sq + (-0.822964594698E-02) * CMIS.b1_sq  ;
+KQS.A45B1   :=  (   0.000000000000E+00) * ona2_b1 + ( 0.113812285983E-01) * CMRS.b1 + ( 0.955159460427E-02) * CMIS.b1 + ( 0.000000000000E+00) * CMRS.b1_sq + ( 0.000000000000E+00) * CMIS.b1_sq  ;
+ KQS.R5B1   :=  (   0.000000000000E+00) * ona2_b1 + ( 0.792323136002E-02) * CMRS.b1 + ( 0.100926247998E-01) * CMIS.b1 + ( 0.000000000000E+00) * CMRS.b1_sq + ( 0.000000000000E+00) * CMIS.b1_sq  ;
+ KQS.L6B1   :=  (   0.000000000000E+00) * ona2_b1 + ( 0.792323136002E-02) * CMRS.b1 + ( 0.100926247998E-01) * CMIS.b1 + ( 0.000000000000E+00) * CMRS.b1_sq + ( 0.000000000000E+00) * CMIS.b1_sq  ;
+KQS.A67B1   :=  (   0.000000000000E+00) * ona2_b1 + (-0.158692136780E-01) * CMRS.b1 + ( 0.106460324212E-01) * CMIS.b1 + (-0.709778694350E-01) * CMRS.b1_sq + ( 0.349381515069E-01) * CMIS.b1_sq  ;
+ KQS.R7B1   :=  (   0.000000000000E+00) * ona2_b1 + (-0.739140462540E-02) * CMRS.b1 + (-0.987710657697E-02) * CMIS.b1 + (-0.549901960504E-02) * CMRS.b1_sq + (-0.185504800255E-01) * CMIS.b1_sq  ;
+ KQS.L8B1   :=  (   0.000000000000E+00) * ona2_b1 + (-0.739140462540E-02) * CMRS.b1 + (-0.987710657697E-02) * CMIS.b1 + (-0.549901960504E-02) * CMRS.b1_sq + (-0.185504800255E-01) * CMIS.b1_sq  ;
+KQS.A81B1   :=  (   0.000000000000E+00) * ona2_b1 + ( 0.241324775639E-01) * CMRS.b1 + (-0.962582146500E-02) * CMIS.b1 + ( 0.000000000000E+00) * CMRS.b1_sq + ( 0.000000000000E+00) * CMIS.b1_sq  ;
+
+
+KQTF.A81B1  :=    0.000000000000E+00 + ( 0.386712697356E-02) * dQx.b1 + ( 0.714404818290E-03) * dQy.b1 + ( 0.000000000000E+00) * dQx.b1_sq + ( 0.000000000000E+00) * dQy.b1_sq  ;
+KQTF.A12B1  :=    0.000000000000E+00 + ( 0.386712697356E-02) * dQx.b1 + ( 0.714404818290E-03) * dQy.b1 + ( 0.000000000000E+00) * dQx.b1_sq + ( 0.000000000000E+00) * dQy.b1_sq  ;
+KQTF.A45B1  :=    0.000000000000E+00 + ( 0.386712697356E-02) * dQx.b1 + ( 0.714404818290E-03) * dQy.b1 + ( 0.000000000000E+00) * dQx.b1_sq + ( 0.000000000000E+00) * dQy.b1_sq  ;
+KQTF.A56B1  :=    0.000000000000E+00 + ( 0.386712697356E-02) * dQx.b1 + ( 0.714404818290E-03) * dQy.b1 + ( 0.000000000000E+00) * dQx.b1_sq + ( 0.000000000000E+00) * dQy.b1_sq  ;
+KQTD.A81B1  :=    0.000000000000E+00 + (-0.675151714265E-03) * dQx.b1 + (-0.361436046931E-02) * dQy.b1 + ( 0.000000000000E+00) * dQx.b1_sq + ( 0.000000000000E+00) * dQy.b1_sq  ;
+KQTD.A12B1  :=    0.000000000000E+00 + (-0.675151714265E-03) * dQx.b1 + (-0.361436046931E-02) * dQy.b1 + ( 0.000000000000E+00) * dQx.b1_sq + ( 0.000000000000E+00) * dQy.b1_sq  ;
+KQTD.A45B1  :=    0.000000000000E+00 + (-0.675151714265E-03) * dQx.b1 + (-0.361436046931E-02) * dQy.b1 + ( 0.000000000000E+00) * dQx.b1_sq + ( 0.000000000000E+00) * dQy.b1_sq  ;
+KQTD.A56B1  :=    0.000000000000E+00 + (-0.675151714265E-03) * dQx.b1 + (-0.361436046931E-02) * dQy.b1 + ( 0.000000000000E+00) * dQx.b1_sq + ( 0.000000000000E+00) * dQy.b1_sq  ;
+KQTF.A78B1  :=    0.761119694300E-03 + ( 0.386712697356E-02) * dQx.b1 + ( 0.714404818290E-03) * dQy.b1 + ( 0.724644665413E-02) * dQx.b1_sq + ( 0.132484451020E-02) * dQy.b1_sq  ;
+KQTF.A23B1  :=   -0.771991245000E-03 + ( 0.386712697356E-02) * dQx.b1 + ( 0.714404818290E-03) * dQy.b1 + ( 0.724644665413E-02) * dQx.b1_sq + ( 0.132484451020E-02) * dQy.b1_sq  ;
+KQTF.A34B1  :=   -0.771991245000E-03 + ( 0.386712697356E-02) * dQx.b1 + ( 0.714404818290E-03) * dQy.b1 + ( 0.724644665413E-02) * dQx.b1_sq + ( 0.132484451020E-02) * dQy.b1_sq  ;
+KQTF.A67B1  :=    0.761119694300E-03 + ( 0.386712697356E-02) * dQx.b1 + ( 0.714404818290E-03) * dQy.b1 + ( 0.724644665413E-02) * dQx.b1_sq + ( 0.132484451020E-02) * dQy.b1_sq  ;
+KQTD.A78B1  :=   -0.933337382400E-03 + (-0.675151714265E-03) * dQx.b1 + (-0.361436046931E-02) * dQy.b1 + (-0.134618584897E-02) * dQx.b1_sq + (-0.720260248926E-02) * dQy.b1_sq  ;
+KQTD.A23B1  :=    0.915761029700E-03 + (-0.675151714265E-03) * dQx.b1 + (-0.361436046931E-02) * dQy.b1 + (-0.134618584897E-02) * dQx.b1_sq + (-0.720260248926E-02) * dQy.b1_sq  ;
+KQTD.A34B1  :=    0.915761029700E-03 + (-0.675151714265E-03) * dQx.b1 + (-0.361436046931E-02) * dQy.b1 + (-0.134618584897E-02) * dQx.b1_sq + (-0.720260248926E-02) * dQy.b1_sq  ;
+KQTD.A67B1  :=   -0.933337382400E-03 + (-0.675151714265E-03) * dQx.b1 + (-0.361436046931E-02) * dQy.b1 + (-0.134618584897E-02) * dQx.b1_sq + (-0.720260248926E-02) * dQy.b1_sq  ;
+
+
+!###################################################################################################################################
+!############################## MISALIGN MCD #######################################################################################
+!###################################################################################################################################
+
+ SELECT, FLAG=ERROR, CLEAR;
+ SELECT,FLAG=ERROR,CLASS=MCO;
+ EOPTION,ADD=TRUE;
+ EALIGN,DX:=-0.00020,DY:=0.000;
+ SELECT, FLAG=ERROR, CLEAR;
+ SELECT,FLAG=ERROR,CLASS=MCD;
+ EOPTION,ADD=TRUE;
+ EALIGN,DX:=-0.00020,DY:=0.000;
+
+if(mylhcbeam==1){
+   match;
+     global, q1=qxinit, q2=qyinit;
+     vary,   name=dQx.b1, step=1.0E-7 ;
+     vary,   name=dQy.b1, step=1.0E-7 ;
+     lmdif,  calls=1000, tolerance=1.0E-21;
+   endmatch;
+   match;
+     global, dq1=dqxinit, dq2=dqyinit;
+     vary,   name=dQpx.b1, step=1.0E-7 ;
+     vary,   name=dQpy.b1, step=1.0E-7 ;
+     lmdif,  calls=1000, tolerance=1.0E-21;
+   endmatch;
+}
+
+!###########################################################################################################################################################################################
+!!###############!! APPLY THE MO
+!###########################################################################################################################################################################################
+if(mylhcbeam==1){
+ kof.a12b1 := MOFpowering;
+ kof.a23b1 := MOFpowering;
+ kof.a34b1 := MOFpowering;
+ kof.a45b1 := MOFpowering;
+ kof.a56b1 := MOFpowering;
+ kof.a67b1 := MOFpowering;
+ kof.a78b1 := MOFpowering;
+ kof.a81b1 := MOFpowering;
+
+ kod.a12b1 := MODpowering;
+ kod.a23b1 := MODpowering;
+ kod.a34b1 := MODpowering;
+ kod.a45b1 := MODpowering;
+ kod.a56b1 := MODpowering*0.775; !THIS MO IS HARDWARE LIMITED IN RUNii;
+ kod.a67b1 := MODpowering;
+ kod.a78b1 := MODpowering;
+ kod.a81b1 := MODpowering;
+}
+if(mylhcbeam>1){
+ kof.a12b2 := MOFpowering;
+ kof.a23b2 := MOFpowering;
+ kof.a34b2 := MOFpowering;
+ kof.a45b2 := MOFpowering;
+ kof.a56b2 := MOFpowering;
+ kof.a67b2 := MOFpowering;
+ kof.a78b2 := MOFpowering;
+ kof.a81b2 := MOFpowering;
+
+ kod.a12b2 := MODpowering;
+ kod.a23b2 := MODpowering;
+ kod.a34b2 := MODpowering;
+ kod.a45b2 := MODpowering;
+ kod.a56b2 := MODpowering;
+ kod.a67b2 := MODpowering;
+ kod.a78b2 := MODpowering;
+ kod.a81b2 := MODpowering;
+}
+
+ !MQS BEAM1
+CMRS.b1      :=   0.000000000000E+00 ;
+CMIS.b1      :=   0.000000000000E+00 ;
+CMRS.b1_sq   :=   0.000000000000E+00 ;
+CMIS.b1_sq   :=   0.000000000000E+00 ;
+ona2_b1      :=   0.000000000000E+00 ;
+ KQS.R1B1   :=  (   0.000000000000E+00) * ona2_b1 + ( 0.266278479040E-01) * CMRS.b1 + (-0.899733288016E-02) * CMIS.b1 + ( 0.000000000000E+00) * CMRS.b1_sq + ( 0.000000000000E+00) * CMIS.b1_sq  ;
+ KQS.L2B1   :=  (   0.000000000000E+00) * ona2_b1 + ( 0.266278479040E-01) * CMRS.b1 + (-0.899733288016E-02) * CMIS.b1 + ( 0.000000000000E+00) * CMRS.b1_sq + ( 0.000000000000E+00) * CMIS.b1_sq  ;
+KQS.A23B1   :=  (   0.000000000000E+00) * ona2_b1 + ( 0.142516736842E-01) * CMRS.b1 + ( 0.848602983914E-02) * CMIS.b1 + ( 0.302173156154E-01) * CMRS.b1_sq + ( 0.109843179604E-01) * CMIS.b1_sq  ;
+ KQS.R3B1   :=  (   0.000000000000E+00) * ona2_b1 + (-0.171205193764E-01) * CMRS.b1 + (-0.807870546221E-02) * CMIS.b1 + (-0.402300957758E-01) * CMRS.b1_sq + (-0.822964594698E-02) * CMIS.b1_sq  ;
+ KQS.L4B1   :=  (   0.000000000000E+00) * ona2_b1 + (-0.171205193764E-01) * CMRS.b1 + (-0.807870546221E-02) * CMIS.b1 + (-0.402300957758E-01) * CMRS.b1_sq + (-0.822964594698E-02) * CMIS.b1_sq  ;
+KQS.A45B1   :=  (   0.000000000000E+00) * ona2_b1 + ( 0.113812285983E-01) * CMRS.b1 + ( 0.955159460427E-02) * CMIS.b1 + ( 0.000000000000E+00) * CMRS.b1_sq + ( 0.000000000000E+00) * CMIS.b1_sq  ;
+ KQS.R5B1   :=  (   0.000000000000E+00) * ona2_b1 + ( 0.792323136002E-02) * CMRS.b1 + ( 0.100926247998E-01) * CMIS.b1 + ( 0.000000000000E+00) * CMRS.b1_sq + ( 0.000000000000E+00) * CMIS.b1_sq  ;
+ KQS.L6B1   :=  (   0.000000000000E+00) * ona2_b1 + ( 0.792323136002E-02) * CMRS.b1 + ( 0.100926247998E-01) * CMIS.b1 + ( 0.000000000000E+00) * CMRS.b1_sq + ( 0.000000000000E+00) * CMIS.b1_sq  ;
+KQS.A67B1   :=  (   0.000000000000E+00) * ona2_b1 + (-0.158692136780E-01) * CMRS.b1 + ( 0.106460324212E-01) * CMIS.b1 + (-0.709778694350E-01) * CMRS.b1_sq + ( 0.349381515069E-01) * CMIS.b1_sq  ;
+ KQS.R7B1   :=  (   0.000000000000E+00) * ona2_b1 + (-0.739140462540E-02) * CMRS.b1 + (-0.987710657697E-02) * CMIS.b1 + (-0.549901960504E-02) * CMRS.b1_sq + (-0.185504800255E-01) * CMIS.b1_sq  ;
+ KQS.L8B1   :=  (   0.000000000000E+00) * ona2_b1 + (-0.739140462540E-02) * CMRS.b1 + (-0.987710657697E-02) * CMIS.b1 + (-0.549901960504E-02) * CMRS.b1_sq + (-0.185504800255E-01) * CMIS.b1_sq  ;
+KQS.A81B1   :=  (   0.000000000000E+00) * ona2_b1 + ( 0.241324775639E-01) * CMRS.b1 + (-0.962582146500E-02) * CMIS.b1 + ( 0.000000000000E+00) * CMRS.b1_sq + ( 0.000000000000E+00) * CMIS.b1_sq  ;
+ !MQS BEAM2
+CMRS.b2      :=   0.000000000000E+00 ;
+CMIS.b2      :=   0.000000000000E+00 ;
+CMRS.b2_sq   :=   0.000000000000E+00 ;
+CMIS.b2_sq   :=   0.000000000000E+00 ;
+ona2_b2      :=   0.000000000000E+00 ;
+KQS.A12B2   :=  (   0.519213248459E-03) * ona2_b2 + ( 0.124458484817E-01) * CMRS.b2 + (-0.207596749726E-01) * CMIS.b2 + ( 0.000000000000E+00) * CMRS.b2_sq + ( 0.000000000000E+00) * CMIS.b2_sq  ;
+ KQS.R2B2   :=  (  -0.952785714091E-03) * ona2_b2 + ( 0.121456161967E-01) * CMRS.b2 + ( 0.397509215137E-02) * CMIS.b2 + ( 0.200920340631E-01) * CMRS.b2_sq + ( 0.712846355318E-02) * CMIS.b2_sq  ;
+ KQS.L3B2   :=  (  -0.952785714091E-03) * ona2_b2 + ( 0.121456161967E-01) * CMRS.b2 + ( 0.397509215137E-02) * CMIS.b2 + ( 0.200920340631E-01) * CMRS.b2_sq + ( 0.712846355318E-02) * CMIS.b2_sq  ;
+KQS.A34B2   :=  (  -0.217284900157E-03) * ona2_b2 + (-0.179258964749E-01) * CMRS.b2 + ( 0.371667140033E-02) * CMIS.b2 + (-0.405188393546E-01) * CMRS.b2_sq + ( 0.212647436428E-01) * CMIS.b2_sq  ;
+ KQS.R4B2   :=  (   0.186949354181E-03) * ona2_b2 + ( 0.128662566995E-01) * CMRS.b2 + ( 0.360504223421E-02) * CMIS.b2 + ( 0.000000000000E+00) * CMRS.b2_sq + ( 0.000000000000E+00) * CMIS.b2_sq  ;
+ KQS.L5B2   :=  (   0.186949354181E-03) * ona2_b2 + ( 0.128662566995E-01) * CMRS.b2 + ( 0.360504223421E-02) * CMIS.b2 + ( 0.000000000000E+00) * CMRS.b2_sq + ( 0.000000000000E+00) * CMIS.b2_sq  ;
+KQS.A56B2   :=  (   0.397814910113E-03) * ona2_b2 + ( 0.148124920807E-01) * CMRS.b2 + ( 0.138602241563E-02) * CMIS.b2 + ( 0.000000000000E+00) * CMRS.b2_sq + ( 0.000000000000E+00) * CMIS.b2_sq  ;
+ KQS.R6B2   :=  (  -0.138800614213E-03) * ona2_b2 + (-0.614988295119E-02) * CMRS.b2 + ( 0.179020058982E-01) * CMIS.b2 + (-0.327505145607E-01) * CMRS.b2_sq + ( 0.624417187334E-01) * CMIS.b2_sq  ;
+ KQS.L7B2   :=  (  -0.138800614213E-03) * ona2_b2 + (-0.614988295119E-02) * CMRS.b2 + ( 0.179020058982E-01) * CMIS.b2 + (-0.327505145607E-01) * CMRS.b2_sq + ( 0.624417187334E-01) * CMIS.b2_sq  ;
+KQS.A78B2   :=  (   0.907653646559E-03) * ona2_b2 + (-0.117937801373E-01) * CMRS.b2 + (-0.447830183312E-02) * CMIS.b2 + (-0.188089814503E-01) * CMRS.b2_sq + (-0.897287738080E-02) * CMIS.b2_sq  ;
+ KQS.R8B2   :=  (   0.173286377641E-03) * ona2_b2 + ( 0.144218507771E-01) * CMRS.b2 + (-0.211979276194E-01) * CMIS.b2 + ( 0.000000000000E+00) * CMRS.b2_sq + ( 0.000000000000E+00) * CMIS.b2_sq  ;
+ KQS.L1B2   :=  (   0.173286377641E-03) * ona2_b2 + ( 0.144218507771E-01) * CMRS.b2 + (-0.211979276194E-01) * CMIS.b2 + ( 0.000000000000E+00) * CMRS.b2_sq + ( 0.000000000000E+00) * CMIS.b2_sq  ;
+
+
+
+if(mylhcbeam==1){
+CMRS.b1      :=   0.0000 ;
+CMIS.b1      :=  -0.0000 ;
+}
+if(mylhcbeam>1){
+CMRS.b2      :=   0.0000 ;
+CMIS.b2      :=  -0.0000 ;
+}
+
+
+
+if(mylhcbeam==1){
+   match;
+     global, q1=qxinit, q2=qyinit;
+     vary,   name=dQx.b1, step=1.0E-7 ;
+     vary,   name=dQy.b1, step=1.0E-7 ;
+     lmdif,  calls=1000, tolerance=1.0E-21;
+   endmatch;
+   match;
+     global, dq1=dqxinit, dq2=dqyinit;
+     vary,   name=dQpx.b1, step=1.0E-7 ;
+     vary,   name=dQpy.b1, step=1.0E-7 ;
+     lmdif,  calls=1000, tolerance=1.0E-21;
+   endmatch;
+}
+
+
+
+twiss, file="twiss.tfs";
+
+if (NRJ<4999.9999) {VRF400:=%rf_vol; LAGRF400.B1=0.5; LAGRF400.B2=0.;};
+if (NRJ>5000.0000) {VRF400:=%rf_vol; LAGRF400.B1=0.5; LAGRF400.B2=0.;};
+! if (NRJ<4999.9999) {VRF400:=8.; LAGRF400.B1=0.5; LAGRF400.B2=0.;};
+! if (NRJ>5000.0000) {VRF400:=16.; LAGRF400.B1=0.5; LAGRF400.B2=0.;};
+
+twiss;
+!emit;
+ value, table(summ,q1);
+ value, table(summ,q2);
+ value, table(summ,dq1);
+ value, table(summ,dq2);
+
+sixtrack, cavall, radius=0.017;
+
+select, flag=twiss, clear;
+if (not_a_mask==1){
+  twiss,file="last_twiss.1";
+  System,"gzip -f last_twiss.1";
+} else {
+  twiss,file="last_twiss.1";
+  System,"gzip -f last_twiss.1";
+};
+stop;
diff --git a/templates/fort.3 b/templates/fort.3
index b409c57..051ad03 100644
--- a/templates/fort.3
+++ b/templates/fort.3
@@ -1,41 +1,55 @@
 GEOME-STRENG TITLE:%Runnam
 PRINTOUT OF INPUT PARAMETERS--------------------------------------------
 NEXT
+
 SETTING----------------------------------------------------------------
 INITIALSTATE text
 FINALSTATE text
 NEXT
-TRACKING PARAMETERS-----------------------------------------------------
-%turnss 0 %nss %ax0s %ax1s 0 %imc
-1 1 %idfor 1 %iclo6
-0 0 1 1 %writebins 50000 2
-NEXT
-INITIAL COORDINATES-----------------------------------------------------
-        2 0. 0. %ratio 0
-        0.
-        0.
-        0. 
-        0.
-        0.
-        %dp1
-        .000001 
-        0.
-        0. 
-        0.
-        0.
-        %dp2
-        %e_0
-        %e_0
-        %e_0
+
+SIMU
+  PARTICLES %nss
+  TURNS %turnss
+  REF_ENERGY %e_0
+  REF_PARTICLE proton 1 1 1
+
+  LATTICE thin 6D
+  6D_CLORB on
+
+  WRITE_TRACKS 1000
+NEXT
+
+RANDOM NUMBERS
+   SEED %sixtrack_seed
+NEXT
+
+DIST
+   FORMAT ID XN YN
+   READ input_dist.txt
+   FILL ZN FLOAT 0.0
+
+   FILL PXN FLOAT 0.0
+   FILL PYN FLOAT 0.0
+   FILL PZN FLOAT 1.666666
+
+   FILL ION_A INT 1
+   FILL ION_Z INT 1
+   FILL CHARGE INT 1
+   FILL MASS FLOAT %pmass
+   EMITTANCE %emit_norm %emit_norm
+   LONGEMIT 1.0 eVs
 NEXT
+
 SYNCHROTRON OSCILLATIONS--------'PLACE AFTER TRACKING PARAMETERS'-------
       35640 .000347 %rf_vol 0. %length %pmass %ition
       1.       1.
 NEXT
-BEAM
-%bunch_charge   %emit_norm_x  %emit_norm_y  %sig_z  %sig_e  1  %ibtype 1 0
 
-NEXT
+/ Beam block is actually already in the fort.3.mad file
+/ BEAM
+/ %bunch_charge   %emit_norm_x  %emit_norm_y  %sig_z  %sig_e  1  %ibtype 1 0
+/ NEXT
+
 ITERATION-ACCURACY------------------------------------------------------
        50 1D-12 1D-15
        10 1D-10 1D-10
@@ -47,6 +61,8 @@ FLUCTUATION
 100000 1 7 3
 NEXT
 
+%FILE: fort.3.mad
+
 %toggle_post/POSTPROCESSING----------------------------------------------------------
 %toggle_post/LHC Dynamic Aperture at Injection Version 5 (1998)
 %toggle_post/20 0 0 1 .08 .08 1 0 0 1. 1.
@@ -62,6 +78,7 @@ NEXT
 %toggle_coll/LINE
 %toggle_coll/ELEMENT 0 %ilin 1 %EI %EI
 %toggle_coll/NEXT
+
 %toggle_coll/COLLIMATION
 %toggle_coll/.TRUE. 
 %toggle_coll/1 6500000
@@ -81,12 +98,6 @@ NEXT
 %toggle_coll/"CollDB.data" 1
 %toggle_coll/.FALSE. .FALSE. WAbsVertLowbcoll 101 1 1.
 %toggle_coll/NEXT
-%toggle_coll/LIMI
-%toggle_coll// DEBUG
-%toggle_coll// PREC 0.001
-%toggle_coll// PRIN ape_dump.dat
-%toggle_coll/LOAD fort3.limi
-%toggle_coll/NEXT
 
 ENDE
 
diff --git a/templates/hl10.mask b/templates/hl10.mask
deleted file mode 100644
index a684cbb..0000000
--- a/templates/hl10.mask
+++ /dev/null
@@ -1,1011 +0,0 @@
-! S. Fartoukh March 2009
-! Ref sample job for tracking with multipole field imperfection in nominal 
-! LHC magnets using V6.503 in collision
-! M. Giovannozzi November 2010
-
-!---------------------------------------------------------------;
-!			FLAGS
-!---------------------------------------------------------------;
-
-mylhcbeam = %lhc_beam ;! LHC beam 1 (clockwise)
-
-correct_for_MCBX=1;
-nominal_MCBXF_signs=0;
-ON_MBH=0;
-thick_lenses=0;
-create_errorfiles_and_quit=0;
-
-!only one of these three:
-chooseMQXF=0;
-chooseMQXFAB=0;
-chooseMQXFbodyends=1;
-
-
-! strengths for MCBX at collision
-A_col_a3 =  20.1200;
-A_col_b3 = -16.6500;
-B_col_a3 = -10.3300;
-B_col_b3 =  17.3700;
-
-A_col_a5 = -3.0400;
-A_col_b5 = -0.3500;
-B_col_a5 = -3.6000;
-B_col_b5 =  2.4900;
-
-A_col_a7 = -3.9800;
-A_col_b7 =  0.9800;
-B_col_a7 = -3.2600;
-B_col_b7 =  0.6200;
-
-A_col_a9 = -0.6200;
-A_col_b9 =  0.0700;
-B_col_a9 = -0.5800;
-B_col_b9 = -0.7500;
-
-A_col_a11 = 0.0200;
-A_col_b11 = 4.3000;
-B_col_a11 = 0.1200;
-B_col_b11 = 3.6000;
-
-! 3th order strengths for MCBX at injection
-A_inj_a3 = A_col_a3;
-A_inj_b3 = A_col_b3;
-B_inj_a3 = B_col_a3;
-B_inj_b3 = B_col_b3;
-A_inj_a5 = A_col_a5;
-A_inj_b5 = A_col_b5;
-B_inj_a5 = B_col_a5;
-B_inj_b5 = B_col_b5;
-A_inj_a7 = A_col_a7;
-A_inj_b7 = A_col_b7;
-B_inj_a7 = B_col_a7;
-B_inj_b7 = B_col_b7;
-A_inj_a9 = A_col_a9;
-A_inj_b9 = A_col_b9;
-B_inj_a9 = B_col_a9;
-B_inj_b9 = B_col_b9;
-A_inj_a11 = A_col_a11;
-A_inj_b11 = A_col_b11;
-B_inj_a11 = B_col_a11;
-B_inj_b11 = B_col_b11;
-
-
-!---|----1----|----2----|----3----|----4----|----5----|----6----|----7----|---8;
-Option, warn, info;
-System, "rm -rf temp";
-System, "mkdir temp";
-System, "ln -fns /afs/cern.ch/eng/lhc/optics/V6.503 db5";
-System, "ln -fns /afs/cern.ch/eng/lhc/optics/HLLHCV1.0 slhc";
-system, "ln -fns /afs/cern.ch/eng/lhc/optics/errors/0705 wise";
-system, "ln -fns /afs/cern.ch/eng/lhc/optics/V6.503/WISE/After_sector_3-4_repair fidel";
-system, "ln -fns /afs/cern.ch/user/f/fvanderv/public/MCBX fre";
-option, -echo, -info;
-
-
-
-!---------------------------------------------------------------;
-!           CALL IN SELECTION, MACRO's, SEQUENCE, BEAM, NRJ and OPTICS FILES;
-!---------------------------------------------------------------;
-call,file="slhc/toolkit/macro.madx";        !Call some macro facilities
-
-if (ON_MBH==1 || thick_lenses==1){
-		!!********************************************************************
-	!!Manual Thin Sequence building
-	!!********************************************************************
-	Option, -echo, -warn, -info;
-	if (mylhcbeam==4){
-	   call, file="db5/beam_four.seq"; bvaux=1;
-	} else {
-	   call, file="db5/V6.5.seq"; bvaux=-1;
-	};
-	Option, -echo, warn, -info;
-
-	!!!Install MBH
-	call, file="slhc/hllhc_sequence.madx";
-
-	if (ON_MBH==1){
-		!!Install MBH
-		mbh_rm_mcs=0; mbh_rm_mcdo=0; !=1 to remove spools
-		mbh_ir2q8=0; mbh_ir2q10=1;
-		mbh_ir7q8=1; mbh_ir7q10=1;
-		mbh_ir1q8=0; mbh_ir1q9=1; mbh_ir1q10=0;
-		mbh_ir5q8=0; mbh_ir5q9=1; mbh_ir5q10=0;
-		call, file="slhc/toolkit/install_mbh.madx";
-	};
-
-	if (thick_lenses==0){
-		! Slice nominal sequence
-		exec, myslice;
-	};
-} else {
-	!********************************************************************
-	!Saved Thin Sequence
-	!missing variables like lhclength and so on
-	!********************************************************************
-	Option, -echo, -warn, -info;
-	if (mylhcbeam==4){
-	   call, file="slhc/hllhc_thinb4.seq"; bvaux=1;
-	} else {
-	   call, file="slhc/hllhc_thin.seq"; bvaux=-1;
-	};
-	Option, -echo, warn, -info;
-};
-
-
-call,file="slhc/errors2/install_mqxf_fringenl.madx"; 	! adding fringe place holder
-call,file="slhc/errors2/install_MCBXAB_errors.madx";     ! adding multipole correctors in IR1/5 (for errors)
-
-!Cycling w.r.t. to IP3 (mandatory to find closed orbit in collision in the presence of errors)
-if (mylhcbeam<3){
-	seqedit, sequence=lhcb1; flatten; cycle,start=IP3; flatten; endedit;
-};
-seqedit, sequence=lhcb2; flatten; cycle,start=IP3; flatten; endedit;
-
-!********************************************
-NRJ = %e_0 * 1e-3; ! [MeV] --> [GeV]
-I_MO = %i_mo; !-570
-
-!General switch to select collision (0/1)
-ON_COLLISION := 0;
-!General switch to install bb lens (0/1)
-ON_BB_SWITCH := 0;
-
-b_t_dist := %b_t_dist; !25 bunch separation [ns]
-emit_norm_x := %emit_norm_x * 1e-6; ! [um] --> [m]
-emit_norm_y := %emit_norm_y * 1e-6; ! [um] --> [m]
-Nb_0 := %bunch_charge;
-sigt := %sig_z; ! bunch length [m] in collision
-sige := %sig_e;
-gamma_rel := NRJ/pmass;
-epsx := emit_norm_x /gamma_rel;
-epsy := emit_norm_y /gamma_rel;
-
-!********************************************
-
-qx00 = %int_tune_x;  !integer h tune
-qy00 = %int_tune_y;  !integer v tune
-tsplit = qx00 - qy00;   !new tune split
-qx0 = %tune_x;
-qy0 = %tune_y;
-qprimex = %q_prime_x;
-qprimey = %q_prime_y;
-
-if (NRJ<4999.9999){
-	!Inj. optics in thin lens
-	call,file="slhc/opt_inj_thin.madx";  ! beta* [m]=6/10/6/10 in IR1/2/5/8
-};
-if (NRJ>5000.0000){
-	!Coll. optics in thin lens for IR1&5 (various possible beta*_x/y in IR1/5)
-	call,file="slhc/opt_round_thin.madx";  ! beta*=15cm in IR1/5, beta*=10 m in IR2, 3m in IR8.
-	!call,file="slhc/opt_sround_thin.madx"; ! beta*=10cm in IR1/5, beta*=10 m in IR2, 3m in IR8.
-	!call,file="slhc/opt_flat_thin.madx";   ! beta*=7.5/30cm in IR1, beta*=30/7.5cm in IR5, beta*=10 m in IR2, 3m in IR8.
-	!call,file="slhc/opt_sflat_thin.madx";  ! beta*=5.0/20cm in IR1, beta*=20/5.0cm in IR5, beta*=10 m in IR2, 3m in IR8.
-	!call,file="slhc/opt_flathv_thin.madx"; ! beta*=30/7.5cm in IR1, beta*=7.5/30cm in IR5, beta*=10 m in IR2, 3m in IR8.
-	!call,file="slhc/opt_sflathv_thin.madx";! beta*=20/5.0cm in IR1, beta*=5.0/20cm in IR5, beta*=10 m in IR2, 3m in IR8.
-	!call,file="slhc/opt_180_180_180_180_thin.madx";
-	!call,file="slhc/opt_220_220_220_220_thin.madx";
-	!call,file="slhc/opt_260_260_260_260_thin.madx";
-	!call,file="slhc/opt_300_300_300_300_thin.madx";
-	!call,file="slhc/opt_340_340_340_340_thin.madx";
-	!call,file="slhc/opt_380_380_380_380_thin.madx";
-	!call,file="slhc/opt_420_420_420_420_thin.madx";
-	!call,file="slhc/opt_presqueeze_thin.madx";
-	!call,file="slhc/opt_presqueeze1000_thin.madx";
-	ON_QPP := 0;  ! Correction of residual Q'' by MO's
-	ON_DISP := 1;  ! Correction of spurious dispersion
-};
-
-if (mylhcbeam<3){
-	Beam, particle=proton, sequence=lhcb1, energy=NRJ, sigt=sigt, bv=1, NPART=Nb_0, sige=sige, ex=epsx, ey=epsy;
-};
-Beam, particle=proton, sequence=lhcb2, energy=NRJ, sigt=sigt, bv=bvaux, NPART=Nb_0, sige=sige, ex=epsx, ey=epsy;
-
-
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++;
-!   Rematch IP1 IP5 phase
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++;
-! nominal round: mux_ip15 = 31.21038468 ; muy_ip15 = 30.37288037 ;
-! nominal injection: mux_ip15 = 31.19538482 ; muy_ip15 = 30.36788031 ;
-
-!mux_ip15 = 31.21038468;
-!muy_ip15 = 31.19538482;
-
-!call,file="slhc/toolkit/make_ip15phase.madx";
-!call,file=slhc/toolkit/make_ip15phase_tm.madx; ! use phase-trombone instead
-!call,file=slhc/toolkit/delete_ip15phase_tm.madx; ! remove phase trombone
-
-! test new optics
-if (mylhcbeam==1){
-	exec, check_ip(b1);
-} else {
-	exec, check_ip(b2);
-}
-mux_ip15_ref=table(twiss,IP1,mux)-table(twiss,IP5,mux);
-muy_ip15_ref=table(twiss,IP1,muy)-table(twiss,IP5,muy);
-value, mux_ip15_ref, muy_ip15_ref;
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++;
-!Set crossing angle and separations knobs
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++;
-phi_IR1 = %phi_IR1; phi_IR5 = %phi_IR5; ! For round and flat optics
-!phi_IR1 =    0.000; phi_IR5 =   90.000; ! For flathv optics
-xing15=295e-6; sep15=2e-3;
-xing8=-250e-6; sep8=-2e-3; on_lhcb=+1;
-xing2= 170e-6; sep2=2e-3; on_alice=1;
-
-on_x1=xing15/(pxIP1b1+pyIP1b1);
-on_sep1=sep15/(xIP1b1+yIP1b1);
-on_x5=xing15/(pxIP5b1+pyIP5b1);
-on_sep5=sep15/(xIP5b1+yIP5b1);
-on_x2=xing2/(pyIP2b1-69.9994e-6); !removing alice=1 from pyIP2b1
-on_sep2=sep2/xIP2b1;
-on_x8=xing8/(pxIP8b1+134.9879e-6); !!removing on_lhcb from pxIP8b1
-on_sep8=sep8/yIP8b1;
-
-value, on_x1, on_sep1, on_x5, on_sep5, on_x2, on_sep2, on_x8, on_sep8;
-
-on_disp=1;
-
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++;
-!                        TEST the optics
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++;
-!/*
-	! Saving crossing scheme knobs
-	on_x1aux=on_x1; on_sep1aux=on_sep1;
-	on_x2aux=on_x2; on_sep2aux=on_sep2;
-	on_x5aux=on_x5; on_sep5aux=on_sep5;
-	on_x8aux=on_x8; on_sep8aux=on_sep8; on_a8aux=on_a8;
-	on_aliceaux=on_alice; on_lhcbaux=on_lhcb;
-	on_dispaux=on_disp;
-
-	on_x1=0; on_x5=0; on_sep1=0; on_sep5=0;
-	on_x2=0; on_x8=0; on_sep2=0; on_sep8=0; on_a8=0;
-	on_lhcb=0; on_alice=0;
-	on_disp=0;
-	system,"rm -f twiss.b1";system,"rm -f twiss.b2";
-	if (mylhcbeam<3){
-		Use, sequence=lhcb1;
-		select, flag=twiss, clear;
-		select, flag=twiss, pattern="IP1",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
-		select, flag=twiss, pattern="IP2",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
-		select, flag=twiss, pattern="IP5",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
-		select, flag=twiss, pattern="IP8",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
-		twiss,sequence=lhcb1,file=twiss.b1;system,"cat twiss.b1";
-	};
-
-	Use, sequence=lhcb2;
-	select, flag=twiss, clear;
-	select, flag=twiss, pattern="IP1",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
-	select, flag=twiss, pattern="IP2",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
-	select, flag=twiss, pattern="IP5",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
-	select, flag=twiss, pattern="IP8",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
-	twiss, sequence=lhcb2,file=twiss.b2;system,"cat twiss.b2";
-
-	system,"rm -f twiss.b1";system,"rm -f twiss.b2";
-	!Switch on Xscheme
-	on_x1=on_x1aux;on_sep1=on_sep1aux;
-	on_x2=on_x2aux;on_sep2=on_sep2aux;
-	on_x5=on_x5aux;on_sep5=on_sep5aux;
-	on_x8=on_x8aux;on_sep8=on_sep8aux; on_a8=an_a8aux;
-	on_alice=on_aliceaux; on_lhcb=on_lhcbaux;
-	on_disp=on_dispaux;
-
-	if (mylhcbeam<3){
-		Use, sequence=lhcb1;
-		select, flag=twiss, clear;
-		select, flag=twiss, pattern="IP1",   column=name,s,x,y,px,py;
-		select, flag=twiss, pattern="IP2",   column=name,s,x,y,px,py;
-		select, flag=twiss, pattern="IP5",   column=name,s,x,y,px,py;
-		select, flag=twiss, pattern="IP8",   column=name,s,x,y,px,py;
-		twiss, sequence=lhcb1, file=twiss.b1;system,"cat twiss.b1";
-	};
-
-	Use, sequence=lhcb2;
-	select, flag=twiss, clear;
-	select, flag=twiss, pattern="IP1",   column=name,s,x,y,px,py;
-	select, flag=twiss, pattern="IP2",   column=name,s,x,y,px,py;
-	select, flag=twiss, pattern="IP5",   column=name,s,x,y,px,py;
-	select, flag=twiss, pattern="IP8",   column=name,s,x,y,px,py;
-	twiss, sequence=lhcb2, file=twiss.b2;system,"cat twiss.b2";
-!*/
-system,"rm -f twiss.b1";system,"rm -f twiss.b2";
-
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++;
-!     Xscheme definition   & Beam-beam encounters installation
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++;
-call,file="slhc/beambeam/macro_bb.madx";                  ! macros for beam-beam
-
-b_t_dist = %b_t_dist;        !bunch separation [ns]
-
-fraction_crab := %fraction_crab;  !Switch on/off the crab of the strong beam (1 means 100%, -2 means -200%,..)
-
-n_insideD1 = 5;    !default value for the number of additionnal parasitic encounters inside D1
-
-nho_IR1= 5;        ! number of slices for head-on in IR1 (between 0 and 201)   
-nho_IR2= 5;        ! number of slices for head-on in IR2 (between 0 and 201)
-nho_IR5= 5;        ! number of slices for head-on in IR5 (between 0 and 201)
-nho_IR8= 5;        ! number of slices for head-on in IR8 (between 0 and 201)
-
-exec DEFINE_BB_PARAM;  !Define main beam-beam parameters
-
-if (ON_COLLISION==1){
-	!Switch on Xscheme in precollision
-	on_x1aux=on_x1; on_sep1aux=1;
-	on_x2aux=on_x2; on_sep2aux=1;
-	on_x5aux=on_x5; on_sep5aux=1;
-	on_x8aux=on_x8; on_sep8aux=1; on_a8aux=on_a8;
-	on_aliceaux=on_alice; on_lhcbaux=on_lhcb;
-	on_disp:=0;
-	halo1=0; halo2=5; halo5=0; halo8=0;  !halo collision at 5 sigma's in Alice
-	!Redefine the on_sep's accordingly
-	exec CALCULATE_XSCHEME($halo1,$halo2,$halo5,$halo8);
-	on_sep1aux=on_sep1;
-	on_sep2aux=on_sep2;
-	on_sep5aux=on_sep5;
-	on_sep8aux=on_sep8;
-};
-
-! Saving crossing scheme knobs
-on_x1aux=on_x1; on_sep1aux=on_sep1;
-on_x2aux=on_x2; on_sep2aux=on_sep2;
-on_x5aux=on_x5; on_sep5aux=on_sep5;
-on_x8aux=on_x8; on_sep8aux=on_sep8; on_a8aux=on_a8;
-on_aliceaux=on_alice; on_lhcbaux=on_lhcb;
-!on_dispaux=on_disp;
-on_disp=0; !see old mask don't know why
-
-!Record the nominal IP position and crossing angle
-if(mylhcbeam==1) {use, sequence=lhcb1;};
-if(mylhcbeam>1) {use, sequence=lhcb2;};
-twiss;
-xnom1=table(twiss,IP1,x); pxnom1=table(twiss,IP1,px); ynom1=table(twiss,IP1,y); pynom1=table(twiss,IP1,py);
-xnom2=table(twiss,IP2,x); pxnom2=table(twiss,IP2,px); ynom2=table(twiss,IP2,y); pynom2=table(twiss,IP2,py);
-xnom5=table(twiss,IP5,x); pxnom5=table(twiss,IP5,px); ynom5=table(twiss,IP5,y); pynom5=table(twiss,IP5,py);
-xnom8=table(twiss,IP8,x); pxnom8=table(twiss,IP8,px); ynom8=table(twiss,IP8,y); pynom8=table(twiss,IP8,py);
-value, xnom1, xnom2, xnom5, xnom8;
-value, ynom1, ynom2, ynom5, ynom8;
-value, pxnom1, pxnom2, pxnom5, pxnom8;
-value, pynom1, pynom2, pynom5, pynom8;
-
-if (ON_BB_SWITCH==1){
-	!Install b-b marker
-	exec INSTALL_BB_MARK(b1);
-	exec INSTALL_BB_MARK(b2);
-
-	!Define bb lenses for both beams in all IR's and calculate # of encounters before D1
-	exec CALCULATE_BB_LENS;
-
-	!Install bb lenses
-	on_ho1 =1; on_lr1l = 1; on_lr1r = 1; npara_1 = npara0_1 + n_insideD1;
-	on_ho5 =1; on_lr5l = 1; on_lr5r = 1; npara_5 = npara0_5 + n_insideD1;
-	on_ho2 =1; on_lr2l = 1; on_lr2r = 1; npara_2 = npara0_2 + n_insideD1;
-	on_ho8 =1; on_lr8l = 1; on_lr8r = 1; npara_8 = npara0_8 + n_insideD1;
-	if(mylhcbeam==1) {exec INSTALL_BB_LENS(b1);};
-	if(mylhcbeam>1) {exec INSTALL_BB_LENS(b2);};
-
-	ON_BB_CHARGE := 0; !Switch off the charge the bb lenses
-
-	/*
-		!Plot b-b separation
-		exec PLOT_BB_SEP(1,$npara_1);
-		exec PLOT_BB_SEP(2,$npara_2);
-		exec PLOT_BB_SEP(5,$npara_5);
-		exec PLOT_BB_SEP(8,$npara_8);
-	*/
-
-	!/*
-		!if (mylhcbeam<3){
-		!	seqedit,sequence=lhcb1;flatten;cycle,start=IP1;endedit;
-		!};
-		!seqedit,sequence=lhcb2;flatten;cycle,start=IP1;endedit;
-
-		if (mylhcbeam<3){
-			use,sequence=lhcb1;
-			select,flag=twiss,clear;
-			select,flag=twiss,class=marker,pattern=PAR.*L1,range=mbxa.4l1..4/IP1.L1,column=s,name,betx,bety,alfx,alfy,mux,muy,x,y,px,py;
-			select,flag=twiss,class=marker,pattern=PAR.*L5,range=mbxa.4l5..4/IP5,column=s,name,betx,bety,alfx,alfy,mux,muy,x,y,px,py;
-			select,flag=twiss,class=marker,pattern=PAR.*R1,range=IP1/mbxa.4r1..1,column=s,name,betx,bety,alfx,alfy,mux,muy,x,y,px,py;
-			select,flag=twiss,class=marker,pattern=PAR.*R5,range=IP5/mbxa.4r5..1,column=s,name,betx,bety,alfx,alfy,mux,muy,x,y,px,py;
-			select,flag=twiss,class=marker,pattern=IP1,column=s,name,betx,bety,alfx,alfy,mux,muy,x,y,px,py;
-			select,flag=twiss,class=marker,pattern=IP5,column=s,name,betx,bety,alfx,alfy,mux,muy,x,y,px,py;
-			twiss,file=twiss_bb.b1;system,"cat twiss_bb.b1";
-		};
-
-		use,sequence=lhcb2;
-		select,flag=twiss,clear;
-		select,flag=twiss,class=marker,pattern=PAR.*L1,range=mbxa.4l1..4/IP1.L1,column=s,name,betx,bety,alfx,alfy,mux,muy,x,y,px,py;
-		select,flag=twiss,class=marker,pattern=PAR.*L5,range=mbxa.4l5..4/IP5,column=s,name,betx,bety,alfx,alfy,mux,muy,x,y,px,py;
-		select,flag=twiss,class=marker,pattern=PAR.*R1,range=IP1/mbxa.4r1..1,column=s,name,betx,bety,alfx,alfy,mux,muy,x,y,px,py;
-		select,flag=twiss,class=marker,pattern=PAR.*R5,range=IP5/mbxa.4r5..1,column=s,name,betx,bety,alfx,alfy,mux,muy,x,y,px,py;
-		select,flag=twiss,class=marker,pattern=IP1,column=s,name,betx,bety,alfx,alfy,mux,muy,x,y,px,py;
-		select,flag=twiss,class=marker,pattern=IP5,column=s,name,betx,bety,alfx,alfy,mux,muy,x,y,px,py;
-		twiss,file=twiss_bb.b2;system,"cat twiss_bb.b2";
-
-		if(mylhcbeam==1) {use,sequence=lhcb1;};
-		if(mylhcbeam>1) {use,sequence=lhcb2;};
-
-		select,flag=twiss,clear;
-		select,flag=twiss,pattern=HO,class=beambeam,column=s,name,betx,bety,alfx,alfy,mux,muy,x,y,px,py;
-		twiss,file=twiss_bb;system,"cat twiss_bb";
-
-		!if (mylhcbeam<3){
-		!	seqedit,sequence=lhcb1;flatten;cycle,start=IP3;endedit;
-		!};
-		!seqedit,sequence=lhcb2;flatten;cycle,start=IP3;endedit;
-	!*/
-
-	!Remove bb markers
-	exec REMOVE_BB_MARKER;
-
-	/*
-		!Make and plot footprint (at 6 sigmas)
-		ON_BB_CHARGE := 1;
-
-		!Switch on Xscheme
-		on_x1=on_x1aux;on_sep1=on_sep1aux;
-		on_x2=on_x2aux;on_sep2=on_sep2aux;
-		on_x5=on_x5aux;on_sep5=on_sep5aux;
-		on_x8=on_x8aux;on_sep8=on_sep8aux; on_a8=an_a8aux;
-		on_alice=on_aliceaux; on_lhcb=on_lhcbaux;
-		on_disp=on_dispaux;
-
-		nsigmax=6;
-
-		if(qx0-qx00<0.3){
-			if(mylhcbeam==1) {exec MAKEFOOTPRINT(b1);exec PLOTFOOTPRINT(b1,0.2795,0.2805,0.3095,0.3105);};
-			if(mylhcbeam>1) {exec MAKEFOOTPRINT(b2);exec PLOTFOOTPRINT(b2,0.2795,0.2805,0.3095,0.3105);};
-		};
-		if(qx0-qx00>0.3){
-			if(mylhcbeam==1) {exec MAKEFOOTPRINT(b1);exec PLOTFOOTPRINT(b1,0.300,0.315,0.310,0.325);};
-			if(mylhcbeam>1) {exec MAKEFOOTPRINT(b2);exec PLOTFOOTPRINT(b2,0.300,0.315,0.310,0.325);};
-		};
-		ON_BB_CHARGE := 0;
-		on_x1=0;on_sep1=0;
-		on_x2=0;on_sep2=0;
-		on_x5=0;on_sep5=0;
-		on_x8=0;on_sep8=0; on_a8=0;
-		on_alice=0; on_lhcb=0;
-		on_disp=0;
-	*/
-
-	!Remove bb lens for both beams
-	!exec REMOVE_BB_LENS;
-};
-
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++;
-!               Prepare nominal twiss tables
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++;
-if (mylhcbeam==1){
-	use,sequence=lhcb1;
-} else {
-	use,sequence=lhcb2;
-};
-
-/*
-	twiss;
-	x.ip1=table(twiss,IP1,x);px.ip1=table(twiss,IP1,px);
-	y.ip1=table(twiss,IP1,y);px.ip1=table(twiss,IP1,py); //irrelevant when recycling the sequences
-*/
-
-twiss, table=nominal;   ! used by orbit correction
-beta.ip1=table(twiss,IP1,betx);value,beta.ip1;
-
-! print nominal optics parameter at the MB, MQS and MSS for
-! b2, b3, b4, b5, a2 and a3 correction
-select, flag=twiss, clear;
-select, flag=twiss, pattern=MB\.   , class=multipole, column=name, k0L, k1L, betx, bety, dx, mux, muy;
-select, flag=twiss, pattern=MBH\.  , class=multipole, column=name, k0L, k1L, betx, bety, dx, mux, muy;
-select, flag=twiss, pattern=MQT\.14, class=multipole, column=name, k0L, k1L, betx, bety, dx, mux, muy;
-select, flag=twiss, pattern=MQT\.15, class=multipole, column=name, k0L, k1L, betx, bety, dx, mux, muy;
-select, flag=twiss, pattern=MQT\.16, class=multipole, column=name, k0L, k1L, betx, bety, dx, mux, muy;
-select, flag=twiss, pattern=MQT\.17, class=multipole, column=name, k0L, k1L, betx, bety, dx, mux, muy;
-select, flag=twiss, pattern=MQT\.18, class=multipole, column=name, k0L, k1L, betx, bety, dx, mux, muy;
-select, flag=twiss, pattern=MQT\.19, class=multipole, column=name, k0L, k1L, betx, bety, dx, mux, muy;
-select, flag=twiss, pattern=MQT\.20, class=multipole, column=name, k0L, k1L, betx, bety, dx, mux, muy;
-select, flag=twiss, pattern=MQT\.21, class=multipole, column=name, k0L, k1L, betx, bety, dx, mux, muy;
-select, flag=twiss, class=MQS                       , column=name, k0L, k1L, betx, bety, dx, mux, muy;
-select, flag=twiss, class=MSS                       , column=name, k0L, k1L, betx, bety, dx, mux, muy;
-select, flag=twiss, class=MCO                       , column=name, k0L, k1L, betx, bety, dx, mux, muy;
-select, flag=twiss, class=MCD                       , column=name, k0L, k1L, betx, bety, dx, mux, muy;
-select, flag=twiss, class=MCS                       , column=name, k0L, k1L, betx, bety, dx, mux, muy;
-twiss,  file='temp/optics0_MB.mad';
-
-! print nominal optics parameter at the D1, MQX and triplet corrector
-! for triplet correction
-select, flag=twiss, clear;
-select, flag=twiss, pattern=MQX  , class=multipole, column=name, betx, bety, x, y;
-select, flag=twiss, pattern=MBX  , class=multipole, column=name, betx, bety, x, y;
-select, flag=twiss, pattern=MBRC , class=multipole, column=name, betx, bety, x, y;
-select, flag=twiss, pattern=MBRD , class=multipole, column=name, betx, bety, x, y;
-select, flag=twiss, pattern=MQSX                  , column=name, betx, bety, x, y;
-select, flag=twiss, pattern=MCQSX                 , column=name, betx, bety, x, y;
-select, flag=twiss, pattern=MCSX                  , column=name, betx, bety, x, y;
-select, flag=twiss, pattern=MCTX                  , column=name, betx, bety, x, y;
-select, flag=twiss, pattern=MCOSX                 , column=name, betx, bety, x, y;
-select, flag=twiss, pattern=MCOX                  , column=name, betx, bety, x, y;
-select, flag=twiss, pattern=MCSSX                 , column=name, betx, bety, x, y;
-select, flag=twiss, pattern=MCDX                  , column=name, betx, bety, x, y;
-select, flag=twiss, pattern=MCDSX                 , column=name, betx, bety, x, y;
-select, flag=twiss, pattern=MCTSX                 , column=name, betx, bety, x, y;
-if (correct_for_MCBX==1){
-	select, flag=twiss, pattern=MCBXF, class=multipole, column=name, betx, bety, x, y;
-};
-twiss, file='temp/optics0_inser.mad';
-
-! disable crossing bumps
-on_x1=0; on_sep1=0;
-on_x2=0; on_sep2=0;
-on_x5=0; on_sep5=0;
-on_x8=0; on_sep8=0; on_a8=0;
-on_alice=0; on_lhcb=0;
-on_disp=0;
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-!                   Align separation magnets
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-/*
-    if(mylhcbeam==1){
-		call, file = "db5/measured_errors/align_D1_slices.b1.madx";
-		call, file = "db5/measured_errors/align_D2_slices.b1.madx";
-		call, file = "db5/measured_errors/align_D3_slices.b1.madx";
-		call, file = "db5/measured_errors/align_D4_slices.b1.madx";
-	};
-	if(mylhcbeam>1){
-		call, file = "db5/measured_errors/align_D1_slices.b2.madx";
-		call, file = "db5/measured_errors/align_D2_slices.b2.madx";
-		call, file = "db5/measured_errors/align_D3_slices.b2.madx";
-		call, file = "db5/measured_errors/align_D4_slices.b2.madx";
-	};
-*/
-call,file="slhc/toolkit/align_sepdip.madx";
-!exec,align_mbxw;  ! V6.503  D1
-!exec,align_mbrc15;! V6.503 D2 in IR15
-exec,align_mbx15; ! HL-LHC D1
-exec,align_mbrd15;! HL-LHC D2 in IR15
-
-exec,align_mbx28; ! V6.503 D1 in IR28
-exec,align_mbrc28;! V6.503 D2 in IR28
-exec,align_mbrs;  ! V6.503 D3 in IR4
-exec,align_mbrb;  ! V6.503 D4 in IR4
-
-if (ON_MBH==1) {
-	call,file="slhc/toolkit/align_mbh.madx"; ! align 11T dipoles
-};
-
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-!         Call error subroutines (nominal machine and new IT/D1)
-!                         and error tables
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-! Error routine and measured error table for nominal LHC
-
-call,file="db5/measured_errors/Msubroutines_new.madx";
-call,file="db5/measured_errors/Msubroutines_MS_MSS_MO_new.madx";
-call,file="db5/toolkit/Orbit_Routines.madx";
-call,file="slhc/errors2/SelectLHCMonCor.madx";
-readtable, file="db5/measured_errors/rotations_Q2_integral.tab";
-
-! Error routine and error table for new IT/D1/D2/Q4/Q5
-
-call,file="slhc/errors2/macro_error.madx";			! macros for error generation in the new IT/D1's
-
-if (chooseMQXF==1){
-	call,file="slhc/errors2/IT_errortable_v66_5";
-} else {
-	if (chooseMQXFAB==1) {
-		call,file="slhc/errors2/ITa_errortable_v65";
-		call,file="slhc/errors2/ITb_errortable_v65";
-	} else {
-		if (chooseMQXFbodyends==1) {
-			call,file="slhc/errors2/ITnc_errortable_v5"; ! target error table for the new IT
-			call,file="slhc/errors2/ITcs_errortable_v5"; ! target error table for the new IT
-			call,file="slhc/errors2/ITbody_errortable_v5"; ! target error table for the new IT
-		};
-	};
-};
-
-call,file="slhc/errors2/D1_errortable_v1_spec";! target error table for the new D1
-call,file="slhc/errors2/D2_errortable_v5_spec";! target error table for the new D2
-call,file="slhc/errors2/Q4_errortable_v2_spec";! target error table for the new Q4 in IR1 and IR5
-call,file="slhc/errors2/Q5_errortable_v0_spec";! target error table for the new Q5 in IR1 and IR5 and IR6
-
-!value,Rr_MQXCD,Rr_MBXAB,Rr_MBRD,Rr_MQYY,Rr_MQYL;
-!a2R_MQXCD_inj:=20.00;a2R_MQXCD_col:=20.00;! up to 3 mrad roll of the new MQX's (1 mrad r.m.s)
-
-!b5M_MQXCD_col  :=  0.0000 ; b5U_MQXCD_col  :=  0.4200 ; b5R_MQXCD_col  :=  0.4200 ; !..Errors as in IT_errortable_v3
-!a5M_MQXCD_col  :=  0.0000 ; a5U_MQXCD_col  :=  0.4300 ; a5R_MQXCD_col  :=  0.4300 ; !..Errors as in IT_errortable_v3
-!a6M_MQXCD_col  :=  0.0000 ; a6U_MQXCD_col  :=  0.3100 ; a6R_MQXCD_col  :=  0.3100 ; !..Errors as in IT_errortable_v3
-
-!call,file="fre/MCBX_errortable_v1";			! target error table for the new MCBX in IR1 and IR5
-call,file="slhc/errors2/MCBXFAB_errortable_v1";			! target error table for the new MCBX in IR1 and IR5
-if (ON_MBH==1){
-	call,file="slhc/errors2/MBH_errortable_v1";
-};
-
-
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-!				Correct signs of MCBXF reference field values
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-if (nominal_MCBXF_signs==1){
-on_x1save=on_x1; on_sep1save=on_sep1;
-on_x5save=on_x5; on_sep5save=on_sep5;
-on_x1=on_x1aux; on_sep1=on_sep1aux;
-on_x5=on_x5aux; on_sep5=on_sep5aux;
-
-sign_MCBXFAH.3L1 = ACBXH3.L1/(abs(ACBXH3.L1));
-sign_MCBXFAH.3R1 = ACBXH3.R1/(abs(ACBXH3.R1));
-sign_MCBXFAH.3L5 = ACBXH3.L5/(abs(ACBXH3.L5));
-sign_MCBXFAH.3L5 = ACBXH3.L5/(abs(ACBXH3.L5));
-sign_MCBXFAV.3L1 = ACBXV3.L1/(abs(ACBXV3.L1));
-sign_MCBXFAV.3R1 = ACBXV3.R1/(abs(ACBXV3.R1));
-sign_MCBXFAV.3L5 = ACBXV3.L5/(abs(ACBXV3.L5));
-sign_MCBXFAV.3L5 = ACBXV3.L5/(abs(ACBXV3.L5));
-
-sign_MCBXFBH.a2L1 = ACBXH1.L1/(abs(ACBXH1.L1));
-sign_MCBXFBH.a2R1 = ACBXH1.R1/(abs(ACBXH1.R1));
-sign_MCBXFBH.a2L5 = ACBXH1.L5/(abs(ACBXH1.L5));
-sign_MCBXFBH.a2L5 = ACBXH1.L5/(abs(ACBXH1.L5));
-sign_MCBXFBV.a2L1 = ACBXV1.L1/(abs(ACBXV1.L1));
-sign_MCBXFBV.a2R1 = ACBXV1.R1/(abs(ACBXV1.R1));
-sign_MCBXFBV.a2L5 = ACBXV1.L5/(abs(ACBXV1.L5));
-sign_MCBXFBV.a2L5 = ACBXV1.L5/(abs(ACBXV1.L5));
-
-sign_MCBXFBH.b2L1 = ACBXH1.L1/(abs(ACBXH1.L1));
-sign_MCBXFBH.b2R1 = ACBXH1.R1/(abs(ACBXH1.R1));
-sign_MCBXFBH.b2L5 = ACBXH1.L5/(abs(ACBXH1.L5));
-sign_MCBXFBH.b2L5 = ACBXH1.L5/(abs(ACBXH1.L5));
-sign_MCBXFBV.b2L1 = ACBXV1.L1/(abs(ACBXV1.L1));
-sign_MCBXFBV.b2R1 = ACBXV1.R1/(abs(ACBXV1.R1));
-sign_MCBXFBV.b2L5 = ACBXV1.L5/(abs(ACBXV1.L5));
-sign_MCBXFBV.b2L5 = ACBXV1.L5/(abs(ACBXV1.L5));
-value,on_x1,on_sep1,on_x5,on_sep5;
-value,sign_MCBXFAH.3L1,
-	  sign_MCBXFAH.3R1,
-	  sign_MCBXFAH.3L5,
-	  sign_MCBXFAH.3L5,
-	  sign_MCBXFAV.3L1,
-	  sign_MCBXFAV.3R1,
-	  sign_MCBXFAV.3L5,
-	  sign_MCBXFAV.3L5,
-	  sign_MCBXFBH.a2L1,
-	  sign_MCBXFBH.a2R1,
-	  sign_MCBXFBH.a2L5,
-	  sign_MCBXFBH.a2L5,
-	  sign_MCBXFBV.a2L1,
-	  sign_MCBXFBV.a2R1,
-	  sign_MCBXFBV.a2L5,
-	  sign_MCBXFBV.a2L5,
-	  sign_MCBXFBH.b2L1,
-	  sign_MCBXFBH.b2R1,
-	  sign_MCBXFBH.b2L5,
-	  sign_MCBXFBH.b2L5,
-	  sign_MCBXFBV.b2L1,
-	  sign_MCBXFBV.b2R1,
-	  sign_MCBXFBV.b2L5,
-	  sign_MCBXFBV.b2L5;
-on_x1=on_x1save; on_sep1=on_sep1save;
-on_x5=on_x5save; on_sep5=on_sep5save;
-};
-
-
-
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-!                Switch ON/OFF some multipole
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-exec ON_ALL_MULT;
-ON_LSF=1;
-
- ON_A1s =  0 ; ON_A1r =  0 ; ON_B1s =  0 ; ON_B1r =  0;
- ON_A2s =  1 ; ON_A2r =  1 ; ON_B2s =  0 ; ON_B2r =  0;
- ON_A3s =  1 ; ON_A3r =  1 ; ON_B3s =  1 ; ON_B3r =  1;
- ON_A4s =  1 ; ON_A4r =  1 ; ON_B4s =  1 ; ON_B4r =  1;
- ON_A5s =  1 ; ON_A5r =  1 ; ON_B5s =  1 ; ON_B5r =  1;
- ON_A6s =  1 ; ON_A6r =  1 ; ON_B6s =  1 ; ON_B6r =  1;
- ON_A7s =  1 ; ON_A7r =  1 ; ON_B7s =  1 ; ON_B7r =  1;
- ON_A8s =  1 ; ON_A8r =  1 ; ON_B8s =  1 ; ON_B8r =  1;
- ON_A9s =  1 ; ON_A9r =  1 ; ON_B9s =  1 ; ON_B9r =  1;
- ON_A10s =  1; ON_A10r =  1; ON_B10s =  1; ON_B10r =  1;
- ON_A11s =  1; ON_A11r =  1; ON_B11s =  1; ON_B11r =  1;
- ON_A12s =  ON_LSF; ON_A12r = ON_LSF; ON_B12s = ON_LSF; ON_B12r =  ON_LSF;
- ON_A13s =  ON_LSF; ON_A13r = ON_LSF; ON_B13s = ON_LSF; ON_B13r =  ON_LSF;
- ON_A14s =  ON_LSF; ON_A14r = ON_LSF; ON_B14s = ON_LSF; ON_B14r =  ON_LSF;
- ON_A15s =  ON_LSF; ON_A15r = ON_LSF; ON_B15s = ON_LSF; ON_B15r =  ON_LSF;
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-! Define the seed number (for statistical error assignment in the new IT/D1)
-! and call the corresponding measured error table for nominal LHC magnets
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-
-if (NRJ<4999.999) {readtable, file="wise/injection_errors-emfqcs-%seed_ran.tfs" ;};
-if (NRJ>5000.000) {readtable, file="wise/collision_errors-emfqcs-%seed_ran.tfs" ;};
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-!             apply field errors to MB magnets
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-call,file="db5/measured_errors/Efcomp_MB.madx";
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-!     correct orbit distortion resulting from MB magnets
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-if((ON_A1S)^2+(ON_A1r)^2+(ON_B1S)^2+(ON_B1r)^2 >0){
-  exec, initial_micado(4);
-  exec, initial_micado(4);
-};
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-!           apply field errors to all other magnets
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-! Nominal LHC magnets
-! Separation Dipoles
-call, file="db5/measured_errors/Efcomp_MBRB.madx";
-call, file="db5/measured_errors/Efcomp_MBRC.madx";
-call, file="db5/measured_errors/Efcomp_MBRS.madx";
-call, file="db5/measured_errors/Efcomp_MBX.madx" ;
-call, file="db5/measured_errors/Efcomp_MBW.madx" ;
-
-! Quadrupoles
-ON_B2Saux=on_B2S; on_B2S=0*ON_B2Saux; call, file="db5/measured_errors/Efcomp_MQW.madx" ; on_B2S=ON_B2Saux;
-ON_B2Saux=on_B2S; on_B2S=0*ON_B2Saux; call, file="db5/measured_errors/Efcomp_MQTL.madx"; on_B2S=ON_B2Saux;
-ON_B2Saux=on_B2S; on_B2S=0*ON_B2Saux; call, file="db5/measured_errors/Efcomp_MQMC.madx"; on_B2S=ON_B2Saux;
-ON_B2Saux=on_B2S; on_B2S=0*ON_B2Saux; call, file="db5/measured_errors/Efcomp_MQX.madx" ; on_B2S=ON_B2Saux;
-ON_B2Saux=on_B2S; on_B2S=0*ON_B2Saux; call, file="db5/measured_errors/Efcomp_MQY.madx" ; on_B2S=ON_B2Saux;
-ON_B2Saux=on_B2S; on_B2S=0*ON_B2Saux; call, file="db5/measured_errors/Efcomp_MQM.madx" ; on_B2S=ON_B2Saux;
-ON_B2Saux=on_B2S; on_B2S=0*ON_B2Saux; call, file="db5/measured_errors/Efcomp_MQML.madx"; on_B2S=ON_B2Saux;
-ON_B2Saux=on_B2S; on_B2S=0*ON_B2Saux; call, file="db5/measured_errors/Efcomp_MQ.madx"  ; on_B2S=ON_B2Saux;
-
-call,file="db5/measured_errors/Set_alignment_errors.madx";
-
-! New IT/D1/D2/Q4/Q5
-
-myseed=%seed_ran;
-
-/*
-eoption,seed=myseed+101;
-if (chooseMQXF==1) {
-	call, file="slhc/errors2/Efcomp_MQXF.madx";	
-} else {
-	if (chooseMQXFAB==1) {
-		call, file="slhc/errors2/Efcomp_MQXFA.madx";
-		call, file="slhc/errors2/Efcomp_MQXFB.madx";
-	} else {
-		if (chooseMQXFbodyends==1) {
-			call, file="slhc/errors2/Efcomp_MQXFbody.madx";
-			call, file="slhc/errors2/Efcomp_MQXFends.madx";
-		};
-	};
-};
-*/
-!eoption,seed=myseed+102;call, file="slhc/errors2/Efcomp_MBXAB.madx";     ! new D1 in IR1/5
-!eoption,seed=myseed+103;call, file="slhc/errors2/Efcomp_MBRD.madx";
-!eoption,seed=myseed+104;call, file="slhc/errors2/Efcomp_MQYY.madx";
-!eoption,seed=myseed+105;call, file="slhc/errors2/Efcomp_MQYL.madx";
-!eoption,seed=myseed+106;call, file="slhc/errors2/Efcomp_MCBXFAB.madx";     ! new triplet correctors in IR1/5
-ON_B2s = 0; ON_B2r = 0;
-if (ON_MBH==1){
-	eoption,seed=myseed+107;call, file="slhc/errors2/Efcomp_MBH.madx";
-};
-
-if (create_errorfiles_and_quit==1){
-	!exec show_error_newHLmagnet;
-	select,flag=error,clear;
-	select,flag=error,pattern=".";
-	esave,file="error_all.tfs";
-};
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-! MO settings
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-brho:=NRJ*1e9/clight;
-if (mylhcbeam==1){
-	KOF.A12B1:=Kmax_MO*I_MO/Imax_MO/brho; KOF.A23B1:=Kmax_MO*I_MO/Imax_MO/brho;
-	KOF.A34B1:=Kmax_MO*I_MO/Imax_MO/brho; KOF.A45B1:=Kmax_MO*I_MO/Imax_MO/brho;
-	KOF.A56B1:=Kmax_MO*I_MO/Imax_MO/brho; KOF.A67B1:=Kmax_MO*I_MO/Imax_MO/brho;
-	KOF.A78B1:=Kmax_MO*I_MO/Imax_MO/brho; KOF.A81B1:=Kmax_MO*I_MO/Imax_MO/brho;
-	KOD.A12B1:=Kmax_MO*I_MO/Imax_MO/brho; KOD.A23B1:=Kmax_MO*I_MO/Imax_MO/brho;
-	KOD.A34B1:=Kmax_MO*I_MO/Imax_MO/brho; KOD.A45B1:=Kmax_MO*I_MO/Imax_MO/brho;
-	KOD.A56B1:=Kmax_MO*I_MO/Imax_MO/brho; KOD.A67B1:=Kmax_MO*I_MO/Imax_MO/brho;
-	KOD.A78B1:=Kmax_MO*I_MO/Imax_MO/brho; KOD.A81B1:=Kmax_MO*I_MO/Imax_MO/brho;
-};
-
-if (mylhcbeam>1){
-	KOF.A12B2:=Kmax_MO*I_MO/Imax_MO/brho; KOF.A23B2:=Kmax_MO*I_MO/Imax_MO/brho;
-	KOF.A34B2:=Kmax_MO*I_MO/Imax_MO/brho; KOF.A45B2:=Kmax_MO*I_MO/Imax_MO/brho;
-	KOF.A56B2:=Kmax_MO*I_MO/Imax_MO/brho; KOF.A67B2:=Kmax_MO*I_MO/Imax_MO/brho;
-	KOF.A78B2:=Kmax_MO*I_MO/Imax_MO/brho; KOF.A81B2:=Kmax_MO*I_MO/Imax_MO/brho;
-	KOD.A12B2:=Kmax_MO*I_MO/Imax_MO/brho; KOD.A23B2:=Kmax_MO*I_MO/Imax_MO/brho;
-	KOD.A34B2:=Kmax_MO*I_MO/Imax_MO/brho; KOD.A45B2:=Kmax_MO*I_MO/Imax_MO/brho;
-	KOD.A56B2:=Kmax_MO*I_MO/Imax_MO/brho; KOD.A67B2:=Kmax_MO*I_MO/Imax_MO/brho;
-	KOD.A78B2:=Kmax_MO*I_MO/Imax_MO/brho; KOD.A81B2:=Kmax_MO*I_MO/Imax_MO/brho;
-};
-
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-!               correction of field errors in MB (compatible with V6.503 & SLHC)
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-select, flag=error, clear;
-select, flag=error, pattern=MB\.,class=multipole;
-if (ON_MBH==1){
-	select, flag=error, pattern=MBH\.,class=multipole;
-};
-esave,  file="temp/MB.errors";
-!system, "slhc/errors2/corr_MB_v3";
-system, "slhc/errors2/corr_MB_ats_v3";
-call,   file="temp/MB_corr_setting.mad";
-!exec reset_MB_corr;
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-! correction of triplet and D1 (only in collision, not compatible V6.503)
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-if (NRJ>5000.0000){
-	option, echo, warn, -info;
-	select, flag=error, clear;
-	select, flag=error, pattern=MQX,  class=multipole;
-	select, flag=error, pattern=MBX,  class=multipole;
-	select, flag=error, pattern=MBRC, class=multipole;
-	select, flag=error, pattern=MBRD, class=multipole;
-	if (correct_for_MCBX==1){
-		select, flag=error, pattern=MCBXF, class=multipole;
-	};
-	esave,  file="temp/tripD1D2.errors";
-	system, "slhc/errors2/corr_tripD1_v4";
-
-	call,file="temp/MCX_setting.mad";
-	!kcsx3.l1:=0 ; kcsx3.r1:=0 ; kcsx3.l5:=0 ; kcsx3.r5:=0 ; !switch off b3  correction in IR1 and IR5
-	!kcox3.l1:=0 ; kcox3.r1:=0 ; kcox3.l5:=0 ; kcox3.r5:=0 ; !switch off b4  correction in IR1 and IR5
-	!kcdx3.l1:=0 ; kcdx3.r1:=0 ; kcdx3.l5:=0 ; kcdx3.r5:=0 ; !switch off b5  correction in IR1 and IR5
-	!kctx3.l1:=0 ; kctx3.r1:=0 ; kctx3.l5:=0 ; kctx3.r5:=0 ; !switch off b6  correction in IR1 and IR5
-	!kqsx3.l1:=0 ; kqsx3.r1:=0 ; kqsx3.l5:=0 ; kqsx3.r5:=0 ; !switch off a2  correction in IR1 and IR5
-	!kcssx3.l1:=0; kcssx3.r1:=0; kcssx3.l5:=0; kcssx3.r5:=0; !switch off a3  correction in IR1 and IR5
-	!kcosx3.l1:=0; kcosx3.r1:=0; kcosx3.l5:=0; kcosx3.r5:=0; !switch off a4  correction in IR1 and IR5
-	!kcdsx3.l1:=0; kcdsx3.r1:=0; kcdsx3.l5:=0; kcdsx3.r5:=0; !switch off a5  correction in IR1 and IR5
-	!kctsx3.l1:=0; kctsx3.r1:=0; kctsx3.l5:=0; kctsx3.r5:=0; !switch off a6  correction in IR1 and IR5
-
-	kcsx3.l2 :=0; kcsx3.r2 :=0; kcsx3.l8 :=0; kcsx3.r8 :=0; !switch off b3  correction in IR2 and IR8
-	kcox3.l2 :=0; kcox3.r2 :=0; kcox3.l8 :=0; kcox3.r8 :=0; !switch off b4  correction in IR2 and IR8
-	kctx3.l2 :=0; kctx3.r2 :=0; kctx3.l8 :=0; kctx3.r8 :=0; !switch off b6  correction in IR2 and IR8
-	kqsx3.l2 :=0; kqsx3.r2 :=0; kqsx3.l8 :=0; kqsx3.r8 :=0; !switch off a2  correction in IR2 and IR8
-	kcssx3.l2:=0; kcssx3.r2:=0; kcssx3.l8:=0; kcssx3.r8:=0; !switch off a3  correction in IR2 and IR8
-	kcosx3.l2:=0; kcosx3.r2:=0; kcosx3.l8:=0; kcosx3.r8:=0; !switch off a4  correction in IR2 and IR8
-
-	!exec reset_MQX_corr;                             !switch off all IT multipole correction in all IR's
-};
-
-if (create_errorfiles_and_quit==1){
-	SELECT,FLAG=twiss,CLEAR;
-	SELECT,FLAG=twiss,COLUMN=name,s,k0l,k0sl,k1l,k1sl,k2l,k2sl,k3l,k3sl,k4l,k4sl,k5l,k5sl,k6l,k6sl,k7l,k7sl,k8l,k8sl,k9l,k9sl,k10l,k10sl,k11l,k11sl,k12l,k12sl,k13l,k13sl,k14l,k14sl,k15l,k15sl,k16l,k16sl,k17l,k17sl,k18l,k18sl,k19l,k19sl,k20l,k20sl,angle,hkick,vkick;
-	TWISS,SEQUENCE=lhcb2,save,centre,FILE=./multipoles.tfs;
-	stop;
-};
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-!    correct orbit distortion resulting from other magnets
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-if((ON_A1S)^2+(ON_A1r)^2+(ON_B1S)^2+(ON_B1r)^2 >0){
-	exec, initial_micado(4);
-	exec, initial_micado(4);
-	exec, initial_micado(4);
-	exec, initial_micado(4);
-};
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-!    final orbit correction before applying crossing scheme
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-if((ON_A1S)^2+(ON_A1r)^2+(ON_B1S)^2+(ON_B1r)^2 >0){
-	exec, final_micado(0.004);
-};
-
-! end of orbit correction, now switch ON crossing scheme
-! restoring crossing angle
-on_x1=on_x1aux; on_sep1=on_sep1aux;
-on_x2=on_x2aux; on_sep2=on_sep2aux;
-on_x5=on_x5aux; on_sep5=on_sep5aux;
-on_x8=on_x8aux; on_sep8=on_sep8aux; on_a8=on_a8aux;
-on_alice=on_aliceaux; on_lhcb=on_lhcbaux;
-on_disp=on_dispaux;
-
-! coguess,x=x.ip1,px=px.ip1,y=y.ip1,py=py.ip1;
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-!   Fine tuning of coupling after CO correction and with Xscheme
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-call,file="slhc/errors2/FineCouplingCorrection.madx";
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-!                 matching of orbit, tune and chromaticity
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-if(ON_COLLISION==0) {ON_BB_CHARGE:=1;};    ! W/o head-on Q and Q' are matched with bb
-
-!Rematch the Xscheme towards specified separation and Xange in IP1/2/5/8
-call,file="slhc/toolkit/rematchCOIP.madx";
-
-!Rematch the CO in the arc for dispersion correction
-if(ON_DISP<>0) {call,file="slhc/toolkit/rematchCOarc.madx";};
-
-if(mylhcbeam==1){
-	kqtf=kqtf.b1;
-	kqtd=kqtd.b1;
-	kqtf.b1:=kqtf;
-	kqtd.b1:=kqtd;
-
-	ksf=ksf.b1;
-	ksd=ksd.b1;
-	ksf.b1:=ksf;
-	ksd.b1:=ksd;
-};
-if(mylhcbeam>1){
-	kqtf=kqtf.b2
-	kqtd=kqtd.b2
-	kqtf.b2:=kqtf
-	kqtd.b2:=kqtd
-
-	ksf=ksf.b2
-	ksd=ksd.b2
-	ksf.b2:=ksf
-	ksd.b2:=ksd;
-};
-/*
-match;
-global, q1=qx0, q2=qy0;
-vary,   name=kqtf, step=1.0E-7 ;
-vary,   name=kqtd, step=1.0E-7 ;
-lmdif,  calls=100, tolerance=1.0E-12;
-endmatch;
-
-match,chrom;
-global, dq1=qprimex, dq2=qprimey;
-vary,   name=ksf;
-vary,   name=ksd;
-lmdif,  calls=100, tolerance=1.0E-12;
-endmatch;
-
-match,chrom;
-global, dq1=qprimex, dq2=qprimey;
-global, q1=qx0, q2=qy0;
-vary,   name=ksf;
-vary,   name=ksd;
-vary,   name=kqtf, step=1.0E-7 ;
-vary,   name=kqtd, step=1.0E-7 ;
-lmdif,  calls=500, tolerance=1.0E-12;
-endmatch;
-*/
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-!                 Check corrector strength
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-call, file="slhc/errors2/corr_value.madx";
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-!                 final twiss before sending to sixtrack
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-ON_BB_CHARGE:=1;
-
-!if (NRJ<4999.9999) {VRF400:=%rf_vol ;LAGRF400.B1=0.5; LAGRF400.B2=0.;};
-!if (NRJ>5000.0000) {VRF400:=%rf_vol ;LAGRF400.B1=0.5; LAGRF400.B2=0.;};
-
-VRF400:=%rf_vol;
-LAGRF400.B1=0.5;
-LAGRF400.B2=0.;
-
-twiss;
-!emit;
-xnom1=table(twiss,IP1,x); pxnom1=table(twiss,IP1,px); ynom1=table(twiss,IP1,y); pynom1=table(twiss,IP1,py);
-xnom2=table(twiss,IP2,x); pxnom2=table(twiss,IP2,px); ynom2=table(twiss,IP2,y); pynom2=table(twiss,IP2,py);
-xnom5=table(twiss,IP5,x); pxnom5=table(twiss,IP5,px); ynom5=table(twiss,IP5,y); pynom5=table(twiss,IP5,py);
-xnom8=table(twiss,IP8,x); pxnom8=table(twiss,IP8,px); ynom8=table(twiss,IP8,y); pynom8=table(twiss,IP8,py);
-
-sixtrack, cavall, radius=0.017;
-
-stop;
diff --git a/templates/init_cartesian_dist.py b/templates/init_cartesian_dist.py
new file mode 100644
index 0000000..2b8739c
--- /dev/null
+++ b/templates/init_cartesian_dist.py
@@ -0,0 +1,41 @@
+#!/usr/bin/env python
+import argparse
+import math
+
+def linspace(a, b, n):
+    '''Numpyless linear spacing function.
+    '''
+    if n < 2:
+        return a
+    diff = (float(b) - a)/(n - 1)
+    return [diff * i + a for i in range(n)]
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser()
+    parser.add_argument('ampx0', type=float)
+    parser.add_argument('ampx1', type=float)
+    parser.add_argument('ampy0', type=float)
+    parser.add_argument('ampy1', type=float)
+    parser.add_argument('number', type=int)
+    args = parser.parse_args()
+
+    num_part_div = int(math.sqrt(args.number))
+    amp_x = linspace(args.ampx0, args.ampx1, num_part_div+1)
+    amp_y = linspace(args.ampy0, args.ampy1, num_part_div+1)
+    offset = (args.ampx1 - args.ampx0) / (2.*float(num_part_div))
+    amp_x = [x + offset for x in amp_x[:-1]]
+    amp_y = [y + offset for y in amp_y[:-1]]
+    with open('input_dist.txt', 'w') as fp:
+        for i, ampx in enumerate(amp_x):
+            for j, ampy in enumerate(amp_y):
+                rx = math.sqrt(ampx)
+                ry = math.sqrt(ampy)
+
+                # particle IDs
+                pid = 2*(i*len(amp_y) + j)
+                pid_twin = pid + 1
+                fp.write('{:d} {:.10f} {:.10f}\n'.format(
+                    pid, rx, ry))
+                fp.write('{:d} {:.10f} {:.10f}\n'.format(
+                    pid_twin, rx+1e-6, ry+1e-6))
diff --git a/templates/init_polar_dist.py b/templates/init_polar_dist.py
new file mode 100644
index 0000000..406786e
--- /dev/null
+++ b/templates/init_polar_dist.py
@@ -0,0 +1,34 @@
+#!/usr/bin/env python
+import argparse
+import math
+
+def linspace(a, b, n):
+    '''Numpyless linear spacing function.
+    '''
+    if n < 2:
+        return a
+    diff = (float(b) - a)/(n - 1)
+    return [diff * i + a for i in range(n)]
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser()
+    parser.add_argument('amp0', type=float)
+    parser.add_argument('amp1', type=float)
+    parser.add_argument('angle', type=float)
+    parser.add_argument('number', type=int)
+    args = parser.parse_args()
+
+    amps = linspace(args.amp0, args.amp1, args.number)
+    with open('input_dist.txt', 'w') as fp:
+        for i, amp in enumerate(amps):
+            rx = math.sqrt(amp * math.cos(args.angle))
+            ry = math.sqrt(amp * math.sin(args.angle))
+
+            # particle IDs
+            pid = 2*i
+            pid_twin = pid + 1
+            fp.write('{:d} {:.10f} {:.10f}\n'.format(
+                pid, rx, ry))
+            fp.write('{:d} {:.10f} {:.10f}\n'.format(
+                pid_twin, rx+1e-6, ry+1e-6))
diff --git a/tests/unit/test_study_params.py b/tests/unit/test_study_params.py
index 027c1c6..0b5cd58 100644
--- a/tests/unit/test_study_params.py
+++ b/tests/unit/test_study_params.py
@@ -1,6 +1,7 @@
 import shutil
 import unittest
 import sys
+from math import pi, cos, sqrt, sin
 from pathlib import Path
 # give the test runner the import access
 pysixdesk_path = str(Path(__file__).parents[2].absolute())
@@ -237,5 +238,21 @@ def test_product_dict(self):
                                {'a': 2, 'b': 4},
                                ])
 
+    def test_da_angle(self):
+        # make sure the angle calculation is correct
+        params = StudyParams(mask_path=self.mask_file,
+                             fort_path=self.fort_file)
+
+        # make sure default angles are ~correct
+        self.assertEqual(len(params['angle']), 7)
+
+        # one angle between ]0, pi/2[ --> pi/4
+        angles = params.da_angles(start=0, end=pi/2, n=1)
+        self.assertEqual(angles, [pi/4])
+
+        # 3 angles between ]0, pi/2[ --> pi/8 increments
+        angles = params.da_angles(start=0, end=pi/2, n=3)
+        self.assertEqual(angles, [(i + 1) * pi / 8 for i in range(3)])
+
     def tearDown(self):
         shutil.rmtree(self.test_folder.parent, ignore_errors=True)
diff --git a/tests/unit/test_utils.py b/tests/unit/test_utils.py
index 6b7ed94..f57deef 100644
--- a/tests/unit/test_utils.py
+++ b/tests/unit/test_utils.py
@@ -56,20 +56,6 @@ def test_compress_buf(self):
 
         # with gzip ...
 
-    def test_merge_dicts(self):
-        a = {'a': 1, 'b': 2, 'c': 3}
-        b = {'a': 10}
-        self.assertEqual(utils.merge_dicts(a, b), {'a': 10, 'b': 2, 'c': 3})
-        # also check to make sure the dicts have not been changed
-        self.assertEqual(a, {'a': 1, 'b': 2, 'c': 3})
-        self.assertEqual(b, {'a': 10})
-
-        b = {'d': 10}
-        self.assertEqual(utils.merge_dicts(a, b), {'a': 1, 'b': 2,
-                                                   'c': 3, 'd': 10})
-        self.assertEqual(a, {'a': 1, 'b': 2, 'c': 3})
-        self.assertEqual(b, {'d': 10})
-
     def test_concatenate_files(self):
         utils.concatenate_files([self.concat_file_in_1, self.concat_file_in_2],
                                 self.concat_file_out)