Distgen Generator¶

Particle generation

In [1]:

# Useful for debugging
%load_ext autoreload
%autoreload 2

# Useful for debugging %load_ext autoreload %autoreload 2

In [2]:

from distgen import Generator
from astra import Astra
from pmd_beamphysics import ParticleGroup
from pmd_beamphysics.plot import marginal_plot
import os

from distgen import Generator from astra import Astra from pmd_beamphysics import ParticleGroup from pmd_beamphysics.plot import marginal_plot import os

In [3]:

# Get an input file
distgen_in = 'templates/dcgun/distgen.yaml'

# Get an input file distgen_in = 'templates/dcgun/distgen.yaml'

In [4]:

# Make object, and show its input (calls __repr__)
D = Generator(distgen_in, verbose=True)

# Change something
D.input['n_particle'] = 10000

D

# Make object, and show its input (calls __repr__) D = Generator(distgen_in, verbose=True) # Change something D.input['n_particle'] = 10000 D

Out[4]:

<disgten.Generator with input: 
n_particle: 100000
output:
  type: null
r_dist:
  max_r:
    units: millimeter
    value: 0.5
  type: radial_uniform
random_type: hammersley
start:
  MTE:
    units: millielectron_volt
    value: 250.0
  type: cathode
t_dist:
  avg_t:
    units: picosecond
    value: 0.0
  n_sigma_cutoff: 3
  sigma_t:
    units: picosecond
    value: 8.5
  type: gaussian
total_charge:
  units: picocoulomb
  value: 100.0
transforms: null

>

In [5]:

# Run, and get particles
D.run()
P = D.particles
P.status

# Run, and get particles D.run() P = D.particles P.status

Distribution format: None
   Warning: no output file specified, defaulting to "None".
Output file: None

Creating beam distribution....
   Beam starting from: cathode
   Total charge: 100 pC.
   Number of macroparticles: 100000.
   Assuming cylindrical symmetry...
   r distribution: radial uniform
      min_r = 0 mm, max_r = 0.5 mm
   theta distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   t distribution: Gaussian
      avg_t = 0 ps, sigma_t = 8.500 ps
      Left n_sigma_cutoff = 3, Right n_sigma_cutoff = -3
   px distribution: Gaussian
      avg_px = 0 eV/c, sigma_px = 357.421 eV/c
   py distribution: Gaussian
      avg_py = 0 eV/c, sigma_py = 357.421 eV/c
   pz distribution: Gaussian
      avg_pz = 0 eV/c, sigma_pz = 357.421 eV/c
   Shifting avg_x = -9.79286E-06 mm -> 0 mm
   Scaling sigma_x = 0.249992 mm -> 0.25 mm
   Shifting avg_y = -9.80995E-07 mm -> 0 mm
   Scaling sigma_y = 0.249999 mm -> 0.25 mm
   Shifting avg_px = -0.0489346 eV/c -> 0 eV/c
   Scaling sigma_px = 357.414 eV/c -> 357.421 eV/c
   Shifting avg_py = -0.0564188 eV/c -> 0 eV/c
   Scaling sigma_py = 357.379 eV/c -> 357.421 eV/c
   Shifting avg_pz = -0.0878522 eV/c -> 0 eV/c
   Scaling sigma_pz = 357.382 eV/c -> 357.421 eV/c
   Shifting avg_t = -0.00102613 ps -> 0 ps
   Scaling sigma_t = 8.38575 ps -> 8.38592 ps
   Cathode start: fixing pz momenta to forward hemisphere
      avg_pz -> 285.197 eV/c, sigma_pz -> 215.435 eV/c
...done. Time Elapsed: 373.021 ms.

   Created particles in .particles: 
   ParticleGroup with 100000 particles with total charge 1.0000000000000003e-10 C

Out[5]:

array([0, 0, 0, ..., 0, 0, 0])

In [6]:

# Check stats from ParticleGroup's calculation
for x in ['x', 'y', 't', 'px', 'py', 'pz']:
    k = 'sigma_'+x
    print(k, ':', P[k], P.units(k))

# Check stats from ParticleGroup's calculation for x in ['x', 'y', 't', 'px', 'py', 'pz']: k = 'sigma_'+x print(k, ':', P[k], P.units(k))

sigma_x : 0.00025 m
sigma_y : 0.00025000000000000006 m
sigma_t : 8.385916336743921e-12 s
sigma_px : 357.42095279935677 eV/c
sigma_py : 357.4209527993567 eV/c
sigma_pz : 215.43483657497615 eV/c

In [7]:

for x in ['x', 'y', 't', 'px', 'py', 'pz']:
    k1 = 'min_'+x
    k2 = 'max_'+x
    print(k1, ':', P[k1], P.units(k1))
    print(k2, ':', P[k2], P.units(k2))

for x in ['x', 'y', 't', 'px', 'py', 'pz']: k1 = 'min_'+x k2 = 'max_'+x print(k1, ':', P[k1], P.units(k1)) print(k2, ':', P[k2], P.units(k2))

min_x : -0.0004997475302364819 m
max_x : 0.0004996360981421189 m
min_y : -0.0004997069647996915 m
max_y : 0.0004997923268076736 m
min_t : -2.548869924080168e-11 s
max_t : 2.5469319141097204e-11 s
min_px : -1629.753693583419 eV/c
max_px : 1504.621099189629 eV/c
min_py : -1537.404841982019 eV/c
max_py : 1448.1372913221817 eV/c
min_pz : 0.0016263845082297446 eV/c
max_pz : 1562.0505583752838 eV/c

