In [1]:

# Useful for debugging
%load_ext autoreload
%autoreload 2
%config InlineBackend.figure_format = 'retina'

# Useful for debugging %load_ext autoreload %autoreload 2 %config InlineBackend.figure_format = 'retina'

In [2]:

from distgen import Generator
from distgen.plot import plot_dist2d, plot_dist1d, plot_current_profile

from matplotlib import pyplot as plt
import yaml

from distgen import Generator from distgen.plot import plot_dist2d, plot_dist1d, plot_current_profile from matplotlib import pyplot as plt import yaml

In [3]:

gen = Generator(verbose=0)
units = {
    "x": "mm",
    "y": "mm",
    "z": "mm",
    "px": "keV/c",
    "py": "keV/c",
    "t": "ps",
    "q": "pC",
    "I": "A",
}

gen = Generator(verbose=0) units = { "x": "mm", "y": "mm", "z": "mm", "px": "keV/c", "py": "keV/c", "t": "ps", "q": "pC", "I": "A", }

Radial Distributions¶

---

Radial Uniform¶

The radial uniform distribution is given by

$\rho_r(r)=\frac{2}{r_R^2-r_L^2}$ for $0 \leq r_L < r < r_R$ and zero otherewise.

Typical usage is shown below for a uniformly filled beam as well as a ring:

In [4]:

gen = Generator("data/rad.uniform.in.yaml", verbose=0)
gen["r_dist:min_r"] = 0
print("\nYAML input:")
print(yaml.dump({"r_dist": gen.input["r_dist"]}))

plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal", title_on=True);

gen = Generator("data/rad.uniform.in.yaml", verbose=0) gen["r_dist:min_r"] = 0 print("\nYAML input:") print(yaml.dump({"r_dist": gen.input["r_dist"]})) plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal", title_on=True);

YAML input:
r_dist:
  max_r:
    units: millimeter
    value: 2
  min_r:
    units: millimeter
    value: 0
  type: radial_uniform

In [5]:

gen = Generator("data/rad.uniform.in.yaml", verbose=0)

print("\nYAML input:")
print(yaml.dump({"r_dist": gen.input["r_dist"]}))

plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal", title_on=True);

gen = Generator("data/rad.uniform.in.yaml", verbose=0) print("\nYAML input:") print(yaml.dump({"r_dist": gen.input["r_dist"]})) plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal", title_on=True);

YAML input:
r_dist:
  max_r:
    units: millimeter
    value: 2
  min_r:
    units: millimeter
    value: 1
  type: radial_uniform

Radial Gaussian (with Truncation)¶

The radial normal distribution including truncation(s) has a probability function given by

$\rho_r(r) = \frac{1}{\sigma^2}\frac{\phi(r/\sigma)}{\phi\left(\frac{r_L}{\sigma}\right)-\phi\left(\frac{r_R}{\sigma}\right)} $ for $0 \leq r_L \leq r \leq r_R$ and zero everywhere else.

In this expresion $\phi(\xi) = \frac{1}{2\pi}\exp\left(-\xi^2/2\right)$ is the canonical raidial normal distirbution (no truncation), and the scale parameter $\sigma$ follows from the product of two normal distributions in $x$ and $y$ when $\sigma=\sigma_x=\sigma_y$. The corresponding CDF is given by

$P(r)= \frac{\phi\left(\frac{r_L}{\sigma}\right)-\phi\left(\frac{r}{\sigma}\right)}{\phi\left(\frac{r_L}{\sigma}\right)-\phi\left(\frac{r_R}{\sigma}\right)} $ for $0 \leq r_L \leq r$.

The corresponding first and second moments are:

$\langle r\rangle = \frac{\frac{r_L}{\sigma}\phi\left(\frac{r_L}{\sigma}\right) -\frac{r_R}{\sigma}\phi\left(\frac{r_R}{\sigma}\right) +\frac{1}{2\sqrt{2\pi}}\left( \text{erf}\left(\frac{r_R}{\sigma\sqrt{2}}\right) - \text{erf}\left(\frac{r_L}{\sigma\sqrt{2}}\right) \right) } {\phi\left(\frac{r_L}{\sigma}\right)-\phi\left(\frac{r_R}{\sigma}\right)}$,

$r_{rms} = \sqrt{ 2\sigma^2 + r_L^2 - \frac{(r_R^2-r_L^2)\phi(r_R/\sigma)}{\phi\left(\frac{r_L}{\sigma}\right)-\phi\left(\frac{r_R}{\sigma}\right)} }$.

Note that in the limits $r_L\rightarrow 0$ and $r_R -> \infty$ the above expressions reduce to the underlying radial normal distribution:

$\rho_r(r)\rightarrow \frac{\phi\left(\frac{r}{\sigma}\right)}{\sigma^2}$, $P(r)\rightarrow 1 - \phi\left(\frac{r}{\sigma}\right)$, $\langle r\rangle\rightarrow \sqrt{\frac{\pi}{2}}\sigma$, and $r_{rms}\rightarrow \sqrt{2}\sigma$. This limiting case is shown first below.

Here three examples are provided: a non-truncated Gaussian with $\sigma_{x,y}$ = 2 mm, followed by a truncated Gaussian using a truncation fraction to specify the maximal radial cut, followed by a truncated Gaussian using n_sigma_cutoff to specify the radial cut.

Example 1: canonical Gaussian with $\sigma_{x,y}$ = 2mm

In [6]:

# Radial gaussian distribution with sigma_xy 2 mm
gen = Generator("data/rad.gaussian.in.yaml", verbose=0)

print("\nYAML input:")
print(yaml.dump({"r_dist": gen.input["r_dist"]}))

# X-Y Plot
plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal", title_on=True);

