Prototype LHC Impedance Database

John Jowett

textDate[]

Introduction

Following the ABP-RLC meeting " 7/11/2006", this notebook takes a quick look at how one could implement the LHC impedance database in Mathematica.

Note that this is a development notebook which contains rather more "code" (trial function definitions etc.) than an eventual application notebook would normally contain.  So there is no need to be put off by this.  A more "graphical" user interface (with buttons, etc.) can also be created.   However full flexibility will be retained with access to the package functions.

Once it is clear what functions will be useful, they can be tidied up to form a Zbase package that will be loaded with a single command.  This is basically just a matter of moving the definition cells to a package.

So far this notebook implements access to the data, reads in some typical examples and prototypes functions for plotting and manipulating them.  It also shows how the impedance database can be linked to the LHC optics to provide the information about numbers of element occurrences, their locations, average β-functions, etc.   Thus there is no need to duplicate this information in the impedance database.

Documentation of the previous implementation:

http://bruening.home.cern.ch/bruening/ZBASE/zbase.html

I have copied key pieces of text from the main documentation where useful in the following sections.

Developing functions in the global context for now.   Later they should be put in a proper package context.

Since this is a brand-new application package, it might be worth switching to Wolfram Workbench as development environment.

Setup

This notebook should work on NICE or CERN Linux systems.

Load packages

Some of the standard Mathematica packages that are useful in this application

Needs["Graphics`Graphics`"]

Needs["Statistics`DataSmoothing`"]

This notebook needs some of my other packages.  The following cell provides access to these on CERN computers (if you don't have it already).

$jowettPackagesDirectory=ToFileName[{If[$OperatingSystem=="Windows","G:\\Users\\j\\jowett","/afs/cern.ch/user/j/jowett"],"public","math"}];
If[Not[MemberQ[$Path,$jowettPackagesDirectory]],PrependTo[$Path,$jowettPackagesDirectory]]

Utility package that facilitates access to CERN file systems and other things:

Needs["JMJUtils`BumpEtc`"]

Check that we have access to AFS system for the data (on Windows, this requires the AFS client):

$AfsRoot

General package with constants, units, dimension-checking, etc.

Needs["Accelerator`ConstantsUnits`"]

This is handy for graphics of mode spectra

Needs["MathMethods`SpectralPlot`"]

Access to LHC optics so that we know the elements installed in the rings, etc.

Needs["Madtomma`LHCapplications`MADinputLHC`"]

Needs["Madtomma`Mfs`OpticsUtilities`"]

Options and other settings

$TextStyle={FontFamily->"Times",FontSize->12}

Working Directory

Define a working directory and go there (this will give a harmless message if the directory already exists).

workDir=ThisNotebookDirectory[]

CreateDirectory[workDir];
SetDirectory[workDir]

Existing Database

This section has a first shot at developing functions to return information about the existing Zbase data structure.

Database data structure

Assume that we keep the structure from the previous implementation for now where impedances are found in a directory tree.

Functions for accessing the data structure

First some things beginning with $ that set up the environment.  They will not change often.

$ZbaseDirectory::usage="$ZbaseDirectory returns the path of the main Zbase main directory.";

We could also consider the possibility of Web access, making this a URL.   It would be easy to import data files but I'm not sure how to implement some of the following functions (we need a FileNames function for the Web?).

ZbaseDataDirectory::usage="ZbaseDataDirectory[] returns the path to the Zbase data directory.";

$ZbaseModes::usage="$ZbaseModes is the list of possible \"modes\" of impedance calculations, occurring as the lowest sub-directory names in the Zbase data structure.";

$ZbaseModes={"long","tran"};

This should give the right path on CERN Linux or Windows (with the AFS client)

$ZbaseDirectory=ToFileName[{$AfsRoot,"cern.ch","eng","lhc","impedance","zbase"}]

ZbaseDataDirectory[]:=ToFileName[{$ZbaseDirectory,"data"}]

Try it out:

ZbaseDataDirectory[]

To avoid having global variables, it is preferable to specify what we are working on as options of a function.  At the moment, I'm not sure what this function should do, but let's call it Zbase.  Maybe it's role is just to supply options to other functions.  The options can be set to default values to facilitate working on a particular machine, a particular item within it, etc.

Zbase::usage="Zbase[] is a function implementing the impedance database.  It's options are used by other functions.";

ZbaseMachine::usage="ZbaseMachine is an option of Zbase defining the machine.";

ZbaseItem::usage="ZbaseItem is an option of Zbase defining the item (it may be a list).";

ZbaseMode::usage="ZbaseMode is an option of Zbase defining the impedance \"mode\", usually either \"long\" or \"tran\".  The possibilities are stored in $ZbaseModes.";

Define the default options.  This is arbitrary for now.

Options[Zbase]={ZbaseMachine->"lhc",ZbaseItem->{"cavities","scavity"},ZbaseMode->First[$ZbaseModes]}

Zbase[]:=Null

The ZbaseItem and ZbaseMode options may not be necessary.  Leave them for now.

Implementation of functions

ZbaseMachines::usage="ZbaseMachines[] returns a list of all machines stored in the Zbase database.";

ZbaseMachines[]:=Module[{fns},
SetDirectory[ZbaseDataDirectory[]];
fns=Cases[{#,FileType[#]}&/@FileNames["*"],{_,Directory}];
ResetDirectory[];
First/@fns]

ZbaseItems::usage="ZbaseItems[itemClass] returns a list of all the items within a class as lists of strings.  With no argument, it returns all the items.";

ZbaseItems[itemClass:_String|{_String..}:".",opts___?OptionQ]:=Module[{machine,items},
machine=ZbaseMachine/.{opts}/.Options[Zbase];
SetDirectory[ToFileName[Flatten[{ZbaseDataDirectory[],machine,itemClass}]]];
items=Cases[{#,FileType[#]}&/@FileNames["*","*",2],{_,Directory}];
ResetDirectory[];
Union[(StringSplit[First[#],$PathnameSeparator]&/@items)/.{a__,Alternatives@@$ZbaseModes}->{a}]
]

ZbaseFiles::usage="ZbaseFiles[itemClass] returns all the impedance data files corresponding to an item class.";

ZbaseFiles[itemClass_String:"."]:=ZbaseFiles[{itemClass}];

ZbaseFiles[itemClass:{_String..},mode:Alternatives@@$ZbaseModes,opts___?OptionQ]:=Module[{machine,items,files},
machine=ZbaseMachine/.{opts}/.Options[Zbase];
SetDirectory[ToFileName[Flatten[{ZbaseDataDirectory[],machine,itemClass,mode}]]];
files=Cases[{#,FileType[#]}&/@FileNames[],{_,File}];
ResetDirectory[];
files
]

?Zbase*

Examples

ZbaseMachines[]

ZbaseItems[]

Length/@%

ZbaseItems[{"experiments"}]

Set an option to be different from the default

ZbaseItems[".",ZbaseMachine->"lep"]//TableForm

FileNames[]

ZbaseFiles[{"bellows-us","collimator"},"long"]//TableForm

Analogues of other Zbase environment variables

Some were specific to the existing implementation.

Do we need ZbaseEditor ?  Possibly, if we are going to launch text editing windows, say, to prepare input files for other programs.

Contents of impedance data files

Only showing how to get the data out of the files.  Presumably we would wrap this information up in suitable functions that would probably return Mathematica rules.

Contents of the .zbase file

Pick an example of this kind of file

zbaseFile=ToFileName[{ZbaseDataDirectory[],"lhc","cavities","scavity","long"},".zbase"]

Importing the file as a Table is an easy way to deal with the E-format numbers.

zbaseFileContents=Import[zbaseFile,"Table"];

TableForm[zbaseFileContents]

There is a small problem with the string in the first line, to be taken care of.

zbase[[1]]//InputForm

Example of how we might represent this information in Mathematica:

zbaseFileRules={"COMMENT"->StringJoin@@ToString[zbaseFileContents[[1]]],NumberOfItems->zbaseFileContents[[2]],
betxav->zbaseFileContents[[3]],
betyav->zbaseFileContents[[4]],
Zovern->zbaseFileContents[[5]],
HOMs->Drop[zbaseFileContents,10]}

InputForm[zbaseFileContents[[8]]]

Show the first few HOMs

Take[HOMs/.zbaseFileRules,5]

Plots, to be wrapped up as functions with proper option inheritance etc.

SpectralListPlot[Drop[HOMs/.zbaseFileRules,1]/.{f_,rq_,_}->{f/(10^9),rq},PlotRange->{{0,Automatic},{0,6}},PlotStyle->Red,AxesLabel->{"f/GHz","R/Q /Ω"}]

SpectralListPlot[HOMs/.zbaseFileRules/.{f_,rq_,_}->{f/(10^9),rq},PlotRange->{{0,Automatic},{0,6}},PlotStyle->Red,AxesLabel->{"f/GHz","R/Q /Ω"}]

SpectralListPlot[HOMs/.zbaseFileRules/.{f_,rq_,q_}->{f,q},PlotRange->{{0,Automatic},{0,Automatic}},PlotStyle->Blue,AxesLabel->{"f/Hz","Q"}]

HOMs/.zbaseFileRules//Last

It will be straightforward to put together functions that display this kind of information, starting from the ZbaseItem and ZbaseMode values.

Contents of the .beam file

exampleBeamFile=ToFileName[{ZbaseDataDirectory[],"lhc"},".beam"]

beam=Import[exampleBeamFile,"Table"]

Length[beam]

We can identify the first 13 of these parameters by comparison with

In Mathematica, we could store the information as a list of rules.  This has the advantage of being self-describing. Something like this:

Rule@@@Transpose@{{Machine,Particle,Kbunch,Nbunch,sigzInjection,sigzCollision,gammaInjection,gammaCollision,sige,slipFactor,frev,QsInjection,QsCollision,xx1,xx2},Flatten[beam]}

However this information is somewhat redundant as we can get much of it from the MAD description of the machine, see below.

Contents of the .loss file

lossFile = 
      ToFileName[
        Flatten[{ZbaseDataDirectory[], ZbaseMachine, ZbaseItem, 
              ZbaseMode} /. Options[
              Zbase]], ".loss"]

lossFileContents = Import[lossFile, "List"]

I don't understand this yet.

ABCI, URMEL, etc. files

Some scheme (rules would be best) for these.  The present system seems to depend on coding bunch lengths into the file name.

Other sources of impedance information

Besides the output from programs like GDFIDL, ABCI, URMEL etc., we may have impedances given by analytical formulas.  Make a scheme for including these - big advantage of Mathematica.

Generic formulas for typical impedances with rules giving specific values of parameters.

Probably good to use ConstantsUnits package to track units and dimensions of symbols appearing in the formulas.

A typical wake

Some of the data files are stored in gzipped archives (the .gz files).  This is not very convenient and somewhat system-dependent.  I think it would be better to change it.

Meanwhile, for demonstration purposes, I extracted a data file containing a wake by hand  (and renamed it)

sampleWakeFile=First[FileNames["*.dat"]]

sampleWakeFileContents=Import[First[FileNames["*.dat"]]]

The first 5 lines of the files are not part of the wake so drop them before plotting.

sampleWakeData=Drop[sampleWakeFileContents,5];

rawPlot=ListPlot[sampleWakeData,PlotStyle->Red,PlotJoined->True,AxesLabel->{"z/m","(Subscript[W, L])(z) / (V/pC/(m^2))"}]

How many points in this wake?

Length[sampleWakeData]

A moving average only takes us so far ...

?MovingAverage

DisplayTogether[rawPlot,ListPlot[MovingAverage[sampleWakeData,100],PlotJoined->True]
]

A moving median with repated smoothing is a bit better perhaps

DisplayTogether[rawPlot,ListPlot[MovingMedian[sampleWakeData,91,RepeatedSmoothing->True],PlotJoined->True]
]

DisplayTogether[rawPlot,ListPlot[ExponentialSmoothing[sampleWakeData,.02],PlotJoined->True]
]

Anyway the point is that we have access to many numerical techniques and can choose or implement what is appropriate for the physics.

How to transform this wake into a continuous function:

sampleWake=Interpolation[MovingAverage[sampleWakeData,100]]

Now we have a function that we can use it at any z.

sampleWake[.3]

Resampling the wake

How much memory does it use?

ByteCount[sampleWake]

This is not a big deal but we can show a standard trick to get a smaller but smooth representation of the wake using the adaptive sampling of Plot.  Restrict the range for now to avoid the noisy stuff at the end

resampledWakeData=Plot[sampleWake[z],{z,0,1.},PlotRange->All][[1,1,1,1]];

Length[resampledWakeData]

Now the interpolation is based on a set of points chosen so that it will be smooth enough.  These points are not evenly distributed:

ListPlot[resampledWakeData,PlotRange->All]

sampleWake1=Interpolation[resampledWakeData]

This is much smaller in memory:

ByteCount[sampleWake1]

Check accuracy

Plot[sampleWake[z]-sampleWake1[z],{z,0,1.},PlotStyle->Red,PlotRange->All]

This can be refined if it is not good enough.

Get rid of the noise

lwak = Interpolation[Drop[lwak1, 5]]

lwakma = Interpolation[MovingAverage[Drop[lwak1, 5], 50]]

Plot[((lwak[s] - lwakma[s])/lwak)[s], {s, 0.1, 1.99}]

Table[((lwak[s] - lwakma[s])/lwak)[s], {s, 0.1, 1.99, 
        0.05}]

lwak[ .8]

Take[lwak, 10]

Take[lwakma, 10]

Elements and how they occur in the LHC rings

The previous implementation contained information about the occurrences of the impedances in the ring.  In principle, each element with an impedance should appear in the MAD description of the LHC.    Therefore we should be able to use this as the basic source of information about the occurrences and not duplicate the information in the impedance database.   The Madtomma interface to the LHC has already been loaded so we have direct access to the information via its functions.

This is somewhat LHC-specific of course but similar things could be done for other machines.

Reminder about the Madtomma interface to the LHC optics

Since all Madtomma functions have usage messages, we have instant documentation to easily remind ourselves what they do.  To remind ourselves which functions exist, we only need to remember that their names begin with a few typical sets of characters.  To see the usage message of a function, just click on its name in the following displays.

?LHC*

?madLHC*

?mfs*

In the following, I shall typically display the usage message before using a function.

?madLHCOpticsVersion

SetOptions[madLHCOptics,madLHCOpticsVersion->"V6.500"];

?madLHCSetOptions

madLHCSetOptions["Injection"]

?LHCtwiss

?mfsMember

Example of the LHC RF cavities

Get the optics at the RF cavities:

?LHCtwiss

?mfsMember

RFcavityOptics=mfsMember[LHCtwiss["LHCB1"],"KEYWORD",{"RFCAVITY"}];

Show a simple table of interesting quantities

With[{cols={"NAME","S","L","BETX","BETY"}},TableForm[Transpose[mfsColumn[RFcavityOptics,cols]],TableHeadings->{{},cols}]]

These cavities belong to classes, the members of which should have the same impedances.   Madtomma has a function that can tell us what they are.

?LHCElementParent

So we can get the list of parent classes of each RF cavity.   The first call of this function builds a database and caches it so future calls will be much faster.

LHCElementParent/@mfsColumn[RFcavityOptics,"NAME"]

The list of different RF cavity types

Union[LHCElementParent/@mfsColumn[RFcavityOptics,"NAME"]]

Now we can count the number of occurrences of each type

?Frequencies

Frequencies[LHCElementParent/@mfsColumn[RFcavityOptics,"NAME"]]//TableForm

If we knew the impedance or wake for each type of cavity, it would now be easy to sum up the contributions using these occurrence counts as weights.

Logically now, one would expect the list of impedance element classe to be a superset of  the bottom level of ZbaseItems[].  This should be arranged, I think.

Optics information for impedance calculations

We can also easily do things like calculate the average beta at each element class.  Get a table of cavity names, lengths and (Subscript[β, y]) values, then select just the cavities belonging to a particular class.

With[{elclass="ACSCA"},eldata=Select[Transpose[mfsColumn[RFcavityOptics,{"NAME","L","BETY"}]],LHCElementParent[First[#]]==elclass&];
((Total[eldata/.{_,l_,b_}->l b])/(Total[eldata/.{_,l_,b_}->l]))
]

It's easy to generalise that

LHCElementClassAverage::usage="LHCElementClassAverage[seq,elclass,func] returns the length-weighted average of the optical function func (e.g. \"BETX\") over the elements of an LHC sequence belonging to the element class elclass.";

LHCElementClassAverage[seq_?LHCsequenceQ,elclass_String,func_String]:=Module[{eldata},
eldata=Select[Transpose[mfsColumn[LHCtwiss[seq],{"NAME","L",func}]],LHCElementParent[First[#]]==elclass&];
((Total[eldata/.{_,l_,b_}->l b])/(Total[eldata/.{_,l_,b_}->l]))
]

It would not be hard to implement expressions like ("BETX"-"BETY")/"DX" instead of simple strings for func using mfsToRules.

Check the previous result

LHCElementClassAverage["LHCB1","ACSCA","BETY"]

And get some more

LHCElementClassAverage["LHCB1","MQM","BETX"]

The caching of the element hierarchy in Mathematica makes these calculations quite fast enough.

LHCElementClassAverage["LHCB1","MB","BETX"]

LHCElementClassAverage["LHCB2","MB","BETX"]

Actually these are only the averages of the values at the exits of the elements.   Probably good enough for instability estimates.  But we could improve ...

With these techniques we could dispense with the need to duplicate optics and layout information in the impedance database.

Beams from the MAD description

The beam parameters are also in Madtomma

madLHCBeams/.Options[madLHCOptics]