Infinite loop with certain twiss parameters

[Laurie Nevay]

The following Twiss parameters result in an infinite loop in CLHEP due to the checks for positive-definiteness on the 6x6 Gauss Twiss matrix.

It seems to be in the diagonalize function inside CLHEP which uses a while loop and is unbound. We should avoid this as we can’t change it but ensure the check, which ultimately ensures that we can invert the matrix as required from generating from it.

‌

beam,   X0=0.0*m,
    Xp0=0.0, 
    Y0=0.0*m, 
    Yp0=0.0,
    alfx=-2.26218442917668e-07, 
    alfy=-1.888047770521548e-08, 
    betx=1.000000085403733*m, 
    bety=0.00160000031965417*m, 
    dispx=4.7831376053358213e-11*m,
    dispxp=-6.896462381454175e-11, 
    dispy=-0.0*m, 
    dispyp=0.0, 
    distrType="gausstwiss", 
    emitx=1e-12*m, 
    emity=1e-12*m, 
    energy=182.5000000007154*GeV, 
!    particle="e-",
    particle="proton",
    sigmaE=0.001;

d1: drift, l=1*m;
l1: line=(d1);
use, l1;

‌

[Laurie Nevay]

I've launched an issue with CLHEP:

https://its.cern.ch/jira/projects/CLHEP/issues/CLHEP-159?filter=allissues

The infinite loop comes in CLHEP::diagonalize. We use this to check if the matrix is positive definite and also CLHEP::RandMultiGauss would use it too even after our checks. There's is nothing obviously wrong to me about it and in any case a clear exception or failure could be easily handled instead of an infinite loop.

I looked at alternative implementations but given what's available in CLHEP I couldn't find an alternative way.

Instinctively, removing the off-diagonal terms (ie the dispersive ones) will likely get around this for now.

I'll update when I hear back from CLHEP.

[Laurie Nevay]

I’ve added print out for this so that it’s at least clear what’s going on and where it’s happening.

BDSBunch::SetEmittances> Geometric (x): 2.5e-08, Normalised (x): 4.89237834186e-05
BDSBunch::SetEmittances> Geometric (y): 2.5e-08, Normalised (y): 4.89237834186e-05
BDSBunchGaussian::CreateMultiGauss> checking 6x6 sigma matrix is positive definite
BDSBunchGaussian::CreateMultiGauss> confirmed: positive definite matrix

(last 2 lines here)

The infinite loop will happen in between them.