# Radial gaussian distribution with sigma_xy 2 mm gen = Generator("data/rad.gaussian.in.yaml", verbose=0) print("\nYAML input:") print(yaml.dump({"r_dist": gen.input["r_dist"]})) # X-Y Plot plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal", title_on=True);

YAML input:
r_dist:
  sigma_xy:
    units: millimeter
    value: 2
  type: radial_gaussian

Example 2: Truncated Gaussian with a pinhole radius of 1 mm, clipped at 50% intensity

In [7]:

# Truncated radial gaussian distribution with pihole diameter of 2 mm, gaussian clipped at 50% intensity
gen = Generator("data/rad.trunc.gaussian.in.yaml", verbose=0)

print("\nYAML input:")
print(yaml.dump({"r_dist": gen.input["r_dist"]}))

# X-Y Plot
plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal", title_on=True);

# Truncated radial gaussian distribution with pihole diameter of 2 mm, gaussian clipped at 50% intensity gen = Generator("data/rad.trunc.gaussian.in.yaml", verbose=0) print("\nYAML input:") print(yaml.dump({"r_dist": gen.input["r_dist"]})) # X-Y Plot plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal", title_on=True);

YAML input:
r_dist:
  truncation_fraction:
    units: dimensionless
    value: 0.5
  truncation_radius:
    units: millimeter
    value: 1
  type: rg

Example 3: truncated Gaussian using underlying $\sigma_{x,y}$ = 2 mm and n_sigma_cutoff = 1

In [8]:

# Radial gaussian distribution with cut using n_sigma_cutoff
gen = Generator("data/rad.gaussian.in.yaml", verbose=0)
gen["r_dist:n_sigma_cutoff"] = 1

print("\nYAML input:")
print(yaml.dump({"r_dist": gen.input["r_dist"]}))

# X-Y Plot
plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal", title_on=True);

# Radial gaussian distribution with cut using n_sigma_cutoff gen = Generator("data/rad.gaussian.in.yaml", verbose=0) gen["r_dist:n_sigma_cutoff"] = 1 print("\nYAML input:") print(yaml.dump({"r_dist": gen.input["r_dist"]})) # X-Y Plot plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal", title_on=True);

YAML input:
r_dist:
  n_sigma_cutoff: 1
  sigma_xy:
    units: millimeter
    value: 2
  type: radial_gaussian

Radial File¶

In [9]:

# Radial distribution file: a radial gaussian with cuts in it
gen = Generator("data/rad.file.in.yaml", verbose=0)

print("\nYAML input:")
print(yaml.dump({"r_dist": gen.input["r_dist"]}))

# X-Y Plot
plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal");

# Radial distribution file: a radial gaussian with cuts in it gen = Generator("data/rad.file.in.yaml", verbose=0) print("\nYAML input:") print(yaml.dump({"r_dist": gen.input["r_dist"]})) # X-Y Plot plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal");

YAML input:
r_dist:
  file: /home/runner/work/distgen/distgen/docs/examples/data/cutgauss.rad.txt
  type: radfile
  units: mm

Radial Tukey¶

In [10]:

# Radial Tukey distribution
gen = Generator("data/tukey.beer.can.in.yaml", verbose=0)

print("\nYAML Input:")
print(yaml.dump({"r_dist": gen.input["r_dist"]}))

plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal", title_on=True);

# Radial Tukey distribution gen = Generator("data/tukey.beer.can.in.yaml", verbose=0) print("\nYAML Input:") print(yaml.dump({"r_dist": gen.input["r_dist"]})) plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal", title_on=True);

YAML Input:
r_dist:
  length:
    units: millimeter
    value: 5
  ratio:
    units: dimensionless
    value: 0.25
  type: radial_tukey

Radial Super Gaussian¶

This implements a radial version of the Super Gaussian function. Here the radial function takes the form:

$2\pi\rho(r;\lambda,p) = \frac{1}{\Gamma\left(1+\frac{1}{p}\right)\lambda^2} \exp\left[-\left(\frac{r^2}{2\lambda^2}\right)^p\right]$.

The first and (rms) second moment of the distribution are given by:

$\langle r\rangle = \frac{2\sqrt{2}}{3}\frac{\Gamma\left(1+\frac{3}{2p}\right)}{\Gamma\left(1+\frac{1}{p}\right)}\lambda$,

$r_{\text{rms}} = \sqrt{\frac{\Gamma\left(1+\frac{2}{p}\right)}{\Gamma\left(1+\frac{1}{p}\right)}}\lambda$.

Example 1: hyper-Gaussian (uniform) limit (power $\rightarrow\infty$): $\sigma_{x,y}$ = 1 mm

In [11]:

gen = Generator("data/SG.beer.can.in.yaml", verbose=0)

print("\nYAML Input:")
print(yaml.dump({"r_dist": gen.input["r_dist"]}))

plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal", title_on=True);

gen = Generator("data/SG.beer.can.in.yaml", verbose=0) print("\nYAML Input:") print(yaml.dump({"r_dist": gen.input["r_dist"]})) plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal", title_on=True);

YAML Input:
r_dist:
  alpha:
    units: dimensionless
    value: 0
  sigma_xy:
    units: millimeter
    value: 1
  type: rsg

Example 2: Gaussian limit (power = 1): $\sigma_{x,y}$ = 1 mm

In [12]:

gen = Generator("data/SG.beer.can.in.yaml", verbose=0)
gen["r_dist:alpha"] = 1.0

print("\nYAML Input:")
print(yaml.dump({"r_dist": gen.input["r_dist"]}))

plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal", title_on=True);

