Basic Example

Data import

The first step of data analysis is import data. In this code, data import is always handled through an instance of the DataLoaders class. In this example, we will work with topas data. For a list of all currently available data importers, see here; for instructions on writing new data importers see here

from pathlib import Path
import sys
sys.path.append('../')  # not necessary when the library is installed
from ParticlePhaseSpace import DataLoaders
from ParticlePhaseSpace import PhaseSpace

test_data_loc = Path(r'../tests/test_data/coll_PhaseSpace_xAng_0.00_yAng_0.00_angular_error_0.0.phsp').absolute()
ps_data = DataLoaders.Load_TopasData(test_data_loc)

Once we have a DataLoader instance, we can pass it to ParticlePhaseSpace for analysis:

PS = PhaseSpace(ps_data)

The underlying data is stored in a pandas data frame at ps_data. You shouldn’t often have to (or want to) interact with this directly, but let’s take a quick look at the data so you know what it looks like:

PS.ps_data.head()

	particle type [pdg_code]	x [mm]	y [mm]	z [mm]	weight	particle id	time [ps]	px [MeV/c]	py [MeV/c]	pz [MeV/c]
0	22	-0.971177	-0.817280	396.0	0.01	0	0	-0.008856	-0.010509	-3.184060
1	22	-0.916198	0.839298	396.0	0.01	1	0	-0.003894	0.008249	-1.560114
2	22	-5.258902	11.318420	396.0	0.01	2	0	-0.193061	0.377104	-1.219742
3	22	15.326971	-25.628899	396.0	0.01	3	0	0.102889	-0.362594	-1.414843
4	22	-0.241948	-0.401026	396.0	0.01	4	0	-0.000192	-0.001470	-0.209232

Calculation of additional quantities

The data show here consists of the required columns - that is, the minimum data required to define a phase space. Many other quantaties can be calculated from this data using various fill methods. For example, to add in Kinetic Energy:

PS.fill.kinetic_E()
PS.ps_data.head()

	particle type [pdg_code]	x [mm]	y [mm]	z [mm]	weight	particle id	time [ps]	px [MeV/c]	py [MeV/c]	pz [MeV/c]	rest mass [MeV/c^2]	p_abs [MeV/c]	Ek [MeV]
0	22	-0.971177	-0.817280	396.0	0.01	0	0	-0.008856	-0.010509	-3.184060	0.0	3.184090	3.184090
1	22	-0.916198	0.839298	396.0	0.01	1	0	-0.003894	0.008249	-1.560114	0.0	1.560141	1.560141
2	22	-5.258902	11.318420	396.0	0.01	2	0	-0.193061	0.377104	-1.219742	0.0	1.291220	1.291220
3	22	15.326971	-25.628899	396.0	0.01	3	0	0.102889	-0.362594	-1.414843	0.0	1.464187	1.464187
4	22	-0.241948	-0.401026	396.0	0.01	4	0	-0.000192	-0.001470	-0.209232	0.0	0.209237	0.209237

Notice that while we only requested Ek, we also got rest mass and absolute momentum. This is because these quantities are required to calculate kinetic energy anyway. The code will always check whether a given quantity is already present in the data before recalculating it.

You can view all the available fill methods with this command:

PS.fill.get_methods()

absolute_momentum
beta_and_gamma
direction_cosines
kinetic_E
relativistic_mass
rest_mass
velocity

If you want to access any of the data inside this data frame, there are two ways: the unit specific way, and the unit agnostic way. The unit specific way is of course how would get data out of any DataFrame:

Ek = PS.ps_data['Ek [MeV]']

You can also access the same quantities without reference to the units like this:

Ek = PS.ps_data[PS.columns['Ek']]

Note that the second method will work regardless of what units the data is in.

Reset the phase space to required columns

In some situations, you may wish to remove all ‘derived’ quantities from the data, and have only the required columns. You may wish to do this to minimise the memory footpring, or because you have carried out some operation that potentially invalidated the calculated quantities (more on that below!) In such situations, you can use reset_phase_space:

PS.reset_phase_space()
PS.ps_data.head()

	particle type [pdg_code]	x [mm]	y [mm]	z [mm]	weight	particle id	time [ps]	px [MeV/c]	py [MeV/c]	pz [MeV/c]
0	22	-0.971177	-0.817280	396.0	0.01	0	0	-0.008856	-0.010509	-3.184060
1	22	-0.916198	0.839298	396.0	0.01	1	0	-0.003894	0.008249	-1.560114
2	22	-5.258902	11.318420	396.0	0.01	2	0	-0.193061	0.377104	-1.219742
3	22	15.326971	-25.628899	396.0	0.01	3	0	0.102889	-0.362594	-1.414843
4	22	-0.241948	-0.401026	396.0	0.01	4	0	-0.000192	-0.001470	-0.209232

Basic analytics

