Dispersion in twiss distribution & optics calculation requires verification

[Laurie Nevay]

Recent optics comparisons show dispersion off at the beginning. There was question over whether this was due to a fault in the generator. But it could also be from the optical function calculation. We need to verify dispersion ASAP.

[Laurie Nevay]

Confirmed. I do not get the correct spatial size (sigma_x) for a given value of dispersion.

 beam, particle="e-",
       energy=1.0*GeV,
       distrType="gausstwiss",
       betx=0.60000*m,
       bety=0.60000*m,
       alfx=-0.000001083824,
       alfy= 0.000001233215,
       emitx=25e-9,
       emity=25e-9,
       sigmaE=0.0005,
       sigmaT=1e-12;
       dispx=1*m;

I generate 100k primaries. There is at least an ambiguity here as sigma_x should be calculated from the fractional momentum spread and not energy. We document what we do correctly:

[http://www.pp.rhul.ac.uk/bdsim/manual/model_description.html#gausstwiss]

However, this will of course be quite wrong for low energy is the fractional energy spread is used. If the fractional momentum spread is 'put in' then the numbers should be as expected.

For the above example, I expect sigma_x = sqrt(0.625e-9 + ((5e-4)^2)1^2)) = 5.1478e-4m. With 100k primaries, I get (from running rebdsimOptics) sigma_x = 1.225e-4 ± 2.747e-7 ; disp_x = -8.743e-6 ±3.144e-4.

If I set the energy spread to 0 and dispx to 0, I get sigma_x = 1.227e-4 ± 2.752e-7m as expected.

[Laurie Nevay]

We should provide sigmaP as an option and be able to calculate between that and sigmaE in the bunch base so we're very clear about what's being used.

[Former user]

The set of input optics that I am using is:

beam,   X0=6.66178249721e-10*m,
        Xp0=-3.32352172158e-11, 
        Y0=1.6884459763e-12*m, 
        Yp0=3.08645629493e-14, 
        alfx=1.15489532801, 
        alfy=-1.73162448532, 
        betx=58.9580610872*m, 
        bety=106.422576724*m, 
        dispx=1.42505125208*m, 
        dispxp=-0.0175639234629, 
        dispy=0.149626369567*m, 
        dispyp=0.00594167483422, 
        distrType="gausstwiss", 
        emitx=3.60873863846e-10*m, 
        emity=3.60873863846e-10*m, 
        energy=6500.0*GeV, 
        particle="proton", 
        sigmaE=0.001;

which comes with its own issues, namely,

ERROR bunch generator sigma matrix is still not positive definite, giving up

