MPL: ImpactXParticleContainer.plot_phasespace()

[ax3l]

Add an interactive plotter to quickly check the current phase space of the particle bunch.

This is heavily inspired by Wake-T's bunch.show() functionality from APTools 💖

• depends on Pandas: ImpactXParticleContainer.to_df() #458
• depend on another pr that adds the beam min/max extent in space and momentum to reduced_beam_characteristics
  • rewrite histogram logic for scalability: histogram locally, reduce to rank 0 on results
• depend on another PR that tracks if the particles are currently in s or in t coordinates Particle Container: Track s or t #497
• add a helper function to generically scale unit prefixes to the closest metric prefix / engineering notation
  • https://github.com/KenKundert/quantiphy ?
  • https://github.com/hgrecco/pint ?
• write user-facing docs
• finalize test for interactive/non-interactive use

[cemitch99]

This looks great. Two comments:

• We might consider the powers of ten used for the axes. I often see plots with x and y in units of [mm], and with px and py in units of [mrad] for non-LPA generated beams.
• It would be good to double-check whether the emittance values are normalized or unnormalized. The notation suggests that the values are normalized.

[ax3l]

Ah yes, mm-mrad also makes sense for LPA beams actually.
What is the unit of our p_t? Is it in p_0 (reference momentum) or also (m)rad?

Yes, exactly, the data is normalized emittance.
Would you use a more specific notation in the plots to properly show this?

[cemitch99]

Ah yes, mm-mrad also makes sense for LPA beams actually. What is the unit of our p_t? Is it in p_0 (reference momentum) or also (m)rad?

Yes, exactly, the data is normalized emittance. Would you use a more specific notation in the plots to properly show this?

The emittance notation looks good. For p_t, the raw units are usually taken to be units of energy. In that case, our p_t is normalized by p_0c. We could write p_t [p_0c] or p_t/p_0c.

[cemitch99]

Regarding the discussion of scatterplots vs. density plots, see the attached screenshot from one of TraceWin's plotting packages (PlotWin). This combines colors (to denote phase space density) and points (to denote particles), resulting in a nice compromise. See, in particular, the particle distribution in the outer low-density halo.

[n01r]

Hmm, this still looks like a 2D histogram density to me, just with a higher resolution.

Here's a little snippet that creates a scatter plot with point coloring according to the density

# Creating a 2D histogram to calculate particle densities
hist, xedges, yedges = np.histogram2d(bunch.x, bunch.y, bins=250)
# hist: 2D array with density values for each bin
# xedges, yedges: bin edge values for x and y dimensions

# Function to get density values for each point from the histogram
def get_density_from_histogram(x, y, hist, xedges, yedges):
    # Find the indices of the bins to which each x and y value belongs
    x_idx = np.searchsorted(xedges, x, side='right') - 1
    y_idx = np.searchsorted(yedges, y, side='right') - 1

    # Ensure indices are within the valid range (0 to number of bins - 1)
    # This prevents index out-of-bounds errors
    x_idx = np.clip(x_idx, 0, hist.shape[0] - 1)
    y_idx = np.clip(y_idx, 0, hist.shape[1] - 1)

    # Return the density value for each (x, y) pair from the histogram
    return hist[x_idx, y_idx]

# Get density values for each point in the scatter plot
colors = get_density_from_histogram(bunch.x, bunch.y, hist, xedges, yedges)

# Plotting the scatter plot with density-based coloring
fig,ax = plt.subplots(1,1)
scatter = ax.scatter(bunch.x, bunch.y, c=colors, cmap='viridis', rasterized=True, s=1)  # Use 'viridis' colormap
plt.colorbar(scatter, ax=ax, label='Density in Bin')  # Add a color bar for density values
ax.set_title('Scatter Plot with Histogram-based Density Coloring')
ax.set_xlabel(r'$x$ (m)')
ax.set_ylabel(r'$y$ (m)')
ax.grid(True)

That produces:

[cemitch99]

Thanks; I guess this is always the case. Scatterplot = histogram with values 0 (no particles per bin) and 1 (> 0 particles per bin). Maybe the important features to balance are:

• if the bin size is too small (as measured by the local average number of particles/bin), the color density will be noisy
• if the bin size is too large, there will be poor visualization of the particles in low-density regions (eg, halo)

[ax3l]

Thank you both!

I agree on the analysis. Thus, the binning is a user-facing parameter and I just showed the default (50x50).

Note that doing a scatter plot is not as easily scalable to arbitrary particle numbers in the bunch as histogramming. (needs more thought to parallelize. We could add an option to scatter plot and let people opt-into the costs of serialization)

For this test, npart = 10000 is simply a bit too few particles and we get holes in phase space quickly.

50x50

100x100

200x200

[ax3l]

@cemitch99 @n01r shall we merge this for now? :)

[n01r]

Yea, totally, let's merge :)

[ax3l]

Just FYI, because it is related. Through CAMPA, @s-sajid-ali is currently working to ensure that ImpactX and Synergia can use https://github.com/ChristopherMayes/openPMD-beamphysics soon :)

[ax3l]

@n01r if you like to contribute the scatter logic (e.g., via a scatter, default: scatter=False) parameter, as a follow-up PR to this then that is definitely welcome! Just would make it default yet due to scalability.

[n01r]

Just FYI, because it is related. Through CAMPA, @s-sajid-ali is currently working to ensure that ImpactX and Synergia can use https://github.com/ChristopherMayes/openPMD-beamphysics soon :)

@ax3l, what's the advantage of hexagons? Why not dots?

You are referring to this here, I guess, https://github.com/ChristopherMayes/openPMD-beamphysics/blob/master/docs/examples/particle_examples.ipynb

[ax3l]

I was remembering this wrong - this is also a scatter plot (of dots) in openPMD-beamphysics.

I thought for a hot minute that was a histogram plot using hexagon bins.