LEP with synchrotron radiation in MAD8

John Jowett

DateString[]

Started this to show how the chromaticity was distributed between low-beta insertions and arcs in LEP, in response to a question from Ferdinand Willeke (around 2006).  
Continued as a revival of radiation calculations in MAD8 asked for by Massimo Giovannozzi et al, with a view to recovering the same functionality in MADX.

Setup

SetDirectory[NotebookDirectory[]]

SystemOpen["C:\\Users\\jowett\\OneDrive\\CERN - OneDrive\\LEPpostfacto\\Radiation"]

MAD program

Needs["Madtomma`RunningMAD`"]

mad8prog=Last[installedMADprograms[]]
mad8dir=DirectoryName[mad8prog]
mad8dic=First[FileNames["*.dic",mad8dir]]

Check the DICT environment variable (had to be set before this Mathematica session, in Windows System settings)

Environment["DICT"]==mad8dic

SetOptions[runMAD,MADprogram->mad8prog]

Some mfs packages  that should work with MAD8

Needs["Madtomma`Mfs`OpticsUtilities`"]

Needs["BeamLineGraphics`"]

Need to use MAD8 for old LEP input files.

Hoping this is not needed later:

SetOptions[MADcommand, MADXsyntax -> False]

SetOptions[MADcommand, MADXsyntax -> False]

We would like to be able to just do this:

Needs["Graphics`ColorListPlot`"]

Auxiliary functions needed

This should probably be moved to the Mfs package.  Note that it can be used to replace the file.

MAD8tfsFileFix::usage="MAD8tfsFileFix[fnold,fnnew] creates a new version of a TFS file created by MAD8 with the lines in the right order (as MADX).";