gen = Generator("data/SG.beer.can.in.yaml", verbose=0) gen["r_dist:alpha"] = 1.0 print("\nYAML Input:") print(yaml.dump({"r_dist": gen.input["r_dist"]})) plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal", title_on=True);

YAML Input:
r_dist:
  alpha:
    units: dimensionless
    value: 1.0
  sigma_xy:
    units: millimeter
    value: 1
  type: rsg

Example 3: sub-Gaussian (power < 1): $\sigma_{x,y} = 1$ mm, $p$ = 0.5

In [13]:

with open("data/SG.beer.can.in.yaml", "r") as fid:
    inputs = yaml.safe_load(fid)

del inputs["r_dist"]["alpha"]
inputs["r_dist"]["p"] = {"value": 0.5, "units": ""}

gen = Generator(inputs, verbose=0)

print("\nYAML Input:")
print(yaml.dump({"r_dist": gen.input["r_dist"]}))

plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal", title_on=True);

with open("data/SG.beer.can.in.yaml", "r") as fid: inputs = yaml.safe_load(fid) del inputs["r_dist"]["alpha"] inputs["r_dist"]["p"] = {"value": 0.5, "units": ""} gen = Generator(inputs, verbose=0) print("\nYAML Input:") print(yaml.dump({"r_dist": gen.input["r_dist"]})) plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal", title_on=True);

YAML Input:
r_dist:
  p:
    units: dimensionless
    value: 0.5
  sigma_xy:
    units: millimeter
    value: 1
  type: rsg

Radial Deformable¶

In [14]:

gen = Generator("data/deformable.in.yaml", verbose=0)

print("\nYAML Input:")
print(yaml.dump({"r_dist": gen.input["r_dist"]}))

plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal", title_on=True);

gen = Generator("data/deformable.in.yaml", verbose=0) print("\nYAML Input:") print(yaml.dump({"r_dist": gen.input["r_dist"]})) plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], axis="equal", title_on=True);

YAML Input:
r_dist:
  alpha:
    units: dimensionless
    value: 0.1
  sigma_xy:
    units: millimeter
    value: 1
  slope_fraction:
    units: dimensionless
    value: -1
  type: dr

Radial Interpolation¶

In [15]:

gen = Generator("data/interp_ex.yaml", verbose=0)
print("\nYAML Input:")
print(yaml.dump({"r_dist": gen.input["r_dist"]}))
PG = gen.run()

gen = Generator("data/interp_ex.yaml", verbose=0) print("\nYAML Input:") print(yaml.dump({"r_dist": gen.input["r_dist"]})) PG = gen.run()

YAML Input:
r_dist:
  Pr:
    pt01: 0.0
    pt02: 1.1
    pt03: 3.2
    pt04: 8.0
    pt05: 9.0
    pt06: 8.0
    pt07: 12.0
    pt08: 3.0
    pt09: 2.0
    pt10: 3.0
    pt11: 5.0
    pt12: 1.0
  method: spline
  type: radial_interpolation

In [16]:

PG.plot("r")

PG.plot("r")

In [17]:

PG.plot("x", "y")

PG.plot("x", "y")

In [18]:

PG["sigma_x"], PG["sigma_y"]

PG["sigma_x"], PG["sigma_y"]

Out[18]:

(np.float64(0.0012), np.float64(0.0012))

In [19]:

from distgen.tools import SUPPORTED_IMAGE_EXTENSIONS

from distgen.tools import SUPPORTED_IMAGE_EXTENSIONS

2D Distributions¶

---

2D File¶

The user can load images for use in generating particles. Supported file formats include:

In [20]:

SUPPORTED_IMAGE_EXTENSIONS

SUPPORTED_IMAGE_EXTENSIONS

Out[20]:

['.jpeg', '.jpg', '.png', '.tiff', '.dicom', '.dcm']

The example below uses .jpeg.

In [21]:

gen = Generator("data/jpeg.image.in.yaml", verbose=0)
gen.run()

print("\nYAML input:")
print(yaml.dump({"xy_dist": gen.input["xy_dist"]}))

gen.particles.plot("x", "y")

gen = Generator("data/jpeg.image.in.yaml", verbose=0) gen.run() print("\nYAML input:") print(yaml.dump({"xy_dist": gen.input["xy_dist"]})) gen.particles.plot("x", "y")

YAML input:
xy_dist:
  file: /home/runner/work/distgen/distgen/docs/examples/data/seaturtle.jpeg
  max_x:
    units: millimeter
    value: 1.0
  max_y:
    units: millimeter
    value: 1.0
  min_x:
    units: millimeter
    value: -1.0
  min_y:
    units: millimeter
    value: -1.0
  threshold: 0.0
  type: file2d

In [22]:

gen = Generator("data/dcm.image.in.yaml", verbose=0)
gen.run()

print("\nYAML input:")
print(yaml.dump({"xy_dist": gen.input["xy_dist"]}))

gen.particles.plot("x", "y")

gen = Generator("data/dcm.image.in.yaml", verbose=0) gen.run() print("\nYAML input:") print(yaml.dump({"xy_dist": gen.input["xy_dist"]})) gen.particles.plot("x", "y")

YAML input:
xy_dist:
  file: /home/runner/work/distgen/distgen/docs/examples/data/image-00319.dcm
  max_x:
    units: millimeter
    value: 1
  max_y:
    units: millimeter
    value: 1
  min_x:
    units: millimeter
    value: -1
  min_y:
    units: millimeter
    value: -1
  threshold: 0.0
  type: file2d

In addition to handling these standard images, distgen also supports the legacy ascii format from the original version of the code. The example below makes use of this legacy format:

In [23]:

# 2D XY distribution: Lase Image Example
# Real laser spot image - using 1 million particles
gen = Generator("data/laser.image.in.yaml", verbose=0)

