Skip to content

4D round q-Gaussian generator #136

Open
elleanor-lamb wants to merge 14 commits into
xsuite:mainfrom
elleanor-lamb:main
Open

4D round q-Gaussian generator #136
elleanor-lamb wants to merge 14 commits into
xsuite:mainfrom
elleanor-lamb:main

Conversation

@elleanor-lamb

@elleanor-lamb elleanor-lamb commented Sep 15, 2025

Copy link
Copy Markdown

Description

Adds 4D q-gaussian generator to create normalised coordinates in 4D.
Beam is round and q is the same in x and y.
Example in xpart/examples/particles_generation/008_generate_q_gaussian.py

Checklist

Mandatory:

  • [x ] I have added tests to cover my changes
  • [x ] All the tests are passing, including my new ones
  • [x ] I described my changes in this PR description

Optional:

  • [x ] The code I wrote follows good style practices (see PEP 8 and PEP 20).
  • [ x] I have updated the docs in relation to my changes, if applicable
  • I have tested also GPU contexts

@elleanor-lamb

elleanor-lamb commented Sep 15, 2025

Copy link
Copy Markdown
Author

Need to add tests and test in GPU context. Would also appreciate some feedback on the numerical stability of generate_radial_distribution inside /transverse_generators/q_gaussian_round.py

@szymonlopaciuk szymonlopaciuk left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm sorry it took so long to answer. I have a couple of comments that may be helpful regarding your numerical stability question. Otherwise the PR is very nice and clean!

if abs(q - 1) < 1e-2:
term2 = -1 / (1 - q)
else:
term2 = gamma(q / (q - 1)) / gamma(1 / (q - 1))

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think you might not need the branch, since $\Gamma(x + 1) = x\Gamma(x)$:

$$ \frac{\Gamma\left(\frac{q}{q - 1}\right)}{\Gamma\left(\frac{1}{q - 1}\right)} = \frac{\Gamma\left(\frac{1 + q - 1}{q - 1}\right)}{\Gamma\left(\frac{1}{q - 1}\right)} = \frac{\Gamma\left(\frac{1}{q - 1} + 1\right)}{\Gamma\left(\frac{1}{q - 1}\right)} = \frac{\frac{1}{q - 1} \Gamma\left(\frac{1}{q - 1}\right)}{\Gamma\left(\frac{1}{q - 1}\right)} = \frac{1}{q - 1} $$

So the first one is not an approximation.

else:
term2 = gamma(q / (q - 1)) / gamma(1 / (q - 1))

term3 = (1 + beta * (q - 1) * F) ** (1 / (1 - q) - 3 / 2)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If I understand right, the stability problem is that for large beta and $q \rightarrow 1$ the base $\rightarrow \infty$ while the exponent $\rightarrow -\infty$? You could try computing with log and see if that helps:

$$ (1 + A)^B = \exp(\log((1 + A)^B)) = \exp(B \cdot \log(1 + A)) = \exp(B \cdot \text{log1p}(A)) $$

np.ndarray: CDF of g(F).
"""
delta_F = np.diff(F, prepend=0)
cdf = np.cumsum(g_F * delta_F) # todo: rewrite with scipy.integrate.quad(g_F, -np.inf, np.inf)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Still todo?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants