PYTHIA8.jl Tutorial
This turorial demostrates a simple application of PYTHIA8.jl to generate some events and present the results in form of histograms. Please refer to the PYTHIA documentation for more details on the parameters and the physics processes.
You can also download this tutorial as a Jupyter notebook and a plain Julia source file.
Table of contents
- Generate events and plot the charged multiplicity distribution
- Multi-threaded version of the same example
- Speed comparison between the two versions
Generate events and plot the charged multiplicity distribution
In this example, we will generate 200 events at a center of mass energy of 8000 GeV and plot the charged multiplicity distribution.
Loading the necessary modules
- We will use the
PYTHIA8module to run the event generation. - We will use the
FHistmodule to create the histograms. - We will use the
Plotsmodule to plot them.
If these modules are not installed, you can install them by running the following commands:
using Pkg
Pkg.add("PYTHIA8")
Pkg.add("FHist")
Pkg.add("Plots")using PYTHIA8
using FHist
using Plots: plot, plot!, themeDefine the theme for the plots. See available themes and attributes in Plots.jl.
theme(:wong2, frame=:box, grid=false, minorticks=true,
guidefontvalign=:top, guidefonthalign=:right,
xlims=(:auto, :auto), lw=1.2, lab="", colorbar=false)Initialize PYTHIA
We will initialize PYTHIA with the following settings:
- Center of mass energy of the collision: 8000 GeV
- All hard QCD processes are enabled
- Minimum transverse momentum of the hard process: 20 GeV
We use the operator << to set the parameters and the operator |> to pipe to the init function. This added interface is quite ergonomic, however alternatively, we can use the functions directly provided by C++ interface, e.g.,
readString(pythia, "Beams:eCM = 8000.")orparm( settings(pythia), "Beams:eCM", 8000.)
pythia = PYTHIA8.Pythia("", false);
pythia << "Beams:eCM = 8000." <<
"HardQCD:all = on" <<
"PhaseSpace:pTHatMin = 20.";The purpose of the next two lines is to reduce the amount of output during the event generation
pythia << "Next:numberShowEvent = 0" <<
"Next:numberShowProcess = 0";Initialize the event generation. The init function returns a boolean value. Alternatively, we can use the direct call to init(pythia) function.
pythia |> init;
*------- PYTHIA Process Initialization --------------------------*
| |
| We collide p+ with p+ at a CM energy of 8.000e+03 GeV |
| |
|------------------------------------------------------------------|
| | |
| Subprocess Code | Estimated |
| | max (mb) |
| | |
|------------------------------------------------------------------|
| | |
| g g -> g g 111 | 1.404e+00 |
| g g -> q qbar (uds) 112 | 1.817e-02 |
| q g -> q g 113 | 1.011e+00 |
| q q(bar)' -> q q(bar)' 114 | 1.100e-01 |
| q qbar -> g g 115 | 8.378e-04 |
| q qbar -> q' qbar' (uds) 116 | 3.698e-04 |
| g g -> c cbar 121 | 5.989e-03 |
| q qbar -> c cbar 122 | 1.225e-04 |
| g g -> b bbar 123 | 5.398e-03 |
| q qbar -> b bbar 124 | 1.160e-04 |
| |
*------- End PYTHIA Process Initialization -----------------------*
*------- PYTHIA Multiparton Interactions Initialization ---------*
| |
| sigmaNonDiffractive = 51.78 mb |
| |
| pT0 = 2.35 gives sigmaInteraction = 258.98 mb: accepted |
| |
*------- End PYTHIA Multiparton Interactions Initialization -----*
*------- PYTHIA Flag + Mode + Parm + Word + FVec + MVec + PVec + WVec Settings (changes only) ------------------*
| |
| Name | Now | Default Min Max |
| | | |
| Beams:eCM | 8000.000 | 14000.000 0.0 |
| HardQCD:all | on | off |
| Next:numberShowEvent | 0 | 1 0 |
| Next:numberShowProcess | 0 | 1 0 |
| PhaseSpace:pTHatMin | 20.00000 | 0.0 0.0 |
| |
*------- End PYTHIA Flag + Mode + Parm + Word + FVec + MVec + PVec + WVec Settings -----------------------------*
-------- PYTHIA Particle Data Table (changed only) ------------------------------------------------------------------------------
id name antiName spn chg col m0 mWidth mMin mMax tau0 res dec ext vis wid
no onMode bRatio meMode products
no particle data has been changed from its default value
-------- End PYTHIA Particle Data Table -----------------------------------------------------------------------------------------Create a 1D histogram
We will create a histogram to store the charged multiplicity distribution. The histogram is defined with 25 bins ranging from -0.5 to 499.5. Note that with ranges we need to specify the number of edges (26) and not the number of bins (25)
mult = Hist1D(binedges=range(-0.5, 499.5, 26));Generate events
We will generate 200 events and fill the histogram with the number of charged particles making use of the count function. We use the push! function to fill the histogram.
for iEvent in 1:200
pythia |> next || continue
# Find number of all final charged particles and fill histogram.
nCharged = count(p -> isFinal(p) && isCharged(p), pythia |> event)
push!(mult, nCharged)
end
-------- PYTHIA Info Listing ----------------------------------------
Beam A: id = 2212, pz = 4.000e+03, e = 4.000e+03, m = 9.383e-01.
Beam B: id = 2212, pz = -4.000e+03, e = 4.000e+03, m = 9.383e-01.
In 1: id = 21, x = 1.534e-02, pdf = 5.499e+00 at Q2 = 4.987e+02.
In 2: id = 21, x = 2.093e-03, pdf = 1.839e+01 at same Q2.
Subprocess g g -> g g with code 111 is 2 -> 2.
It has sHat = 2.054e+03, tHat = -1.202e+03, uHat = -8.525e+02,
pTHat = 2.233e+01, m3Hat = 0.000e+00, m4Hat = 0.000e+00,
thetaHat = 1.742e+00, phiHat = 3.620e-01.
alphaEM = 7.695e-03, alphaS = 1.673e-01 at Q2 = 4.987e+02.
Impact parameter b = 5.084e-01 gives enhancement factor = 2.369e+00.
Max pT scale for MPI = 2.233e+01, ISR = 2.233e+01, FSR = 2.233e+01.
Number of MPI = 10, ISR = 31, FSRproc = 82, FSRreson = 0.
-------- End PYTHIA Info Listing ------------------------------------Print the statistics and plot charged multiplicity distribution
Print the generation statistics using the stat function from PYTHIA.
pythia |> PYTHIA8.stat
*------- PYTHIA Event and Cross Section Statistics -------------------------------------------------------------*
| |
| Subprocess Code | Number of events | sigma +- delta |
| | Tried Selected Accepted | (estimated) (mb) |
| | | |
|-----------------------------------------------------------------------------------------------------------------|
| | | |
| g g -> g g 111 | 802 112 112 | 2.007e-01 1.035e-02 |
| g g -> q qbar (uds) 112 | 10 1 1 | 3.338e-03 3.338e-03 |
| q g -> q g 113 | 598 77 77 | 1.372e-01 8.658e-03 |
| q q(bar)' -> q q(bar)' 114 | 63 7 7 | 1.614e-02 2.981e-03 |
| q qbar -> g g 115 | 2 1 1 | 4.818e-04 4.818e-04 |
| q qbar -> q' qbar' (uds) 116 | 1 1 1 | 3.409e-04 3.409e-04 |
| g g -> c cbar 121 | 3 0 0 | 0.000e+00 0.000e+00 |
| q qbar -> c cbar 122 | 0 0 0 | 0.000e+00 0.000e+00 |
| g g -> b bbar 123 | 3 1 1 | 8.686e-04 8.686e-04 |
| q qbar -> b bbar 124 | 0 0 0 | 0.000e+00 0.000e+00 |
| | | |
| sum | 1482 200 200 | 3.592e-01 1.426e-02 |
| |
*------- End PYTHIA Event and Cross Section Statistics ----------------------------------------------------------*
*------- PYTHIA Error and Warning Messages Statistics ----------------------------------------------------------*
| |
| times message |
| |
| 0 no errors or warnings to report |
| |
*------- End PYTHIA Error and Warning Messages Statistics ------------------------------------------------------*Plot the histogram with the plot function from Plots.jl. We use the stepbins seriestype to plot the histogram as a step plot.
img = plot(mult, xlabel="counts", ylabel="#particles",
title="Charged Multiplicity Distribution", seriestype=:stepbins,
c=2, lc=1, fill=0);
display(img)
Multi-threaded version of the same example
In this example we repeat the same example as before but using multi-threading to generate the events faster.
Initialize Parallel PYTHIA
Create an instance of PythiaParallel with the same settings as before. The PythiaParallel class is used to generate events concurrently using multiple threads. The settings are copied for each Pythia instance.
pythia_mt = PYTHIA8.PythiaParallel("", false);PythiaParallel reads settings the same way as the normal Pythia does.
pythia_mt << "Beams:eCM = 8000." <<
"HardQCD:all = on" <<
"PhaseSpace:pTHatMin = 20.";The maximum degree of parallelism. If set to 0 (default), the program will use the maximum number of threads supported by the hardware.
pythia_mt << "Parallelism:numThreads = 4";This defines the number of events generated by PythiaParallel::run.
pythia_mt << "Main:numberOfEvents = 200";The next is a user provided function to initialize each underlying Pythia instance (one per thread). Please note that the function is called concurrently, therefore the use needs to avoid calling thread-unsafe functions. This is the reason we use the Core.println function instead of println.
function w_init(pythia)::CxxBool
Core.println("Initializing Pythia with index $(mode(pythia |> settings, "Parallelism:index")).")
return pythia |> init
end;Initialize the event generation. The init function returns a boolean value to indicate the success. The user provided function w_init is called for each thread.
init(pythia_mt, w_init);Initializing Pythia with index 1.Initializing Pythia with index 3.Initializing Pythia with index 0.Initializing Pythia with index 2.Create a 1D histogram
mult_mt = Hist1D(binedges=range(-0.5, 499.5, 26));Generate events
Define a user function that will be called for each generated event. This function will be called concurrently by the threads. Therefore, we need to use thread-safe functions. Filling the histogram with atomic_push is thread-safe. The function takes a Pythia object as the first argument and returns Nothing.
function analyze(pythiaNow)::Nothing
nCharged = count(p -> isFinal(p) && isCharged(p), pythiaNow |> event)
atomic_push!(mult_mt, nCharged)
return
endanalyze (generic function with 1 method)Run the event generation. The number of events has been set in the Main:numberOfEvents parameter.
PYTHIA8.run(pythia_mt, analyze);Print the statistics and plot charged multiplicity distribution
Print the statistics using the stat function from PYTHIA8.jl
pythia_mt |> PYTHIA8.stat PYTHIA Warning in PythiaParallel::init: experimental feature, please send feedback to authors@pythia.org
*------- PYTHIA Error and Warning Messages Statistics ----------------------------------------------------------*
| |
| times message |
| |
| 4 Warning in Pythia::init: be aware that successive calls to init() do not clear previous settings |
| 1 Warning in PythiaParallel::init: experimental feature, please send feedback to authors@pythia.org |
| |
*------- End PYTHIA Error and Warning Messages Statistics ------------------------------------------------------*Plot the result histogram
img = plot(mult_mt, xlabel="counts", ylabel="#particles",
title="Charged Multiplicity Distribution (multi-threaded)", seriestype=:stepbins,
c=2, lc=1, fill=0);
display(img)
Speed comparison between the two versions
We will generate 1000 events and compare the time taken by the two versions.
empty!(mult);
empty!(mult_mt);Single-threaded version
Generate 1000 events and measure the time taken. This is done using the @elapsed macro.
st_elap = @elapsed begin
for iEvent in 1:1000
pythia |> next || continue
nCharged = count(p -> isFinal(p) && isCharged(p), pythia |> event)
atomic_push!(mult, nCharged)
end
end4.898925607Multi-threaded version
Generate 1000 events and measure the time taken. This is done using the @elapsed macro.
pythia_mt << "Main:numberOfEvents = 1000"
mt_elap = @elapsed begin
PYTHIA8.run(pythia_mt, analyze)
end2.002769466Print the speedup factor (4 threads)
println("Speedup is $(st_elap/mt_elap).")Speedup is 2.44607564183825.Make a plot of the speedup factor
We will generate events using 1, 2, 4, 8, 16, and 32 threads and plot the speedup factor. The speedup factor is defined as the ratio of the time taken by the single-threaded version to the time taken by the multi-threaded version.
n_threads = [2^n for n in 0:5]
speedup = Float64[]
w_init(pythiaNow)::CxxBool = pythiaNow |> init
for n_thread in n_threads
empty!(mult_mt)
_pythia = PYTHIA8.PythiaParallel("", false);
_pythia << "Beams:eCM = 8000." <<
"HardQCD:all = on" <<
"PhaseSpace:pTHatMin = 20." <<
"Parallelism:numThreads = $n_thread" <<
"Main:numberOfEvents = 1000"
init(_pythia, w_init)
elap = @elapsed begin
PYTHIA8.run(_pythia, analyze)
end
push!(speedup, st_elap/elap)
end;Plot the speedup factor vs the number of threads and overlay the line for the ideal scaling.
img = plot(n_threads, speedup, xlabel="#threads", ylabel="speedup",
title="Speedup vs #threads (#cores = $(Sys.CPU_THREADS))",
seriestype=:scatter, legend=false,
xscale=:log10, yscale=:log10);Add a line to the same plot
plot!(n_threads, n_threads, seriestype=:line)
display(img)
The plot is generated in a CI node using the cores available in the node. The speedup factor may not scale as expected.
This page was generated using Literate.jl.