Perhaps these (Laurie's and this one) are two manifestations of the same underlying problem, but I cannot even run BDSIM in the first place with this input distribution.

[Laurie Nevay]

I would suggest tackling basic spatial dispersion and P vs E first - at least you can run with this.

The non-positive definite matrix problem may be a consequence of incorrect usage, but that can be investigated afterwards.

[Laurie Nevay]

Indeed the typo in my original example fails to set dispersion:

   sigmaT=1e-12;
   dispx=1*m;

This results in a variable called 'dispx' being defined independent of the beam structure. This is valid syntax so we can't protect against it.

However, now I also hit the not positive definite problem with only dispx set.

  beam, particle="e-",
       energy=1.0*GeV,
       distrType="gausstwiss",
       betx=0.60000*m,
       bety=0.60000*m,
       alfx=-0.000001083824,
       alfy= 0.000001233215,
       emitx=25e-9,
       emity=25e-9,
      sigmaE=0.0005,
      sigmaT=1e-12,
      dispx=1*m;

I used seed=123

BDSRandom::SetSeed> Random number generator's seed = 123
BDSBunch::CreateMultiGauss> WARNING bunch generator sigma matrix is not positive definite

1.5e-08       2.70956e-14   0             0             0             2.5e-07       
2.70956e-14   4.16667e-08   0             0             0             0             
0             0             1.5e-08       -3.08304e-14  0             0             
0             0             -3.08304e-14  4.16667e-08   0             0             
0             0             0             0             1e-24         0             
2.5e-07       0             0             0             0             2.5e-07       

BDSBunch::CreateMultiGauss> adding a small error to zero diagonal elements
BDSBunch::CreateMultiGauss> WARNING bunch generator sigma matrix is still not positive definite

1.5e-08       2.70956e-14   0             0             0             2.5e-07       
2.70956e-14   4.16667e-08   0             0             0             0             
0             0             1.5e-08       -3.08304e-14  0             0             
0             0             -3.08304e-14  4.16667e-08   0             0             
0             0             0             0             1e-24         0             
2.5e-07       0             0             0             0             2.5e-07       

BDSBunch::CreateMultiGauss> adding a small error to all elements
BDSBunch::CreateMultiGauss> ERROR bunch generator sigma matrix is still not positive definite, giving up

1.5e-08       2.70956e-14   1e-50         1e-50         1e-50         2.5e-07       
2.70956e-14   4.16667e-08   1e-50         1e-50         1e-50         1e-50         
1e-50         1e-50         1.5e-08       -3.08304e-14  1e-50         1e-50         
1e-50         1e-50         -3.08304e-14  4.16667e-08   1e-50         1e-50         
1e-50         1e-50         1e-50         1e-50         1e-24         1e-50         
2.5e-07       1e-50         1e-50         1e-50         1e-50         2.5e-07       

[Former user]

I implemented a conversion based on gaussmatrix, but the error persists:

beam,   X0=6.66178249721e-10*m,
        Xp0=-3.32352172158e-11, 
        Y0=1.6884459763e-12*m, 
        Yp0=3.08645629493e-14, 
        distrType="gaussmatrix", 
        energy=6500.0*GeV, 
        particle="proton", 
        sigma11=2.12764233093e-08, 
        sigma12=-4.16771539351e-10, 
        sigma13=6.1625189346e-15, 
        sigma14=1.19476093059e-16, 
        sigma15=-7.98983365008e-10, 
        sigma16=0.0, 
        sigma22=1.42847534663e-11, 
        sigma23=-8.60489094556e-17, 
        sigma24=-2.75846258446e-18, 
        sigma25=2.2579947028e-11, 
        sigma26=0.0, 
        sigma33=3.84051264625e-08, 
        sigma34=6.24898018746e-10, 
        sigma35=-2.40895345912e-10, 
        sigma36=0.0, 
        sigma44=1.35588013217e-11, 
        sigma45=-4.34140643085e-12, 
        sigma46=0.0, 
        sigma55=3.9411427877e-11, 
        sigma56=0.0, 
        sigma66=0.0;

giving:

BDSBunch::CreateMultiGauss> ERROR bunch generator sigma matrix is still not positive definite, giving up

2.12764e-08   -4.16772e-10  6.16252e-15   1.19476e-16   -7.98983e-10  1e-50         
-4.16772e-10  1.42848e-11   -8.60489e-17  -2.75846e-18  2.25799e-11   1e-50         
6.16252e-15   -8.60489e-17  3.84051e-08   6.24898e-10   -2.40895e-10  1e-50         
1.19476e-16   -2.75846e-18  6.24898e-10   1.35588e-11   -4.34141e-12  1e-50         
-7.98983e-10  2.25799e-11   -2.40895e-10  -4.34141e-12  3.94114e-11   1e-50         
1e-50         1e-50         1e-50         1e-50         1e-50         1e-50

I really do not know what to make of this now, this is very peculiar to me.

[Laurie Nevay]

Just of note - this last matrix had no dispersion specified at all. Only coupling between x,xp,y,yp and t.

[Andrey Abramov]

Solved the problem by introducing all longitudinal-transverse and transverse-transverse couplings resulting from having an initial dispersion and a finite energy spread.

Beam is now generated successfully in cases where failure was previously observed and optical functions agree well.

[Former user]

Fixed in a349b80 (just for reference).

[Andrey Abramov]

Also, just for comparison, here is the beam matrix that actually works. Interestingly, it looks nothing like the MADX one Stuart posted earlier.