Run Astra with Distgen¶

In [8]:

# Input template file 
ASTRA_IN = 'templates/dcgun/astra.in'

# Input template file ASTRA_IN = 'templates/dcgun/astra.in'

In [9]:

# Make an Astra object
A = Astra(input_file=ASTRA_IN, initial_particles=P)

# Change some inputs
A.input['newrun']['zemit']  = 1000
A.input['newrun']['zphase'] = 20
A.input['newrun']['phases'] = True
A.input['newrun']['zstop']  = 1

# Special flag
A.verbose = False

# Make an Astra object A = Astra(input_file=ASTRA_IN, initial_particles=P) # Change some inputs A.input['newrun']['zemit'] = 1000 A.input['newrun']['zphase'] = 20 A.input['newrun']['phases'] = True A.input['newrun']['zstop'] = 1 # Special flag A.verbose = False

In [10]:

A.run()

A.run()

Plot¶

In [11]:

import matplotlib
import matplotlib.pyplot as plt

%matplotlib inline
%config InlineBackend.figure_format = 'retina'
matplotlib.rcParams['figure.figsize'] = (8,5)

import matplotlib import matplotlib.pyplot as plt %matplotlib inline %config InlineBackend.figure_format = 'retina' matplotlib.rcParams['figure.figsize'] = (8,5)

In [12]:

plt.plot(A.stat('mean_z'), A.stat('sigma_x'))
plt.scatter(A.particle_stat('mean_z'), A.particle_stat('sigma_x'))

plt.plot(A.stat('mean_z'), A.stat('sigma_x')) plt.scatter(A.particle_stat('mean_z'), A.particle_stat('sigma_x'))

Out[12]:

<matplotlib.collections.PathCollection at 0x1520c4fa0>

Compare with AstraGenerator¶

In [13]:

from astra import AstraGenerator

from astra import AstraGenerator

In [14]:

GENERATOR_IN = 'templates/dcgun/generator.in'
# Make generator object
G = AstraGenerator(input_file=GENERATOR_IN)
G.input['ipart'] = 10000
G.run()
P2 = G.output['particles']

GENERATOR_IN = 'templates/dcgun/generator.in' # Make generator object G = AstraGenerator(input_file=GENERATOR_IN) G.input['ipart'] = 10000 G.run() P2 = G.output['particles']

In [15]:

A2 = Astra(input_file=ASTRA_IN, initial_particles=P2)

# Change some inputs
A2.input['newrun']['zemit']  = 1000
A2.input['newrun']['zphase'] = 20
A2.input['newrun']['phases'] = True
A2.input['newrun']['zstop']  = 1
A2.run()

A2 = Astra(input_file=ASTRA_IN, initial_particles=P2) # Change some inputs A2.input['newrun']['zemit'] = 1000 A2.input['newrun']['zphase'] = 20 A2.input['newrun']['phases'] = True A2.input['newrun']['zstop'] = 1 A2.run()

In [16]:

k1 = 'mean_z'
k2 = 'norm_emit_x'

u1 = A.units(k1)
u2 = A.units(k2)
plt.xlabel(k1+f' ({u1})')
plt.ylabel(k2+f' ({u2})')
plt.yscale('log')
plt.plot(A.stat(k1), A.stat(k2), color='red', label='distgen')
plt.plot(A2.stat(k1), A2.stat(k2), color='green', label='generator')
plt.scatter(A.particle_stat(k1), A.particle_stat(k2), color='red', label='distgen'),
plt.scatter(A2.particle_stat(k1), A2.particle_stat(k2), color='green', marker='x', label='generator')
plt.legend()

k1 = 'mean_z' k2 = 'norm_emit_x' u1 = A.units(k1) u2 = A.units(k2) plt.xlabel(k1+f' ({u1})') plt.ylabel(k2+f' ({u2})') plt.yscale('log') plt.plot(A.stat(k1), A.stat(k2), color='red', label='distgen') plt.plot(A2.stat(k1), A2.stat(k2), color='green', label='generator') plt.scatter(A.particle_stat(k1), A.particle_stat(k2), color='red', label='distgen'), plt.scatter(A2.particle_stat(k1), A2.particle_stat(k2), color='green', marker='x', label='generator') plt.legend()

Out[16]:

<matplotlib.legend.Legend at 0x152203fd0>

In [17]:

# distgen
A.initial_particles.plot('x', 'y', bins=80)

# distgen A.initial_particles.plot('x', 'y', bins=80)

In [18]:

# Generator
A2.initial_particles.plot('x', 'y', bins=80)

# Generator A2.initial_particles.plot('x', 'y', bins=80)

In [19]:

# Distgen initial particles in time, z-momentum
A.initial_particles.plot( 't', 'pz', bins=200)

# Distgen initial particles in time, z-momentum A.initial_particles.plot( 't', 'pz', bins=200)

In [20]:

# Generator initial particles in time, z-momentum
# The spike at pz=0 is a known limitation
A2.initial_particles.plot('t', 'pz', bins=200)

# Generator initial particles in time, z-momentum # The spike at pz=0 is a known limitation A2.initial_particles.plot('t', 'pz', bins=200)