I made a workaround to add a correlation using "cov(t,p)", "cov(t,pz)" "cov(t,energy)" when in gaussian mode.

I haven't had a chance to look at it. It looks weird. It may be the alpha sign. I note that in longitudinal phase space things are different because faster particles arrive earlier, i.e. with lower t. On the other hand particles with higher px arrive with larger x). I haven't really had a chance to dig into the details. If it is working okay then maybe this is a false alarm...

Looks like I was generating a beam using pz as longitudinal variable, then calculating longitudinal optical functions using energy - which doesn't work if transverse emittance is non-zero