MAD8tfsFileFix[fnold_String,fnnew_String]:=Module[{f,fchar1,ib1,ib2,f1},
f=Import[fnold,"Lines"];
fchar1=StringTake[#1,1]&/@f;
ib1=Position[fchar1,"@"][[1,1]];
ib2=Position[fchar1,"@"][[-1,1]];
f1=Join[f[[ib1;;ib2]],f[[{1,ib1-1}]],f[[ib2+1;;]]];
Export[fnnew,f1,"Lines"]
]

LEP optics database

LEPopticsDatabase=ToFileName[{If[$OperatingSystem=="Windows","\\\\cernbox-smb.cern.ch","\\"],
"eos","project","l","lhcions","JMJ","slath","machines"}]

However we are going to run into the old problem with MAD8’s line length of 80 characters.

StringLength[LEPopticsDatabase]

So we map the following folder to network drive Z: in Windows explorer.

"\\\\cernbox-smb.cern.ch\\eos\\project\\l\\lhcions\\JMJ\\slath\\machines"

SystemOpen["\\\\cernbox-smb.cern.ch\\eos\\project\\l\\lhcions\\JMJ\\slath\\machines"]

LEPopticsDatabase="Z:"

FileNames["*",LEPopticsDatabase]

Compute a LEP optics, basic, without radiation

A typical LEP optics from around 1998 http://cds.cern.ch/record/702448/files/sl-note-98-020.pdf  with arc cell phase advances (Subscript[μ, x])=102°,  (Subscript[μ, y])=90°.

It seems that we have to change the backslashes to slashes in filenames.  I don’t remember that.  I suppose the old MADXsyntax option being False dealt with that.

seqFile=StringReplace[FileNameJoin[{LEPopticsDatabase,"lep98","lep98_90.seq"}] ,"\\"->"/"]

FileType[seqFile]

strengthFile=StringReplace[FileNameJoin[{LEPopticsDatabase,"lep98","FBv1","F05150B98v1.lep"}] ,"\\"->"/"]

FileType[strengthFile]

MADcall[strengthFile]//MADPrint

Construct some commands manually since MadInput has been tuned to MAD-X lately.  Remember how to break long lines. 
Fix the output TFS file straight away.

Options[runMAD]

Construct the MAD8 run.

Export["LEPoptics.mou",
runMAD[{MADcall[seqFile],
MADcall[strengthFile],
MADuse["LEP"],
"SELECT,OPTICS,FULL",
"OPTICS,CENTRE,FILENAME=\"LEPoptics.tfs\",COLUMNS=NAME,KEYWORD,S,L,K0L,K1L, K2L,&","VOLT,ALFX, ALFY, BETX, BETY, DPX, DX,DY, MUX, MUY, PX, PY, X, Y" ,
},runMADTemporaryInputFile->"LEPoptics.madx"],
"Lines"]
LEPoptics=tfsRead[MAD8tfsFileFix["LEPoptics.tfs","LEPoptics.tfs"]]

Look back at contents of the MAD file

FilePrint["LEPoptics.madx"]

Basic checks on the optics

mfsKeyValue[LEPoptics]//TableForm

mfsColumnNames[LEPoptics]

mfsColumn[LEPoptics,{"S","BETX"}]//Transpose//ListPlot

mfsColumn[LEPoptics,{"S","BETY"}]//Transpose//ListPlot

All the types of elements in the machine

Union[mfsColumn[LEPoptics,"KEYWORD"]]

Chromaticity contributions

Functions that will take an optics and return each elements contribution to the chromaticity (using simple perturbation formulas).  The horizontal is not quite right since I forgot to include the bend strength so it does not have the weak forcusing due to the bends.

dQxp[opt_mfs]:=mfsColumn[opt,{"S","K1L","BETX","K2L","DX"}]/.{s_,k1l_,bx_,k2l_,dx_}->-((k1l bx)/(4 π))+k2l dx (bx/(4 π));

The natural part, ignoring sextupoles

dQxpNatural[opt_mfs]:=mfsColumn[opt,{"S","K1L","BETX","K2L","DX"}]/.{s_,k1l_,bx_,k2l_,dx_}->-((k1l bx)/(4 π));

dQyp[opt_mfs]:=mfsColumn[opt,{"S","K1L","BETX","BETY","K2L","DX"}]/.{s_,k1l_,bx_,by_,k2l_,dx_}->((k1l by)/(4 π))-k2l dx (by/(4 π));

dQypNatural[opt_mfs]:=mfsColumn[opt,{"S","K1L","BETX","BETY","K2L","DX"}]/.{s_,k1l_,bx_,by_,k2l_,dx_}->((k1l by)/(4 π));

ListPlot[dQxp[LEPoptics],PlotRange->All]

ListPlot[dQyp[LEPoptics],PlotRange->All]

Sum up to get the total chromaticity of the ring

Plus@@dQxp[LEPoptics]

Plus@@dQyp[LEPoptics]

mfsKeyNames[LEPoptics]

Not too different from the MAD values

mfsKeyValue[LEPoptics,"XIX"]

mfsKeyValue[LEPoptics,"XIY"]

mfsKeyValue[LEPoptics,"QX"]

mfsKeyValue[LEPoptics,"QY"]

Show how the chromaticity evolves around the ring, big contributions from insertions, over-compensated arcs

ListPlot[FoldList[Plus,0,dQxp[LEPoptics]],Joined->True]

ListPlot[FoldList[Plus,0,dQyp[LEPoptics]],Joined->True]

Insertions and their chromaticity

mfsColumn[mfsRemoveDrifts[LEPoptics],"NAME"]

Get a list of elements delimiting IRs from the other parts

limitels=mfsColumn[mfsSelect[LEPoptics,StringMatchQ[mfsColumnValue[LEPoptics,#,"NAME"],"PU.QS11*"]||StringMatchQ[mfsColumnValue[LEPoptics,#,"NAME"],"PU.QL11*"]||StringMatchQ[mfsColumnValue[LEPoptics,#,"NAME"],"IP*"]&],"NAME"]

LEPsections= Transpose[{limitels,RotateLeft[limitels,1]}]

Optics plots of the sections

SetOptions[BeamLineElementGraphic,{BeamLineGraphicWidth->40,BeamLineGraphicDot->0.3`,BeamLineGraphicLabelLength->2`}]

SetOptions[mfsOpticsClassicBeamLinePlot,DyFactor->10]

mfsOpticsClassicBeamLinePlot[mfsOpticsRange[LEPoptics,LEPsections[[4]]]]

mfsOpticsClassicBeamLinePlot[mfsOpticsRange[LEPoptics,LEPsections[[5]]]]

mfsOpticsClassicBeamLinePlot[mfsOpticsRange[LEPoptics,#],PlotLabel->#]&/@LEPsections

Chromaticity of sections

Plot the natural chromaticity contributions of each piece of the ring

ListPlot[dQxpNatural[mfsOpticsRange[LEPoptics,#]],Joined->True,PlotRange->All,PlotLabel->#]&/@LEPsections

Sum up the net contributions to the natural chromaticity from each piece of the ring

{#,(Plus@@dQxpNatural[mfsOpticsRange[LEPoptics,#]])}&/@LEPsections

Show same as a table for clarity

TableForm[%]

In the vertical plane

{#,(Plus@@dQypNatural[mfsOpticsRange[LEPoptics,#]])}&/@LEPsections

TableForm[%]

This shows that the contributions from the odd-numbered straights are small while the contributions from the experimental low-beta insertions in the 4 even-numbered IPS are typically (Subscript[ΔQ, y])'≈-20 compared with a contribution from each of the 8 arce of (Subscript[ΔQ, y])'≈-12.

The contributions of each piece of the ring, including the correction from the sextupoles

{#,(Plus@@dQxp[mfsOpticsRange[LEPoptics,#]])}&/@LEPsections//TableForm

{#,(Plus@@dQyp[mfsOpticsRange[LEPoptics,#]])}&/@LEPsections//TableForm

Optics with RF system and radiation

Setup and compute optics and beam parameters using 6D formalism

Build the 6D orbits and optics with radiation, storing parts of the input expression in functions for easy and consistent re-use later.  The initial value of the LAG parameter is necessary so that RF cavities are close enough to the stable phase and the closed-orbit finder can find a solution.  The LAG is then matched so that the Taylor maps for tracking are computed by expansion around the 6D closed orbit.

This section was first developed with zero betatron coupling. To introduce a moderate level of betatron coupling I inserted

MADset["KQT1.L2",0.001]

Feed that matched value of the cavity LAG back into the definition so we don’t need to keep re-matching it.

MAD8LEPoptics[]:=MADtag["LEPoptics",{MADcall[seqFile],
MADcall[strengthFile],
MADset["LAG",0.3609269],
MADset["VRFC",2.2],MADset["VRFSC",42.],
MADset["LRFC","LAG"],MADset["LRFSC","LAG"],
MADset["KQT1.L2",0.001],
MADuse["LEP"]}]

MAD8LEPbeam1[]:=MADtag["LEPBeam1",{"BEAM,PARTICLE=POSITRON,ENERGY=94.,RADIATE=.TRUE."}]

MAD8OPTICS[fn_String]:=MADtag["OPTICStable",{"SELECT,OPTICS,FULL",
"OPTICS,FILENAME=\""<>fn<>"\",&",
"COLUMNS=NAME,KEYWORD,S,L,K0L,K1L, K2L,VOLT,&",
"ALFX, ALFY, BETX, BETY, DPX, DX,DY, MUX, MUY, PX, PY, X, Y"}]

MAD8EMIT[fn_String]:=MADtag["EMIT",
{"ASSIGN,PRINT=\""<>fn<>"\"",
"EMIT",
"ASSIGN,PRINT=TERM"
}]

MAD8TWISS3[fn_String]:=MADtag["TWISS3",{
"SELECT,TWISS3,FULL",
"TWISS3,SAVE=TWISS3",
"ARCHIVE,TABLE=TWISS3,FILENAME=\""<>fn<>"\""
}]

Export["LEPRFoptics.mou",
runMAD[{"ASSIGN,PRINT=TERM",
MAD8LEPoptics[],
MADset["LAG",0.4],
MAD8LEPbeam1[],
MADtag["MatchLAG",
{"CELL",
"CONSTRAINT,RANGE=IP1,T=0.0",
"VARY,LAG,STEP=1.E-4",
"LMDIF,CALLS=100,TOLERANCE=1.E-6",
"ENDMATCH"}],
MAD8OPTICS["LEPRFoptics.tfs"] ,
MAD8EMIT["LEPRFemit.mou"],
MAD8TWISS3["LEPRFTwiss3.tfs"]
},runMADTemporaryInputFile->"LEPRFoptics.madx"],
"Lines"]

LEPRFoptics=tfsRead[MAD8tfsFileFix["LEPRFoptics.tfs","LEPRFoptics.tfs"]]

LEPRFemit=MADEMITRead["LEPRFemit.mou"]

LEPRFTwiss3=tfsRead[MAD8tfsFileFix["LEPRFTwiss3.tfs","LEPRFTwiss3.tfs"]]

SystemOpen["LEPRFoptics.mou"]

Orbit in 6D, energy sawtooth

mfsRow[LEPRFoptics,{"NAME","VOLT"}]

ListPlot[Transpose[{mfsColumn[LEPRFoptics,"S"],Accumulate[mfsColumn[LEPRFoptics,"VOLT"]]}]]

Total[mfsColumn[LEPRFoptics,"VOLT"]]

That agrees with the old SL Note.

The TWISS3 optics contains the functions

mfsKeyValue[LEPRFTwiss3]//TableForm

mfsColumnNames[LEPRFTwiss3]

Select[mfsColumnNames[LEPRFTwiss3],StringMatchQ[#,"*CO"]&]

Plots of the closed-orbit components, including the energy sawtooth (Subscript[p, tc]):

ListPlot[mfsRow[LEPRFTwiss3,{"S",#}],AxesLabel->{"s/m",#},PlotRange->All]&/@Select[mfsColumnNames[LEPRFTwiss3],StringMatchQ[#,"*CO"]&]

Generalised dispersion functions

Check that the s values in the OPTICS and TWISS3  files are the same.

ListPlot[mfsRow[LEPRFTwiss3,"S"]-mfsRow[LEPRFoptics,"S"],PlotRange->All]

To get the “horizontal dispersion” amplitude function (ie the coupling of mode 3 into mode 1), we need the function (Subscript[γ, t3]) which is close to being a constant around the ring.

ListPlot[mfsRow[LEPRFTwiss3,{"S","GAMT3"}]]

Then we can reconstruct the dispersion, to be understood as an amplitude function (Sqrt[((Subscript[β, x3])/(Subscript[γ, t3]))]), not the shift of closed orbit with momentum

ListPlot[mfsRow[LEPRFTwiss3,{"S","BETX3","GAMT3"}]/.{s_,bx3_,gt3_}:>{s,(Sqrt[(bx3/gt3)])}]

Comparing with the dispersion function in the basic “TWISS” calculation

ListPlot[mfsRow[LEPRFoptics,{"S","DX"}],PlotRange->All]

With[{s=mfsRow[LEPRFTwiss3,"S"],
Dx=mfsRow[LEPRFoptics,"DX"],DxMR=mfsRow[LEPRFTwiss3,{"BETX3","GAMT3"}]/. {bx3_,gt3_}:>(Sqrt[(bx3/gt3)])},
{ ListPlot[Transpose[{Dx,DxMR}],AspectRatio->1],
ListPlot[Transpose[{s,Dx-DxMR}],PlotRange->All]}]

The vertical dispersion

ListPlot[mfsRow[LEPRFTwiss3,{"S","BETY3","GAMT3"}]/.{s_,by3_,gt3_}:>{s,(Sqrt[(by3/gt3)])},PlotRange->All,Joined->True]

Bassetti component of the orbit

The Bassetti component of the horizontal closed orbit is (Subscript[x, c])-(Subscript[D, x])(Subscript[p, tc])

With[{s=mfsRow[LEPRFTwiss3,"S"],DxMR=mfsRow[LEPRFTwiss3,{"BETX3","GAMT3"}]/. {bx3_,gt3_}:>(Sqrt[(bx3/gt3)]),xc=mfsRow[LEPRFTwiss3,"XCO"],ptc=mfsRow[LEPRFTwiss3,"PTCO"]},
GraphicsGrid[Transpose[{{
ListPlot[Transpose[{s,xc}],AxesLabel->{"s/m","(Subscript[x, c])/m"},Joined->True],
ListPlot[Transpose[{s,DxMR ptc}],AxesLabel->{"s/m","(Subscript[D, x])(Subscript[p, tc])/m"},Joined->True],
ListPlot[Transpose[{s,xc-DxMR ptc}],AxesLabel->{"s/m","((Subscript[x, c])-(Subscript[D, x])(Subscript[p, tc]))/m"},Joined->True]}}],
ImageSize->600]]

Generalised β-functions

Names of the β-functions

beta3=Select[mfsColumnNames[LEPRFTwiss3],StringMatchQ[#,"BET*"]&]

ListPlot[mfsRow[LEPRFTwiss3,{"S",#}],AxesLabel->{"s/m",#},PlotRange->All]&/@beta3

mfsRow[LEPRFoptics,"BETX"]

mfsRow[LEPRFTwiss3,"BETX1"]

Small effect of radiation on horizontal optics but not on vertical

With[
{s=mfsRow[LEPRFTwiss3,"S"],bx=mfsRow[LEPRFoptics,"BETX"],bx1=mfsRow[LEPRFTwiss3,"BETX1"]},
{ListPlot[bx-bx1],ListPlot[Transpose[{bx,bx1}]]}
]

With[
{s=mfsRow[LEPRFTwiss3,"S"],by=mfsRow[LEPRFoptics,"BETY"],by2=mfsRow[LEPRFTwiss3,"BETY2"]},
{ListPlot[by-by2],ListPlot[Transpose[{by,by2}]]}
]

Beam parameters

Use  functions for working with Associations.

Keys[LEPRFemit]

Values[LEPRFemit]

LEPRFemit["Emittance"]

LEPRFemit["DampingPartition"]

LEPRFemit["DampingTime"]

Tracking with radiation damping

Seems we can’t use MADrun since it doesn’t support the FPRINT option.  Fix later maybe.     LEP tracking in 1990s was done using the binary file output.

Options[MADrun]

Export["LEPRFoptics.mou",
runMAD[{"ASSIGN,PRINT=TERM",
MAD8LEPoptics[],
MAD8LEPbeam1[],
MAD8EMIT["TrackEMIT.mou"],
"ASSIGN,PRINT=\""<>"track.mou"<>"\"",
MADtrack[{Table[MADstart[0.005 j,0,0.01,0,0,0],{j,1,10}],
 
"RUN,TURNS=200,FPRINT=1"
},Damp->True],
"ARCHIVE,TABLE=TRACK,FILE=\"track.tfs\"",
"ASSIGN,PRINT=TERM"
},runMADTemporaryInputFile->"LEPRFdtrack.madx"],
"Lines"]

SystemOpen["LEPRFoptics.mou"]

Timing[tr=Import["track.mou","Table"];]

tr[[1]]

tr[[2]]

tr[[3]]

Length[tr]

tr[[9;;19]]//TableForm

tr[[14]]

Length[tr[[9;;]]]

tr[[-10;;-1]]//TableForm

Take[Cases[tr,{_,"x",__}|{"y",__}|{"t",__}],20]//TableForm

Take[Cases[tr,{n_Integer,"x",q_,p_,__}|{"y",q_,p_,__}|{"t",q_,p_,__}->{q,p}],20]//TableForm

Take[Partition[Cases[tr,{_Integer,"x",q_,p_,__}|{"y",q_,p_,__}|{"t",q_,p_,__}],3],3]//Timing

tr1=Partition[Cases[tr,{_Integer,"x",q_,p_,__}|{"y",q_,p_,__}|{"t",q_,p_,__}],3]

tr2=Cases[tr1,{{n_Integer,"x",x_,px_,__},{"y",y_,py_,__},{"t",t_,pt_,__}}:>{n,x,px,y,py,t,pt}]

tr21=Cases[tr2,{2,qp__}->{qp}];

ColorListPlot[tr21[[All,1;;2]],PlotRange->All]

FindPeaks[Re[Fourier[tr21[[All,1]]]]]/. {f_,a_}->{N[((f-1)/200)],a}

MapIndexed[{First[#2],#1}&,Re[Fourier[tr21[[All,1]]]]]/. {f_,a_}->{N[((f-1)/200)],a}

ListPlot[%,Joined->True,PlotRange->{0,All}]

ColorListPlot[tr21[[All,1]],PlotRange->All,Joined->True]

ColorListPlot[tr21[[All,3;;4]],PlotRange->All]

ColorListPlot[tr21[[All,3]],PlotRange->All]

ColorListPlot[tr21[[All,5;;6]],PlotRange->All,Joined->True]

ColorListPlot[tr21[[All,5]],PlotRange->All,Joined->True]

Tracking with quantum fluctuations

Seems we can’t use MADrun since it doesn’t support the FPRINT option.  Fix later maybe.     LEP tracking in 1990s was done using the binary file output.

Options[MADtrack]

{OnePass->False,Damp->False,Quantum->True,Deltap->0}

MAD8LEPbeam1[]//MADcommand//InputForm

Export["LEPRFoptics.mou",
runMAD[{"ASSIGN,PRINT=TERM",
MAD8LEPoptics[],MADset["VRFSC",50.],
MAD8LEPbeam1[],
MAD8EMIT["TrackEMIT.mou"],
"ASSIGN,PRINT=\""<>"qtrack.mou"<>"\"",
MADtrack[{MADstart[0,0,0,0,0,0], 
"RUN,TURNS=1000,TABLE=TRACK,FPRINT=1"
},Damp->True,Quantum->True],
"ARCHIVE,TABLE=TRACK,FILE=\"qtrack.tfs\"",
"ASSIGN,PRINT=TERM"
},runMADTemporaryInputFile->"LEPqtrack.madx"],
"Lines"]

SystemOpen["LEPRFoptics.mou"]

tr=Import["qtrack.mou","Table"];
tr1=Partition[Cases[tr,{_Integer,"x",q_,p_,__}|{"y",q_,p_,__}|{"t",q_,p_,__}],3] ;
tr2=Cases[tr1,{{n_Integer,"x",x_,px_,__},{"y",y_,py_,__},{"t",t_,pt_,__}}:>{n,x,px,y,py,t,pt}];
tr21=Cases[tr2,{1,qp__}->{qp}];

ColorListPlot[tr21[[All,1;;2]],PlotRange->All]

ColorListPlot[tr21[[All,1]],PlotRange->All,Joined->True]

ColorListPlot[tr21[[All,3;;4]],PlotRange->All]

ColorListPlot[tr21[[All,3]],PlotRange->All]

ColorListPlot[tr21[[All,5;;6]],PlotRange->All,Joined->False]

ColorListPlot[tr21[[All,5]],PlotRange->All,Joined->True]

It seems there is a bug in the quantum tracking output so you can only get the first 33 turns on the terminal output and nothing in the archived TRACK table.  In the 1990s, I had to use the binary output files in EPIO format with PAW.