print("\nYAML input:")
print(yaml.dump({"xy_dist": gen.input["xy_dist"]}))

# X-Y Plot
plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], nbins=100, axis="equal");

# 2D XY distribution: Lase Image Example # Real laser spot image - using 1 million particles gen = Generator("data/laser.image.in.yaml", verbose=0) print("\nYAML input:") print(yaml.dump({"xy_dist": gen.input["xy_dist"]})) # X-Y Plot plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], nbins=100, axis="equal");

YAML input:
xy_dist:
  file: /home/runner/work/distgen/distgen/docs/examples/data/laser.image.txt
  type: file2d

2D Image¶

If the user has an image dataset already loaded into python, they can directly pass this to Distgen using the Image2d class. An example of how to use this in practice is given below:

In [24]:

from distgen.tools import read_2d_file

xs, ys, Pxy, xstr, ystr = read_2d_file("data/laser.image.txt")

inputs = {
    "n_particle": 200000,
    "species": "electron",
    "total_charge": {"value": 1, "units": "pC"},
    "random_type": "hammersley",
    "start": {"MTE": {"units": "meV", "value": 150}, "type": "cathode"},
    "xy_dist": {
        "type": "image2d",
        "min_x": {"value": -2.0, "units": "mm"},
        "max_x": {"value": +2.0, "units": "mm"},
        "min_y": {"value": -2.0, "units": "mm"},
        "max_y": {"value": +2.0, "units": "mm"},
        "P": Pxy,
    },
}

gen = Generator(inputs, verbose=0)

gen["xy_dist"]
# plot_dist2d(gen.beam(), 'x', units['x'], 'y', units['y'], nbins=100, axis='equal');

from distgen.tools import read_2d_file xs, ys, Pxy, xstr, ystr = read_2d_file("data/laser.image.txt") inputs = { "n_particle": 200000, "species": "electron", "total_charge": {"value": 1, "units": "pC"}, "random_type": "hammersley", "start": {"MTE": {"units": "meV", "value": 150}, "type": "cathode"}, "xy_dist": { "type": "image2d", "min_x": {"value": -2.0, "units": "mm"}, "max_x": {"value": +2.0, "units": "mm"}, "min_y": {"value": -2.0, "units": "mm"}, "max_y": {"value": +2.0, "units": "mm"}, "P": Pxy, }, } gen = Generator(inputs, verbose=0) gen["xy_dist"] # plot_dist2d(gen.beam(), 'x', units['x'], 'y', units['y'], nbins=100, axis='equal');

Out[24]:

{'type': 'image2d',
 'min_x': {'value': -2.0, 'units': 'millimeter'},
 'max_x': {'value': 2.0, 'units': 'millimeter'},
 'min_y': {'value': -2.0, 'units': 'millimeter'},
 'max_y': {'value': 2.0, 'units': 'millimeter'},
 'P': {'value': array([[0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         ...,
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.]]),
  'units': '1 / millimeter ** 2'}}

Note that this allows the user to manipulate the 2d distribution and length scales:

In [25]:

gen = Generator(inputs, verbose=0)

P2 = Pxy.magnitude
P2 = P2.transpose()

# Cut out a rectangle in the center
P2[236 - 50 : 225 + 50, 236 - 50 : 225 + 50] = 0

gen["xy_dist:P"] = P2
gen["xy_dist:min_x:value"], gen["xy_dist:max_x:value"] = -5000, +5000
gen["xy_dist:min_x:units"], gen["xy_dist:max_x:units"] = "um", "um"
gen["xy_dist:min_y:units"], gen["xy_dist:max_y:units"] = "mm", "mm"

plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], nbins=100, axis="equal");

gen = Generator(inputs, verbose=0) P2 = Pxy.magnitude P2 = P2.transpose() # Cut out a rectangle in the center P2[236 - 50 : 225 + 50, 236 - 50 : 225 + 50] = 0 gen["xy_dist:P"] = P2 gen["xy_dist:min_x:value"], gen["xy_dist:max_x:value"] = -5000, +5000 gen["xy_dist:min_x:units"], gen["xy_dist:max_x:units"] = "um", "um" gen["xy_dist:min_y:units"], gen["xy_dist:max_y:units"] = "mm", "mm" plot_dist2d(gen.beam(), "x", units["x"], "y", units["y"], nbins=100, axis="equal");

Uniform Laser Speckle¶

This distribution provides a very simplistic method of simulating laser speckle. The input distribution parameter sigma controls the size of the speckles. Sigma is the size of a blurring Gaussian which is used to scale the FFT of the image before transforming back. Smaller $\sigma$ gives larger speckles. In addition the user can specify a pixel_threshold [0 1] that will zero a given pixel if the resulting speckle pattern amplitude is below the threshold.

In [26]:

gen = Generator("data/uniform_speckle.in.yaml")

print("\nYAML input:")
print(yaml.dump({"xy_dist": gen.input["xy_dist"]}))

gen.run()
gen.particles.plot("x", "y")

gen = Generator("data/uniform_speckle.in.yaml") print("\nYAML input:") print(yaml.dump({"xy_dist": gen.input["xy_dist"]})) gen.run() gen.particles.plot("x", "y")

YAML input:
xy_dist:
  max_x:
    units: millimeter
    value: 1
  max_y:
    units: millimeter
    value: 2
  min_x:
    units: millimeter
    value: -1
  min_y:
    units: millimeter
    value: -2
  pixel_threshold: 0.25
  random_seed: 42
  sigma: 100
  type: uniform_laser_speckle

In [27]:

gen["xy_dist:sigma"] = 10
gen.run()
gen.particles.plot("x", "y")

