Correct Dipole tapering

[swhite2401]

This PR fixes issues with tapering method used until now that is making of polynomB[0].

[swhite2401]

@lfarv , I had a quick look and it seems the C part is ready, could you please confirm?

On the python side:

• update with the master
• adapt the interface to using the scaling attribute instead of PolynomB[0] for dipole
• benchmarking and tests
• Anything else?

I the c part is confirmed to be ok I can take care of adapting the python interface.

[simoneliuzzo]

I will be happy to try it with the latest FCC lattice

[lfarv]

This branch provides a scaling of magnet strength by changing the reference momentum. So it's exact and includes all effects. There are 2 possibilities:

• Acting on dipoles through a new scaling dipole attribute. scaling > 1 means increasing the dipole field. I just noticed that the modification is not done on the "RadPass" variant. I think it has to be done (I'll do it),
• Acting on a whole section of a lattice by inserting markers with a ChangePRefPass passmethod at both ends of the section. This element needs a single scaling attribute: use scaling=scale at entrance and scaling=1.0/scale at exit. There is no new element for that, but the Marker element is good enough.

Ok for you to modify the python interface to act on dipoles. The second option is more difficult to handle (it needs inserting elements in the lattice), but it avoids acting on single elements.

Do you think that a scaling attribute on multipoles also would make things easier? It's easy to implement!

[swhite2401]

Ok I implement the python interface, I don't think this is necessary for multipoles.
Can you let me know when this is done for RadPass methods? If I don't activate radiation in dipole there will not be be energy variations so tapering does not make much sense...

[lfarv]

@swhite2401: The correct scaling is now done in BndMPoleSymplectic4Pass and BndMPoleSymplectic4RadPass. For straight magnets, it must be done using PolynomB, while a scaling parameter could easily be added. For other bending magnet integrators (E2, …), it could be done also, if needed…

For info: I had once again to modify the test accuracy, in the Matlab tests this time, and for the Windows platform. I still do not understand what is the reason for those unexpected changes which already occurred several times.

[swhite2401]

This accuracy issue is really a mystery... it is true that I had to change these several times too...
Anyway, thanks for the update!

[swhite2401]

Tapering function has been adapted, results are identical for the tests I have done with the FCChh lattice. This is ready for review.

[simoneliuzzo]

I gave a quick look at this discussion and at the code. I admit that the name "Scaling" for a new Bending magnet attribute seems to me a bit too open to interpretation. May be it could be EnergyTaperingScaling? or TaperingScale, ... or whatever explains in a more immediate and concise way what is this parameter for.

[lfarv]

@simoneliuzzo : no, it represents a magnetic field scaling, you can use it for any purpose you like, not only for tapering !

[lfarv]

Practically, this Scaling parameter could also be used for steering. Then you get the full effect of changing the magnetic field, on focusing and higher field components AND on the intrinsic focusing effect due to the curvature.
The KickAngle parameter usually used for steering acts only on PolynomB[0]: it adds some pure dipole component and neglects the intrinsic focusing.

[swhite2401]

@lfarv, it is true that Scaling could need a bit more explanation. Maybe FieldScaling would be more self-explanatory?

[lfarv]

Ok for FieldScaling !

[lfarv]

I make the modifications…

[simoneliuzzo]

Also FieldScaling is not so clear to me, but I imagine that I may be did not understand correctly what this scaling is doing to the field. It changes Brho? Why not EnergyScaling?

Also,

• I see in the pull request files very little lines for documentation. For example several sentences in this discussion would be a nice contribution to documentation in the code
• is there a matlab AT tapering(ring) function ?

[lfarv]

@simoneliuzzo : the attribute changes the magnetic field of the dipole. For straight magnets, this is just a scaling of PolynomA/B. For curved magnets (sector dipole), the magnetic field is defined by both BendingAngle and PolynomA/B, so it's more tricky. For a sector dipole, the exact computation needs to temporarily change the reference momentum (not the momentum itself). For a rectangular dipole, this would be done differently. But in the end, there is:

• no change in energy
• no change in BRho

So FieldScaling is the exact definition of the effect on the dipole.

Then the tapering function scales the magnet strengths to accommodate the local energy.

@swhite2401 : in the doctoring of tapering, the 1st sentence:

"Scales magnet strengths with local energy to cancel the closed orbit and optics errors due to synchrotron radiations."

is very clear and describes exactly what happens. But the 2nd one:

"The dipole bending angle is changed by changing the reference momentum."

is may be responsible for the confusion reported by @simoneliuzzo. Is it really necessary to enter computational details?

@simoneliuzzo again: you are right, we need the Matlab equivalent. Any volunteer? otherwise, I'll do it, if there is not too much urgency.

[swhite2401]

"The dipole bending angle is changed by changing the reference momentum."

Ok to remove this part. If you could do the matlab that would be nice... did you change the attribute name in the python tapering function as well?

[lfarv]

Yes, change is done in the tapering function. I prepare the Matlab function, but I would like a test lattice…

[lfarv]

The Matlab function attapering is added.

@swhite2401 : from my experience with the Matlab function:

• the multipole scaling should be included if the iteration loop,
• iteration is useful only if multipoles are scaled

But this is based on results of a single lattice, so maybe it's not significant…

[swhite2401]

The Matlab function attapering is added.

@swhite2401 : from my experience with the Matlab function:

* the multipole scaling should be included if the iteration loop,

* iteration is useful only if multipoles are scaled

But this is based on results of a single lattice, so maybe it's not significant…

Yes it is well possible, should we integrate the multipole in the iteration loop? It is probably the right time... default is 1 anyway so it won't impact existing codes

[lfarv]

should we integrate the multipole in the iteration loop?

I would say so. As it is now, if you don't ask for multipoles, niter is useless, the convergence is immediate. On the other hand, if you ask for multipoles, they do change the orbit, and you have to explicitly call again tapering to compensate. Again, this is based on the few tests I made, but I think it makes sense. Ideally, we could even set the default for niter at 1 without multiples and 3 with multipoles.

[swhite2401]

Yes I agree! Please go ahead!
I would still keep the dipole in the iterative loop just in case. There may be examples for which it is needed such as combined function dipole.

[lfarv]

There is a problem with the multipoles: unlike for dipoles, the scalings accumulate on multipoles, so that they cannot be in the iteration loop. So for now, they are taken out, but set before the loop on dipoles, so that they are taken into account.

However it is a problem if the function is called twice on the same lattice: then it's wrong. The easiest solution would be to add the FieldScaling attribute on multipoles too.

[lfarv]

All pass methods now accept FieldScaling

[lfarv]

Here is a notebook checking that the tapering is now correct for dipoles, using either the FieldScaling attribute on all magnets, of scaling the whole ring with elements with the new ChangePRefPass passmethod. It also shows the present behaviour of the master branch (slightly wrong).

testscaling.pdf

[lfarv]

A full test of the FieldScaling attribute is added, the branch is ready for merging

[simoneliuzzo]

Dear @lfarv and @swhite2401

I am very happy that a test/example is added. I hope it will find a place in AT itself (as for collective effects).
However I looked at the example and I did not fully understand it.

To my understanding "tapering" is needed to adjust the magnets to the local energy of the beam, between two RF cavities.
I would have expected to see a test as

load lattice
disable 6D
get optics
enable 6D
get optics 6D (distorted)
enable tapering (OLD)
get optics 6D (back to the reference ones)
enable tapering (this PR)
get optics 6D (BETTER back to the reference ones)
plot H/V orbit, beta, dispersion vs reference for 6D, 6D + tapering, 6D + tapering new
plot gradient of magnets without tapering and with tapering (not periodic anymore).

(I attach a partial script I have for this, that needs modification)

Could this be added? It would be in my opinion a more intuitive example for "tapering users".

If FieldScaling has more potential than just correct optics due to large radiation effects, then examples of real use cases could be detailed for example via dedicated notebook tests (very appreciated) as the one proposed by @lfarv, or in the function help.

I think that mixing tapering (energy loss due to radiation, s dependent) and off-momentum (global offset of energy) may lead to non immediate understanding for the users.
When I set delta=value in get_optics, I would expect to see a change. Why should I cancel it with FieldScaling? What I would like to do instead is to look at on- and off- energy optics with tapering.

I also have an operational comment.
Let say that we have a lattice with radiation and tapering set with "FieldScaling". Where/how does a user extracts the gradients of the magnets to be converted to PS currents?
For example, we could think to set "tapering" corrections in the EBS-SR lattice to test if they have a positive impact on beam properties. The individual magnet gradients would be Polynom(A/B)*FieldScaling? If this is the case, it should be reported in the help of AT at the level of each element creator that allows the FieldScaling attribute and in the lattice itself with a sentence such "operational individual magnets gradients are Polynom(A/B)*FieldScaling"
If this is the case, I would consider this also a big conceptual change for AT, worth at least a minor release (0.5.0).

thank you for your help

best regards
Simone

import at
import matplotlib.pyplot as plt

# load lattice
lattice_file = '/machfs/liuzzo/FCC/v61_t/Errors/FCCee_hfd61_ttbar_diag_cor.mat'
ring_var_name = 'r'
# load starting, ideal lattice
r0 = at.load_lattice(lattice_file, mat_key=ring_var_name)
ind_bpms = range(len(r0))

# compute optics
r0.disable_6d()
_, _, opt0 = r0.linopt6(ind_bpms)

# compute optics 6D
r0.enable_6d()
_, _, optrad = r0.linopt6(ind_bpms)

# apply tapering
r0.tapering()

# compute optics 6D and tapering
_, _, optradtap = r0.linopt6(ind_bpms)

s = at.get_s_pos(r0, ind_bpms)
# plot beta-beating dispersion deviation and Delta COD
fig, ax = plt.subplots(4, 2)
fig.subplots_adjust(hspace=0.5, wspace=0.5)
ax[0, 0].plot(s, [h[0]*1e6 for h in optrad.closed_orbit], label='6D')
ax[0, 0].plot(s, [h[0]*1e6 for h in optradtap.closed_orbit], label='6D + tapering')
# ax[0, 0].plot(s, [h[0]*1e6 for h in opt0.closed_orbit], label='4D')
ax[0, 0].set_ylabel('hor. closed orbit [$\mu$m]')
ax[0, 0].legend()
ax[0, 1].plot(s, [h[2]*1e6 for h in optrad.closed_orbit], label='6D')
ax[0, 1].plot(s, [h[2]*1e6 for h in optradtap.closed_orbit], label='6D + tapering')
# ax[0, 1].plot(s, [h[2]*1e6 for h in opt0.closed_orbit], label='4D')
ax[0, 1].set_ylabel('ver. closed orbit [$\mu$m]')

ax[1, 0].plot(s, [(h[0]-h0[0])*1e3 for h, h0 in zip(optrad.dispersion, opt0.dispersion)], label='6D')
ax[1, 0].plot(s, [(h[0]-h0[0])*1e3 for h, h0 in zip(optradtap.dispersion, opt0.dispersion)], label='6D + tapering')
ax[1, 0].legend()
ax[1, 0].set_ylabel('hor. $\Delta\eta$ [mm]')
ax[1, 1].plot(s, [(h[2]-h0[2])*1e3 for h, h0 in zip(optrad.dispersion, opt0.dispersion)], label='6D')
ax[1, 1].plot(s, [(h[2]-h0[2])*1e3 for h, h0 in zip(optradtap.dispersion, opt0.dispersion)], label='6D + tapering')
ax[1, 1].set_ylabel('ver. $\Delta\eta$ [mm]')


# ax[2, 0].plot([h[0] for h in opt0.beta], label='no rad')
ax[2, 0].plot(s, [(h[0]-h0[0])/h0[0] *100 for h, h0 in zip(optrad.beta, opt0.beta)], label='6D')
ax[2, 0].plot(s, [(h[0]-h0[0])/h0[0] *100 for h, h0 in zip(optradtap.beta, opt0.beta)], label='6D + tapering')
ax[2, 0].set_ylabel('hor. $\Delta\\beta / \\beta_0$ [%]')
ax[2, 0].legend()
ax[2, 0].set_xlabel('s [m]')
# ax[2, 1].plot([h[1] for h in opt0.beta], label='no rad')
ax[2, 1].plot(s, [(h[1]-h0[1])/h0[1] *100 for h, h0 in zip(optrad.beta, opt0.beta)], label='6D')
ax[2, 1].plot(s, [(h[1]-h0[1])/h0[1] *100 for h, h0 in zip(optradtap.beta, opt0.beta)], label='6D + tapering')
ax[2, 1].set_ylabel('ver. $\Delta\\beta / \\beta_0$')
ax[2, 1].set_xlabel('s [m]')

ax[3, 1].plot(s, [h0[5] for h0 in optrad.closed_orbit], label='6D')
ax[3, 1].plot(s, [h[5] for h in optradtap.closed_orbit], label='6D + tapering')
ax[3, 1].set_ylabel('ct')
ax[3, 1].legend()
ax[3, 1].set_xlabel('s [m]')

ax[3, 0].plot(s, [h0[4]*100 for h0 in optrad.closed_orbit], label='6D')
ax[3, 0].plot(s, [h[4]*100 for h in optradtap.closed_orbit], label='6D + tapering')
ax[3, 0].set_ylabel('delta [%]')
ax[3, 0].legend()
ax[3, 0].set_xlabel('s [m]')

# plt.savefig('tapering_effect.png')

plt.show()

[simoneliuzzo]

Dear @lfarv and @swhite2401,

I run 2 cases, only in python for now. I removed all ExactPass methods

1. FCC lattice with the (exact_multipole_pass) AT
   

2. FCC lattice with this branch (dipole_tapering)
   

I see a few differences. Most remarkably in this branch (dipole_tapering) the vertical dispersion with radiation seems slowly diverging.

[simoneliuzzo]

Dear @swhite2401 and @lfarv
for the same lattice as above I tried the function attapering in matlab AT.
I expected to get identical results to those in python.
Here below the output.

[lfarv]

@simoneliuzzo : concerning tapering, this all looks good. Now for the differences between Matlab and python, they mostly appear without applying tapering (horizontal beta and eta), so it does not look related to this branch. It could interesting to look at that independently of this PR…

In the vertical plane, closed orbit and dispersion seems to be in the numerical noise. Are there any elements in your lattice generating vertical orbit or dispersion?

[simoneliuzzo]

Dear @lfarv and @swhite2401 ,

here the two figures updated as requested.

I see differences that are not related to this pull request so I do not feel confortable approving it. If for you those differences are irrelevant please approve, I will not oppose.

[lfarv]

@simoneliuzzo : Where do you see problematic differences? For me the tapering effect looks correct. The difference with the master branch is the correction of the dipole focusing effect when applying tapering. This focusing effect in dipoles appears for small bending radii, so I'm not surprised that with your lattice, it does not make a significant difference…
Again, in the vertical plane, is there any reason why the orbit and dispersion are not zero ?

[swhite2401]

Hello both, I am also more concerned about this lattice giving strange v orbit and dispersion... did you get a chance to check with MADX how it looks like? do you see the same with the exact passmethods? It would be nice to understand what the problem is.

If you have the madx lattice available I can take a look.

Tapering looks fine to me.

[simoneliuzzo]

Dear @swhite2401 and @lfarv,

to my understanding tapering is supposed to have no effect in the vertical plane. So I expected no change in the vertical plane, while I see some (very small) changes. This is what (lightly) worries me. The modification done in this branch seem to have an effect on parameters that were supposed to be unaffected.

We can check with MADX (M.Hofer at CERN has the files, but I bet there are small differences here and there, and it will take ages to get to exactly corresponding lattices and tracking). I personally do not see the interest here (may be for FCC studies) as it is unrelated to this pull request.

For me the problem is simply: vertical should have not changed at all. Why does it change?

Not to mention that I do not clearly see the added value of such large code changes (21 files changed! all passmethods!). Seen the changes in final optics provided by this pull request compared to the the previous "wrong" tapering function it does not seem really motivated. Let's also remember that FCC is unlikely to be realized and this effect has applications only for this lattice.

Concerning the lattice:

• The file is this: /machfs/liuzzo/FCC/v61_t/Errors/FCCee_hfd61_ttbar_diag_cor.mat'
• It works only in the exact_multipole_pass branch, unless you switch back all passmethods to the standard ones (as I do to make the figures above).
• There are no skew quad
• There are no vertical dipoles

For completeness I run the same figures for EBS (no vertical dipoles, no skew quadrupoles)

in this case tapering reduces a bit the computed (extremely small) vertical dispersion. Also in this case, the computed vertical dispersion changes in dipole_tapering branch. By the way, this appears to me as the most visible change, thanks to the scales.

[simoneliuzzo]

Dear @lfarv and @swhite2401,

@swhite2401 convinced me that small vertical dispersion changes could appear.

I approve for the sake of moving forward to more relevant topics such as ExactPass.

However I strongly encourage NOT to merge, as it is an overkill of code changes for no net final effect.

Apparently the approximations done in the first version of the tapering functions (just few lines of easy to read code) where rather good approximations!

[lfarv]

@simoneliuzzo, @swhite2401 : If I understand, both lattices should give exactly zero vertical closed orbit and vertical dispersion. So what we get is just numerical noise (and its range is acceptable). You can also notice that it gets larger with the number of elements in the lattice. I do not see anything suspicious in these results. The problem seems to be in the numpy.linalg.solve function: when the solution of the linear system is an exact zero, its solution is a very small but non-zero vector. When propagated along a long lattice, you may get what appears here. In my experience, in Matlab I instead always got an exact zero (but that does not seem to be true in @simoneliuzzo's example).

Concerning the number of modifications, it's obviously due to the FieldScaling parameter necessary for the tapering: all "magnet" passmethods must get it. I recall that in the master branch, calling the "tapering" function twice will cumulate the field modifications (see #623 (comment)). Also there is no way to check if tapering was applied, and no easy way back. All this is now avoided with the FieldScaling parameter.

So if there is no other problem, I will indeed merge this branch to allow further work on C integrators.

[swhite2401]

Somehow I was the author of this branch but in reality I am not so please go ahead

[lfarv]

@swhite2401 : why ? Any change since #623 (comment)?

If the tapering is not useful, we have at least to remove the wrong one from the master branch. But I will anyway keep all the cleaning done in the the C integrators.

[swhite2401]

Well you obviously took over! All I am saying is tapering is useful and used and this branch is fine for me, so please go ahead an merge it as is if you are happy with it

[swhite2401]

I am not so please go ahead

Ah this is was missleading it ment -> I am not author, please go ahead with the merge

[lfarv]

@swhite2401 : complete misunderstanding ! Sorry, I understood that you were "not pleased" with the way the PR turned out…