Once we have an instance of ParticlePhaseSpace, we can perform some analysis.

PS.print_energy_stats()

===================================================
                 ENERGY STATS
===================================================
total number of particles in phase space:  311489
number of unique particle species:  3
     308280 gammas
        mean energy:  1.91 MeV
        median energy:  1.20 MeV
        Energy spread IQR:  2.03 MeV
        min energy  0.01 MeV
        max energy  10.35 MeV
     2853 electrons
        mean energy:  3.16 MeV
        median energy:  3.39 MeV
        Energy spread IQR:  2.95 MeV
        min energy  0.02 MeV
        max energy  9.44 MeV
     356 positrons
        mean energy:  2.66 MeV
        median energy:  2.39 MeV
        Energy spread IQR:  2.39 MeV
        min energy  0.08 MeV
        max energy  8.46 MeV

You can see that this phase space consists of three different particle species. We can also generate plots for energy, and particle positions:

PS.plot.energy_hist_1D()

PS.plot.position_hist_1D()

PS.plot.particle_positions_scatter_2D(beam_direction='z')

Note that in the plot above, particles with different weights are plotted with different colors. Another way to visualise the particle positions is the plot_beam_intensity method, which will produce an image of accumulated intensity:

PS.plot.particle_positions_hist_2D(xlim=[-10, 10], ylim=[-10,10], grid=True)

This plot is actually a lot more illuminating than the scatter plot in this instance! This phase space was scored at the exit of a novel type of x-ray collimator; you can see that the x-rays (gammas) have been shaped to a square, which is as expected. Meanwhile the electrons and positrons are randomly scattered, which is also as expected.

as with the fill methods, you can get all plotting methods like this:

PS.plot.get_methods()

energy_hist_1D
momentum_hist_1D
n_particles_v_time
particle_positions_hist_2D
particle_positions_scatter_2D
position_hist_1D
transverse_trace_space_hist_2D
transverse_trace_space_scatter_2D

Seperating, adding, and subtracting phase space objects

You can easily add, subtract, and seperate phase space objects:

electron_PS = PS('electrons')  # this generates a phase space with only the electrons
# you can also use pdg codes:
electron_PS = PS(11)  # this is the same
electron_PS.plot.energy_hist_1D()

You can also pass a list of particles to get multiple phase spaces at once:

gamma_PS, positron_PS = PS(['gammas', 'positrons'])  # you can also pass a list

You can subtract one phase space from another. For example, the following will produce a new phase space object where all the gamma particles are removed:

no_gamma_PS = PS - gamma_PS
no_gamma_PS.plot.energy_hist_1D()

You can also add phase spaces together:

Original_PS = no_gamma_PS + gamma_PS
Original_PS.plot.energy_hist_1D()

However: you cannot add identical particles into an existing phase space. If you really want to do this, you will need to update the particle_ID field first.

Twiss parameters

Twiss parameters or courant-snyder parameters are commonly used to describe the distribution of particles in a transverse phase space. You can calculate, save, and plot the twiss parameters:

PS.print_twiss_parameters(beam_direction='z')  # calculate the RMS twiss parameters

===================================================
                 TWISS PARAMETERS
===================================================

gammas:
                x         y
epsilon  2.525901  2.673361
alpha    0.619585  0.525727
beta     8.976938  9.488084
gamma    0.154160  0.134526

electrons:
                x         y
epsilon  2.208669  2.357446
alpha    1.159332  1.009164
beta     4.971102  4.634599
gamma    0.471535  0.435510

positrons:
                x         y
epsilon  1.881200  1.178092
alpha    1.141023  2.283977
beta     3.888975  5.672388
gamma    0.591913  1.095932

PS.plot.transverse_trace_space_scatter_2D()

The above plot shows what is sometimes called a ‘trace space’ which is plot of position (x) versus divergence (x’ = px/pz). The Red ellipse represents the space enclosed by the RMS twiss parameters, which should enclode 37% of the beam. As an alternative to the scatter plot, you can also visualise the intensity (defined as the sum of weights) in trace space:

PS.plot.transverse_trace_space_hist_2D(xlim=[-2,2], ylim=[-1, 1], plot_twiss_ellipse=True, grid=False)

Again, this is quite illuminating; it shows that the gamma particles are fairly evenly distributed spread between -1 and 1 mm, have almost no divergence (note the highly intesense line in the middle, and the use og a log scale). Again, since this data is from an x-ray collimator, this makes a lot of sense!

Exporting the data

Similarly to the data import stage we use an instance of ParticlePhaseSpace.DataLoaders to export the data. You can see the currently available data exporters here, and you see a demonstration of writing a new data exporter here.

from ParticlePhaseSpace import DataExporters

DataExporters.CSV_Exporter(PS, output_location=Path('.'), output_name='test_export')

<ParticlePhaseSpace.DataExporters.CSV_Exporter at 0x7f943666dca0>