gen["xy_dist:sigma"] = 10 gen.run() gen.particles.plot("x", "y")

Product 2D¶

It is possible to multiply two or more 2d or radial distributions. This is useful for creating particles from a laser and folding in a cathod QE map. The example shows a truncated Gaussian laser profile combined with a simple checker board QE map:

In [28]:

import numpy as np

mask = np.kron([[1, 0] * 4, [0, 1] * 4] * 4, np.ones((50, 50)))

ys = np.linspace(-4, +4, mask.shape[0])
xs = np.linspace(-2, +2, mask.shape[1])

fig, ax = plt.subplots(1, 3, constrained_layout=True, figsize=(8, 4))

ax[0].imshow(mask, extent=[xs.min(), xs.max(), ys.min(), ys.max()])
ax[0].set_aspect("equal", adjustable="box")
ax[0].set_title("QE Map")
ax[0].set_xlabel("x (mm)")
ax[0].set_ylabel("y (mm)")

laser_inputs = {
    "n_particle": 200000,
    "species": "electron",
    "total_charge": {"value": 1, "units": "pC"},
    "random_type": "hammersley",
    "start": {"MTE": {"units": "meV", "value": 150}, "type": "cathode"},
    "r_dist": {
        "type": "rg",
        "truncation_fraction": {"units": "dimensionless", "value": 0.5},
        "truncation_radius": {"units": "mm", "value": 2},
    },
}

gen = Generator(laser_inputs)
beam1 = gen.beam()

qe_and_laser_inputs = {
    "n_particle": 200000,
    "species": "electron",
    "total_charge": {"value": 1, "units": "pC"},
    "random_type": "hammersley",
    "start": {"MTE": {"units": "meV", "value": 150}, "type": "cathode"},
    "xy_dist": {
        "type": "product",
        "dists": {
            "QE": {
                "type": "image2d",
                "x": {"value": xs, "units": "mm"},
                "y": {"value": ys, "units": "mm"},
                "P": mask,
            },
            "Laser": {
                "type": "rg",
                "truncation_fraction": {"units": "dimensionless", "value": 0.5},
                "truncation_radius": {"units": "mm", "value": 2},
            },
        },
    },
}

gen = Generator(qe_and_laser_inputs)
beam2 = gen.beam()

plot_dist2d(beam1, "x", units["x"], "y", units["y"], nbins=100, axis="equal", ax=ax[1])
ax[1].set_title("Laser Profile")

plot_dist2d(beam2, "x", units["x"], "y", units["y"], nbins=100, axis="equal", ax=ax[2])
ax[2].set_title("Laser & QE Profile");

import numpy as np mask = np.kron([[1, 0] * 4, [0, 1] * 4] * 4, np.ones((50, 50))) ys = np.linspace(-4, +4, mask.shape[0]) xs = np.linspace(-2, +2, mask.shape[1]) fig, ax = plt.subplots(1, 3, constrained_layout=True, figsize=(8, 4)) ax[0].imshow(mask, extent=[xs.min(), xs.max(), ys.min(), ys.max()]) ax[0].set_aspect("equal", adjustable="box") ax[0].set_title("QE Map") ax[0].set_xlabel("x (mm)") ax[0].set_ylabel("y (mm)") laser_inputs = { "n_particle": 200000, "species": "electron", "total_charge": {"value": 1, "units": "pC"}, "random_type": "hammersley", "start": {"MTE": {"units": "meV", "value": 150}, "type": "cathode"}, "r_dist": { "type": "rg", "truncation_fraction": {"units": "dimensionless", "value": 0.5}, "truncation_radius": {"units": "mm", "value": 2}, }, } gen = Generator(laser_inputs) beam1 = gen.beam() qe_and_laser_inputs = { "n_particle": 200000, "species": "electron", "total_charge": {"value": 1, "units": "pC"}, "random_type": "hammersley", "start": {"MTE": {"units": "meV", "value": 150}, "type": "cathode"}, "xy_dist": { "type": "product", "dists": { "QE": { "type": "image2d", "x": {"value": xs, "units": "mm"}, "y": {"value": ys, "units": "mm"}, "P": mask, }, "Laser": { "type": "rg", "truncation_fraction": {"units": "dimensionless", "value": 0.5}, "truncation_radius": {"units": "mm", "value": 2}, }, }, }, } gen = Generator(qe_and_laser_inputs) beam2 = gen.beam() plot_dist2d(beam1, "x", units["x"], "y", units["y"], nbins=100, axis="equal", ax=ax[1]) ax[1].set_title("Laser Profile") plot_dist2d(beam2, "x", units["x"], "y", units["y"], nbins=100, axis="equal", ax=ax[2]) ax[2].set_title("Laser & QE Profile");

In [29]:

inputs = """
n_particle: 1000000
species: electron
xy_dist:
  type: product
  dists:
    speckle:
      max_x: +2 mm 
      min_x: -2 mm
      max_y: +2 mm 
      min_y: -2 mm
      sigma: 100
      pixel_threshold: 0.25
      random_seed: 42
      type: uniform_laser_speckle
    uniform:
      max_r: 2 mm 
      min_r: 0 mm
      type: radial_uniform     
random:
  type: hammersley
start:
  MTE: 150 meV
  type: cathode
t_dist:
  avg_t: 0.5 ns
  sigma_t: 2 ns
  type: gaussian
total_charge: 10 pC
"""

gen = Generator(inputs)
gen.run()
gen.particles.plot("x", "y", figsize=(5, 5))

