import json import copy import numpy as np from PyQt5 import QtWidgets,QtGui from onlinemodel.core import Facility from onlinemodel.madx import CMadX converttwiss= {'betx':'betx','bety':'bety','alfx':'alfx','alfy':'alfy', 'etax':'dx','etay':'dy','etapx':'dpx','etapy':'dpy','mux':'mux','muy':'muy', 'x':'x','y':'y','px':'px','py':'py','energy':'energy'} maxgradient= {'QFF':50.,'QFM':50,'QFD':20,'QFS':20,'QFDM':20,'QFA':31.75,'HFA':37.,'HFB':1445} class Model: def __init__(self, phase=0, parent=None): print('Initializing online model ...') self.phase = phase # current planned future self.parent=parent self.om = Facility(init=1, alt = phase) self.order = None self.madx = CMadX() self.startTwiss = None self.startEnergy = None self.energyReference ='SINLH02.MBND100' # flag to enfore new lattice self.forceLat=True # hook up events self.eventHandling() def eventHandling(self): self.parent.UITrack.clicked.connect(self.track) self.parent.UIReportMagnetStrength.clicked.connect(self.checkMagnetLimit) def forceLatticeUpdate(self): self.forceLat=True def getLatticeVersion(self): return self.om.Version def getInitialEnergy(self): return self.om.EnergyAt(self.energyReference)[0] def initializeMagnets(self): print('Initializing all magnets to zero') self.om.setRegExpElement('.*', 'MQUA', 'k1', 0) self.om.setRegExpElement('.*', 'MQSK', 'k1', 0) self.om.setRegExpElement('.*', 'MSEX', 'k2', 0) def updateFromMachine(self,machine): pol = ['LH', 'LV+', 'LV-', 'C+', 'C-', 'ZL'] mag = machine['Magnet'] for key in mag: keyom=key.replace('-','.') self.updateElement(keyom,[mag[key]]) # needs to be a list und = machine['Undulator'] for key in und: keyom=key.replace('-','.') self.updateElement(keyom,[und[key][0],pol[int(und[key][1])]]) rf = machine['RF'] for key in rf: keyom=key.replace('-','.') self.updateElement(keyom,rf[key]) energy = machine['Energy'] for key in energy: self.updateEnergy(energy[key]) self.forceLatticeUpdate() def updateEnergy(self,E0): if isinstance(E0,list): E0=E0[0] self.om.forceEnergyAt(self.energyReference, E0*1e6) def updateElement(self,name,val): if 'MQUA' in name or 'MQSK' in name: L = self.om.getRegExpElement(name[0:7], name[8:15],'Length')[0] self.om.setRegExpElement(name[0:7], name[8:15], 'k1', float(val[0])/L) if 'MSEX' in name: L = self.om.getRegExpElement(name[0:7], name[8:15], 'Length')[0] self.om.setRegExpElement(name[0:7], name[8:15], 'k2', float(val[0])/L) if 'MBND' in name: self.om.setRegExpElement(name[0:7], 'MBND', 'angle', float(val[0])) if 'UMOD' in name: self.om.setRegExpElement(name[0:7], 'UMOD', 'K', float(val[0])) if 'UIND' in name: self.om.setRegExpElement(name[0:7], name[8:15], 'K', float(val[0])) if 'SATUN' in name: kx = 0 ky = 1 if 'LV' in val[1]: kx = 1 ky = 0 if 'C' in val[1]: kx = 0.5 ky=0.5 self.om.setRegExpElement(name[0:7], name[8:15], 'kx', kx) self.om.setRegExpElement(name[0:7], name[8:15], 'ky', ky) if 'RSYS' in name: grad = float(val[0]) phase = float(val[1]) if 'CB' in name[0:7]: grad = grad/ 4. elif 'XB' in name[0:7] or 'SINSB03' in name or 'SINSB04' in name: grad = grad/2 L = self.om.getRegExpElement(name[0:7], 'RACC', 'Length')[0] self.om.setRegExpElement(name[0:7], 'RACC', 'Gradient', grad/L) self.om.setRegExpElement(name[0:7], 'RACC', 'Phase', phase) def getElements(self): return self.om.listElement('*', 1) def getSettings(self): elements = self.getElements() quadrupoles={} sextupoles={} dipoles={} rf={} undulators={} kicker={} energy={'location': self.energyReference, 'energy' : 1e-6*self.om.EnergyAt(self.energyReference)[0]} for ele in elements: if 'MQUA' in ele.Name or 'MQSK' in ele.Name: quadrupoles[ele.Name]={'k1':ele.k1,'k1L':ele.k1*ele.Length} elif 'MSEX' in ele.Name: sextupoles[ele.Name]={'k2':ele.k2,'k2L':ele.k2*ele.Length} elif 'MBND' in ele.Name: if 'SINLH' in ele.Name or 'SINBC' in ele.Name or 'S10BC' in ele.Name or 'SATMA' in ele.Name or 'SATUN' in ele.Name: dipoles[ele.Name]={'angle':ele.angle} elif 'UIND' in ele.Name or 'UMOD' in ele.Name: undulators[ele.Name]={'K':ele.K,'kx':ele.kx,'ky':ele.ky} elif 'RACC' in ele.Name: rf[ele.Name]={'Gradient':ele.Gradient*ele.Length,'Phase':ele.Phase} elif 'MKAC' in ele.Name or 'MKDC' in ele.Name: kicker[ele.Name] = {'cory': ele.cory,'design_kick':ele.design_kick} return {'Quadrupole':quadrupoles,'Sextupole':sextupoles,'Dipole':dipoles,'RF':rf,'Undulator':undulators, 'Kicker':kicker,'Energy':energy, 'InitialCondition':self.startTwiss,'Phase':self.phase} def loadSettingsGroup(self,group,fields,normalized=False): for key in group.keys(): ele = self.om.getElement(key) for field in fields: ele.__dict__[field]=group[key][field] if normalized: ele.__dict__[field]/=ele.Length def loadSettings(self,settings): if not 'Phase' in settings.keys(): return False if not settings['Phase'] == self.phase: return False self.loadSettingsGroup(settings['Quadrupole'],['k1']) self.loadSettingsGroup(settings['Sextupole'], ['k2']) self.loadSettingsGroup(settings['Dipole'], ['angle']) self.loadSettingsGroup(settings['RF'], ['Gradient'],normalized=True) self.loadSettingsGroup(settings['RF'], ['Phase']) self.loadSettingsGroup(settings['Undulator'], ['K','kx','ky']) self.loadSettingsGroup(settings['Kicker'], ['cory','design_kick']) self.startEnergy = settings['Energy']['energy'] self.energyReference = settings['Energy']['location'] self.startTwiss = settings['InitialCondition'] self.updateEnergy(self.startEnergy) print('Settings loaded (Reference Energy:',self.om.EnergyAt(self.energyReference)[0],')') self.forceLat=True return True def updateModelFromMatching(self,var): for magm in var.keys(): mag0 = magm[0:15] val = var[magm] ele = self.om.getElement(mag0) if 'mqua' in magm or 'mqsk' in magm: if not ele is None: ele.k1 = val print('Updating %s to k1: %8.4f' % (mag0, val)) magnets=self.parent.reference.getDependence(mag0) if not magnets is None: for magd in magnets: ele = self.om.getElement(magd) if not ele is None: ele.k1 = val print('Updating %s to k1: %8.4f' % (magd, val)) elif 'mkac' in magm or 'mkdc' in magm: if not ele is None: ele.design_kick = val ele.cory = val print('Updating %s to design_kick: %8.4f' % (mag0, val)) ################## # tracking def track(self): start = str(self.parent.UITrackStart.text()).upper() end = str(self.parent.UITrackEnd.text()).upper() if len(start)>7: start = start[0:7] if len(end)>7: end = end[0:7] refloc, twiss0 = self.parent.reference.getReference() if 'SWISSFEL' in refloc.upper(): refloc = 'SINLH01$START' if refloc.upper() == 'START': refloc = start.upper() if self.startEnergy is None: self.startEnergy = self.om.EnergyAt(self.energyReference)[0]*1e-6 twiss0['energy'] = self.startEnergy*1e-3 # convert to GeV for madx start, end = self.checkRange(start, end, refloc[0:7]) if start is None: return print('Tracking from',start,'to',end) self.doTrack(start,end,refloc,twiss0) def doTrack(self, start,end, refloc, twiss_in, plot = True): print(twiss_in) twiss0 = {converttwiss[key]:twiss_in[key] for key in twiss_in.keys()} # follows the madx format self.setBranch(end.upper()) print("Tracking:",refloc[0:7],start) if not refloc[0:7] == start: twiss0 = self.doBackTrack(refloc,start,twiss0) # track backwards to get new initial twiss values twiss = self.madx.track('swissfel',start+'$START',end+'$END',twiss0) energy = self.calcEnergyProfile(twiss) if plot: self.parent.plot.newData(twiss,energy) def doBackTrack(self,start=None,end=None,twiss0=None): print('Backtracking 1:',twiss0) twiss0['alfx'] = -twiss0['alfx'] # revert particle trajectories twiss0['alfy'] = -twiss0['alfy'] twiss0['dpx'] = -twiss0['dpx'] twiss0['dpy'] = -twiss0['dpy'] twiss0['px'] = -twiss0['px'] twiss0['py'] = -twiss0['py'] # self.madx.updateVariables(twiss0) # print(start,end) twiss = self.madx.track('invswissfel', start, end+'$START',twiss0) twiss1 = {} twiss1['betx'] = twiss.betx[-1] # revert trajectory back twiss1['bety'] = twiss.bety[-1] twiss1['alfx'] = -twiss.alfx[-1] twiss1['alfy'] = -twiss.alfy[-1] twiss1['dx'] = twiss.dx[-1] twiss1['dy'] = twiss.dy[-1] twiss1['dpx'] = -twiss.dpx[-1] twiss1['dpy'] = -twiss.dpy[-1] twiss1['x'] = twiss.x[-1] twiss1['y'] = twiss.y[-1] twiss1['px'] = -twiss.px[-1] twiss1['py'] = -twiss.py[-1] return twiss1 def calcEnergyProfile(self,twiss): energy = np.array([0. for i in range(len(twiss.betx))]) e0 = 0. for i, name in enumerate(twiss.name): if len(name) > 15: elename = name[0:15] ele = self.om.getElement(elename) if not ele is None: erg= self.om.EnergyAt(ele) if e0 == 0.: energy[:i]+=erg[0]*1e-6 e0 = (erg[0]+erg[1])*1e-6 energy[i]=e0 return energy def setBranch(self,end): destination = 'ARAMIS' if 'SPO' in end: destination = 'PORTHOS' elif 'SAT' in end: destination = 'ATHOS' elif 'S10BD' in end or 'SIN' in end: destination = 'INJECTOR' self.om.setBranch(destination,'SINLH01') self.order=self.om.getBranchElements() self.madx.updateLattice(self.om,destination,self.forceLat) self.forceLat=False def checkRange(self,start,end,ref): if self.om.getSection(start) is None: print('Invalid staring point for tracking. Setting to SINLH01') startsec='SINLH01' else: startsec = self.om.getSection(start).Element[0].Name if self.om.getSection(end) is None: print('Invalid staring point for tracking. Setting to starting section') endsec = startsec else: endsec = self.om.getSection(end).Element[0].Name if self.om.getSection(ref) is None: print('Invalid section for reference. Aborting tracking') return None,None refsec = self.om.getSection(ref).Element[0].Name if not self.om.isUpstream(startsec,endsec): return self.checkRange(end,start,ref) if not self.om.isUpstream(startsec,refsec): start = ref print('Reference point is upstream the tracking range. Extending range') if not self.om.isUpstream(refsec,endsec): end = ref print('Reference point is downstream the tracking range. Extending range') return start,end ########################################3 def checkMagnetLimit(self): # using magnet calibration self.parent.UIReportMagnetResult.clear() self.updateEnergy(140.) quad={} quad['QFA'] = [0.15, 0.173, 150.0, 91.1, 0.0, 0.74521, -0.00813, -0.03542, 22.5e-3] # From MEASUREMENT. Documentation: https://intranet.psi.ch/pub/Swiss_FEL/FinQuadrupoles/QFA.pdf quad['QFB'] = [0.08, 0.11, 10.0, 4.9, 0.0, 0.03815, -1.8e-4, -0.7e-3, 22.5e-3] # From MEASUREMENT(?). Documentation(Implementaion): the matlab version Injector on-line model, Curr2KLquad.m quad['QFC'] = [0.08, 0.11, 10.0, 10.1, 0.0, 0.5, 0.0, 0.0, 1.0] # From DESIGN(?). Documentation: Injector wiki, Subsystems (as of 14th Aug 2014) quad['QFD'] = [0.15, 0.1628, 10.0, 5.6, 0.0, 0.23313, -27.6e-4, 15.5e-4, 0.011] # From Measurement. Documentation: FEL-SS88-007-7.pdf (found in Alfresco: Company Home > Projects > SwissFEL > Facility > 8850 Magnets) quad['QFDM'] = quad['QFD'] # same type but weaker windings for corrector quad['QFF'] = [0.08, 0.0875, 10.0, 2.9, 0.0, 0.32897, -0.5e-4, -42.9e-4, 0.006] # From Measurement. Documentation: FEL-SS88-008-8.pdf (found in Alfresco) quad['QFM'] = [0.3, 0.311, 50.0, 21.1, 0.0, 0.64541, -0.00296, -0.00617, 0.011] # From Measurement k1brho={} for key in quad.keys(): [PhysicalLength, MagneticLength, Imax, Ilin, b0, b1, b2, b3, R] = quad[key] # quick conversion for quadrupoles and sextupoles if abs(Imax) < Ilin: B = b0 + b1 else: B = b0 + b1 + b2 +b3 k1brho[key]= B * MagneticLength/R/PhysicalLength #k1 = val/brho elements=self.getElements() for ele in elements: if 'MQUA' in ele.Name.upper(): bg = ele.Baugruppe Energy = self.om.EnergyAt(ele.Name)[0]*1e-6 gamma = (Energy + 0.511) / 0.511 # Energy in units of MeV since 29.06.2015 beta = np.sqrt(1 - 1 / gamma ** 2) p = 511e3 * beta * gamma / 1e9 # momentum in GeV/c brho = p / 0.299792 # P[GeV/c]=0.3*Brho[T.m] if bg in k1brho.keys(): k1max=k1brho[bg]/brho rat=np.abs(ele.k1/k1max)*100. label = '%s: %5.1f %%' % (ele.Name.replace('.','-'),rat) #print(label) listitem = QtWidgets.QListWidgetItem(label) color = QtGui.QColor(255,255,255) # white if rat>95: color = QtGui.QColor(255,255,0) if rat>100: color = QtGui.QColor(255,0,0) listitem.setBackground(color) self.parent.UIReportMagnetResult.addItem(listitem) else: print('%s: Unsupported type %s' % (ele.Name.replace('.','-'),bg))