inputs = """ n_particle: 1000000 species: electron xy_dist: type: product dists: speckle: max_x: +2 mm min_x: -2 mm max_y: +2 mm min_y: -2 mm sigma: 100 pixel_threshold: 0.25 random_seed: 42 type: uniform_laser_speckle uniform: max_r: 2 mm min_r: 0 mm type: radial_uniform random: type: hammersley start: MTE: 150 meV type: cathode t_dist: avg_t: 0.5 ns sigma_t: 2 ns type: gaussian total_charge: 10 pC """ gen = Generator(inputs) gen.run() gen.particles.plot("x", "y", figsize=(5, 5))

In [30]:

gen["xy_dist:dists:speckle:sigma"] = 5
gen.run()
gen.particles.plot("x", "y", figsize=(5, 5))

gen["xy_dist:dists:speckle:sigma"] = 5 gen.run() gen.particles.plot("x", "y", figsize=(5, 5))

Superposition 2D¶

TODO

1D Distributions¶

---

Uniform¶

The uniform distirbuition is defined by a probability distribution function:

$\rho(x) = \frac{1}{b-a}$ for $a\leq x\leq b$ and zero elsewhere.

The corresponding CDF is

$P(x) = \frac{x-a}{b-a}$ for $a\leq x\leq b$ and zero elsewhere.

The first and second moments of this distribution are:

$\langle x \rangle = \frac{1}{2}(a+b)$ and $\sigma_x = \frac{b-a}{\sqrt{12}}$.

It is possible to set a uniform distribution using the min_{var} and max_{var} input variables:

In [31]:

gen = Generator("data/x.uniform.y.file.t.gaussian.in.yaml", verbose=0)

print("\nYAML Input:")
print(yaml.dump({"x_dist": gen.input["x_dist"]}))

plot_dist1d(gen.beam(), "x", units["y"], nbins=100, title_on=True)

gen = Generator("data/x.uniform.y.file.t.gaussian.in.yaml", verbose=0) print("\nYAML Input:") print(yaml.dump({"x_dist": gen.input["x_dist"]})) plot_dist1d(gen.beam(), "x", units["y"], nbins=100, title_on=True)

YAML Input:
x_dist:
  max_x:
    units: millimeter
    value: 2
  min_x:
    units: millimeter
    value: -1
  type: uniform

The user may also specify the distribution using the sigma_{var} and avg_{var} input params:

In [32]:

import numpy as np

gen.input["t_dist"] = {
    "type": "uniform",
    "sigma_t": {"value": float(1 / np.sqrt(12)), "units": "ps"},
    "avg_t": "2 ps",
}

gen = Generator(gen.input, verbose=0)

print("\nYAML Input:")
print(yaml.dump({"t_dist": gen.input["t_dist"]}))

plot_dist1d(gen.beam(), "t", units["t"], nbins=100, title_on=True)

import numpy as np gen.input["t_dist"] = { "type": "uniform", "sigma_t": {"value": float(1 / np.sqrt(12)), "units": "ps"}, "avg_t": "2 ps", } gen = Generator(gen.input, verbose=0) print("\nYAML Input:") print(yaml.dump({"t_dist": gen.input["t_dist"]})) plot_dist1d(gen.beam(), "t", units["t"], nbins=100, title_on=True)

YAML Input:
t_dist:
  avg_t:
    units: picosecond
    value: -0.5
  sigma_t:
    units: picosecond
    value: 1.0
  type: gaussian

Linear¶

The linear distirbuition is defined by a probability distribution function:

$\rho(x) = \frac{p_b-p_a}{b-a}(x-a) + p_a$ for $a\leq x\leq b$ and zero elsewhere.

The corresponding CDF is

$P(x) = \frac{1}{2}\frac{p_b-p_a}{b-a}(x-a)$ for $a\leq x\leq b$ and zero elsewhere.

It is possible to set a uniform distribution using the min_{var} and max_{var}(equivalent to $a$ and $b$ above) input variables as well as the slope_fraction $\in[-1,1]$.

1D File¶

Disgten supports importing a 1D PDF saved in column form in. The input form of the file should have space separated headers such as $x$ and $Px$, with corresponding column data below it. The PDF is normalized numerically using the numpy.trapz numerical integration routine. The CDF is computed using the scipy.cumtrapz cumulative numerical intgration routine.

The following example shows a gaussian PDF with cuts added to it.

In [33]:

gen = Generator("data/x.uniform.y.file.t.gaussian.in.yaml", verbose=0)

print("\nYAML Input:")
print(yaml.dump({"y_dist": gen.input["y_dist"]}))
plot_dist1d(gen.beam(), "y", units["y"], nbins=1000)

gen = Generator("data/x.uniform.y.file.t.gaussian.in.yaml", verbose=0) print("\nYAML Input:") print(yaml.dump({"y_dist": gen.input["y_dist"]})) plot_dist1d(gen.beam(), "y", units["y"], nbins=1000)

YAML Input:
y_dist:
  file: /home/runner/work/distgen/distgen/docs/examples/data/cutgauss.1d.txt
  type: file1d
  units: mm

Gaussian (with Truncation)¶

If not specified, the Gaussian is truncated at $\mu\pm5\sigma$:

In [34]:

gen = Generator("data/x.uniform.y.file.t.gaussian.in.yaml", verbose=0)

print("\nYAML Input:")
print(yaml.dump({"t_dist": gen.input["t_dist"]}))

plot_dist1d(gen.beam(), "t", units["t"], nbins=1000)

gen = Generator("data/x.uniform.y.file.t.gaussian.in.yaml", verbose=0) print("\nYAML Input:") print(yaml.dump({"t_dist": gen.input["t_dist"]})) plot_dist1d(gen.beam(), "t", units["t"], nbins=1000)

YAML Input:
t_dist:
  avg_t:
    units: picosecond
    value: -0.5
  sigma_t:
    units: picosecond
    value: 1
  type: gaussian

The user can specify n_sigma_cutoff to symmetrically truncate at $ \mu\pm N_{\sigma~\text{cutoff}}\times\sigma$:

In [35]:

gen = Generator("data/x.uniform.y.file.t.gaussian.in.yaml", verbose=0)
gen["t_dist:avg_t"] = 1
gen["t_dist:n_sigma_cutoff"] = 2

print("\nYAML Input:")
print(yaml.dump({"t_dist": gen.input["t_dist"]}))

plot_dist1d(gen.beam(), "t", units["t"], nbins=1000)

gen = Generator("data/x.uniform.y.file.t.gaussian.in.yaml", verbose=0) gen["t_dist:avg_t"] = 1 gen["t_dist:n_sigma_cutoff"] = 2 print("\nYAML Input:") print(yaml.dump({"t_dist": gen.input["t_dist"]})) plot_dist1d(gen.beam(), "t", units["t"], nbins=1000)

YAML Input:
t_dist:
  avg_t:
    units: picosecond
    value: 1
  n_sigma_cutoff: 2
  sigma_t:
    units: picosecond
    value: 1
  type: gaussian

Super Gaussian¶

In additional to the regular Gaussian function, it is also possible to sample a super-Gaussian distribution defined by

$\rho(x; \lambda, p) = \frac{1}{2\sqrt{2}\Gamma\left(1+\frac{1}{2p}\right)\lambda } \exp\left[-\left(\frac{(x-\mu)^2 }{2\lambda^2}\right)^{p}\right]$

Here $\lambda$ is the length scale and $p$ is the power of the super-Gaussian. Note when $p=1$ this expression reduces to a Normal distirbution, in which case $\sigma_x=\lambda$. As $p\rightarrow\infty$ the distribution reduces to a flat-top (uniform). The full range of powers is given by $p\in\left(0,\infty\right]$.

The first and second moments of the distribution are given by:

$\langle x\rangle = \mu$, and $\sigma_x = \left(\frac{2\Gamma\left(1+\frac{3}{2p}\right)}{3\Gamma\left(1+\frac{1}{2p}\right)}\right)^{1/2}\lambda$.

Often, it is convenient to scan the distribution from the uniform limit to the Gaussian limit. To do some, the input $p$ can be parameterized by $\alpha\in[0,1]$ where $p = 1/\alpha$. Here $\alpha=0$ corresponds to a flat-top (uniform) and $\alpha=1$ corresponds to a Gaussian. Examples of both types of usage are shown below.\

Example 1: specifying $\alpha = 1/p \in[0,1]$:

In [36]:

# Super Gaussian distribution
gen = Generator("data/SG.beer.can.in.yaml", verbose=0)

print("\nYAML Input:")
print(yaml.dump({"t_dist": gen.input["t_dist"]}))

plot_current_profile(gen.beam(), units["t"], units["I"], title_on=True)

# Super Gaussian distribution gen = Generator("data/SG.beer.can.in.yaml", verbose=0) print("\nYAML Input:") print(yaml.dump({"t_dist": gen.input["t_dist"]})) plot_current_profile(gen.beam(), units["t"], units["I"], title_on=True)

YAML Input:
t_dist:
  alpha:
    units: dimensionless
    value: 0.5
  sigma_t:
    units: picosecond
    value: 8
  type: sg

Tukey¶

In [37]:

# Tukey distribution
gen = Generator("data/tukey.beer.can.in.yaml", verbose=0)

print("\nYAML Input:")
print(yaml.dump({"t_dist": gen.input["t_dist"]}))

plot_current_profile(gen.beam(), units["t"], units["I"], nbins=500, title_on=True)

# Tukey distribution gen = Generator("data/tukey.beer.can.in.yaml", verbose=0) print("\nYAML Input:") print(yaml.dump({"t_dist": gen.input["t_dist"]})) plot_current_profile(gen.beam(), units["t"], units["I"], nbins=500, title_on=True)

YAML Input:
t_dist:
  length:
    units: picosecond
    value: 5.0
  ratio:
    units: dimensionless
    value: 0.75
  type: tukey

Superposition 1D¶

This distribution allows the user to superimpose an arbitrary number of 1D PDFs. The general form is:

$\rho(x) = \sum_i w_i \frac{\rho_i(x)}{\max(\rho_i)}$.

Here the $w_i$ are user specified weights. If no weight is specified for a given $\rho_i$, then the weight(s) will default to one.

Example 1: equal weighting

In [38]:

gen = Generator("data/multi.gaussian.in.yaml", verbose=0)
print("\nYAML Input:")
print(yaml.dump({"z_dist": gen.input["z_dist"]}))
plot_dist1d(gen.beam(), "z", units["z"], title_on=True)
print(gen["z_dist:weights"])

gen = Generator("data/multi.gaussian.in.yaml", verbose=0) print("\nYAML Input:") print(yaml.dump({"z_dist": gen.input["z_dist"]})) plot_dist1d(gen.beam(), "z", units["z"], title_on=True) print(gen["z_dist:weights"])

YAML Input:
z_dist:
  dists:
    d1:
      avg_z:
        units: millimeter
        value: -2
      sigma_z:
        units: millimeter
        value: 1.0
      type: gaussian
    d2:
      avg_z:
        units: millimeter
        value: 2
      sigma_z:
        units: millimeter
        value: 1
      type: gaussian
  type: superposition
  weights:
    d1: 1
    d2: 1

{'d1': 1, 'd2': 1}

Example 2: with user specified weighting:

In [39]:

gen["z_dist:weights"] = {"d1": 2, "d2": 1}
print("\nYAML Input:")
print(yaml.dump({"z_dist": gen.input["z_dist"]}))
plot_dist1d(gen.beam(), "z", units["z"], title_on=True)

gen["z_dist:weights"] = {"d1": 2, "d2": 1} print("\nYAML Input:") print(yaml.dump({"z_dist": gen.input["z_dist"]})) plot_dist1d(gen.beam(), "z", units["z"], title_on=True)

YAML Input:
z_dist:
  dists:
    d1:
      avg_z:
        units: millimeter
        value: -2
      sigma_z:
        units: millimeter
        value: 1.0
      type: gaussian
    d2:
      avg_z:
        units: millimeter
        value: 2
      sigma_z:
        units: millimeter
        value: 1
      type: gaussian
  type: superposition
  weights:
    d1: 2
    d2: 1

Product 1D¶

In [40]:

gen = Generator("data/product.gaussian.in.yaml", verbose=0)
print("\nYAML Input:")
print(yaml.dump({"z_dist": gen.input["z_dist"]}))
plot_dist1d(gen.beam(), "z", units["z"], title_on=True)

gen = Generator("data/product.gaussian.in.yaml", verbose=0) print("\nYAML Input:") print(yaml.dump({"z_dist": gen.input["z_dist"]})) plot_dist1d(gen.beam(), "z", units["z"], title_on=True)

YAML Input:
z_dist:
  dists:
    d1:
      max_z:
        units: millimeter
        value: 3
      min_z:
        units: millimeter
        value: -3
      type: uniform
    d2:
      avg_z:
        units: millimeter
        value: 0
      sigma_z:
        units: millimeter
        value: 2
      type: gaussian
  type: product

Deformable¶

In [41]:

gen = Generator("data/deformable.in.yaml", verbose=0)
print("\nYAML Input:")
print(yaml.dump({"t_dist": gen.input["t_dist"]}))
plot_current_profile(gen.beam(), units["t"], units["I"], title_on=True)

gen = Generator("data/deformable.in.yaml", verbose=0) print("\nYAML Input:") print(yaml.dump({"t_dist": gen.input["t_dist"]})) plot_current_profile(gen.beam(), units["t"], units["I"], title_on=True)

YAML Input:
t_dist:
  alpha:
    units: dimensionless
    value: 0.25
  avg_t:
    units: picosecond
    value: 0
  sigma_t:
    units: picosecond
    value: 1
  slope_fraction:
    units: dimensionless
    value: 1
  type: deformable

Interpolation¶

For the 1d interpolation distribution the user provides N points of the PDF for the desired variable and interpolation is employed to make a smooth curve. The PDF knots can be passed in as a dictionary or as an array. Passing in a corresponding coordinate vector is optional. If a coordinate vector is not input, then the user must use a transform operator to scale the resulting distribution (shown below).

In [42]:

gen = Generator("data/interp_ex.yaml")
print("\nYAML Input:")
print(yaml.dump({"t_dist": gen.input["t_dist"]}))

gen = Generator("data/interp_ex.yaml") print("\nYAML Input:") print(yaml.dump({"t_dist": gen.input["t_dist"]}))

YAML Input:
t_dist:
  Pt:
    pt01: 0.1
    pt02: 0.3
    pt03: 4.0
    pt04: 2.9
    pt05: 8.5
    pt06: 2.0
    pt07: 0.1
    pt08: 1.0
    pt09: 5.0
    pt10: 10.0
    pt11: 6.0
    pt12: 1.0
  avg_t:
    units: nanosecond
    value: 0.0
  method: spline
  sigma_t:
    units: nanosecond
    value: 0.08287285705
  type: interp

In [43]:

ax = plot_dist1d(gen.beam(), "t", units["t"], title_on=True)

ax = plot_dist1d(gen.beam(), "t", units["t"], title_on=True)

Special Temporal Distributions¶

---

Pulse stacking temporal crystals¶

In [44]:

gen = Generator("data/t.crystals.in.yaml", verbose=0)
print("\nYAML Input:")
print(yaml.dump({"t_dist": gen.input["t_dist"]}))
plot_current_profile(gen.beam(), units["t"], units["I"], title_on=True)

gen = Generator("data/t.crystals.in.yaml", verbose=0) print("\nYAML Input:") print(yaml.dump({"t_dist": gen.input["t_dist"]})) plot_current_profile(gen.beam(), units["t"], units["I"], title_on=True)

YAML Input:
t_dist:
  crystal_angle_1:
    units: degree
    value: 0.6
  crystal_angle_2:
    units: degree
    value: 1.8
  crystal_angle_3:
    units: degree
    value: -0.9
  crystal_angle_4:
    units: degree
    value: -0.5
  crystal_length_1:
    units: millimeter
    value: 15.096
  crystal_length_2:
    units: millimeter
    value: 7.548
  crystal_length_3:
    units: millimeter
    value: 3.774
  crystal_length_4:
    units: millimeter
    value: 1.887
  type: crystals

$sech^2$ Pulse Laser Soliton¶

Similar to Gaussian in shape, but more accurate for certain types of pulsed laser shapes. The example below shows setting up the distribution using the pulse duraction $\tau$. It is possible to also use $\sigma_t$ instead.

In [45]:

gen = Generator("data/t-sech2.in.yaml", verbose=0)
print("\nYAML Input:")
print(yaml.dump({"t_dist": gen.input["t_dist"]}))
plot_current_profile(gen.beam(), units["t"], units["I"], title_on=True)

gen = Generator("data/t-sech2.in.yaml", verbose=0) print("\nYAML Input:") print(yaml.dump({"t_dist": gen.input["t_dist"]})) plot_current_profile(gen.beam(), units["t"], units["I"], title_on=True)

YAML Input:
t_dist:
  n_tau_cutoff: 4
  tau:
    units: picosecond
    value: 3
  type: sech2

In [ ]: