MCBooster: a library for fast Monte Carlo generation of phase-space decays on massively parallel platforms

TL;DR

MCBooster is a high-performance C++ library that leverages parallel computing on CUDA and OpenMP platforms to efficiently generate large samples of phase-space Monte Carlo events for physics simulations.

Contribution

It introduces a header-only, platform-independent library that significantly accelerates Monte Carlo phase-space generation using parallel processing techniques.

Findings

Achieves high performance on CUDA and OpenMP devices

Supports large-scale phase-space event generation

Demonstrates efficiency across various computing environments

Abstract

MCBooster is a header-only, C++11-compliant library that provides routines to generate and perform calculations on large samples of phase space Monte Carlo events. To achieve superior performance, MCBooster is capable to perform most of its calculations in parallel using CUDA- and OpenMP-enabled devices. MCBooster is built on top of the Thrust library and runs on Linux systems. This contribution summarizes the main features of MCBooster. A basic description of the user interface and some examples of applications are provided, along with measurements of performance in a variety of environments

Figures

Tables

Table 1. Table 1: The time spent, in seconds, by different NVidia GPU models to process 10 million events. The row labeled B 0 → K − ​ π + ​ J / ψ → superscript 𝐵 0 superscript 𝐾 superscript 𝜋 𝐽 𝜓 B^{0}\to K^{-}\pi^{+}J/\psi is the time taken on the device to generate this three-body decay. The row labeled J / ψ → μ + ​ μ − → 𝐽 𝜓 superscript 𝜇 superscript 𝜇 J/\psi\to\mu^{+}\mu^{-} is the time taken on the device to generate this two-body decay. The row labeled Dataset is the time taken to do the calculations illustrated in 1 .

	Time spent per task for different cards (s)
	K2200	GTX 970	GTX TITAN Z	K40c
$B^0\to K^{-}{\pi }^{+}J/\psi$	0.034	0.016	0.020	0.021
$J/\psi \to {\mu }^{+}{\mu }^{-}$	$12\times {10}^{-06}$	$8.8\times {10}^{-06}$	$9.5\times {10}^{-06}$	$14\times {10}^{-06}$
Dataset	0.18	0.10	0.078	0.049

Table 2. Table 2: The time spent, in seconds, to generate 10 million B 0 → K − ​ π + ​ J / ψ → superscript 𝐵 0 superscript 𝐾 superscript 𝜋 𝐽 𝜓 B^{0}\to K^{-}\pi^{+}J/\psi events under OpenMP using a 24 physical-core, 48 logical-core Xeon ES-2680 @2.5 GHz. ROOT :: TGenPhaseSpace runs in a single thread.

	Time spent per number of threads (s)
	#1	#12	#24	#48
$B^0\to K{\pi }^{+}J/\psi$	3.65	0.369	0.218	0.161
$J/\psi \to {\mu }^{+}{\mu }^{-}$	2.39	0.232	0.152	0.136
Dataset	1.69	0.33	0.33	0.271
$B^0\to K{\pi }^{+}J/\psi$ w/ TGenPhaseSpace	4.68	-	-	-

Full text

1 Introduction

Phase space Monte Carlo simulates the decay kinematics of a mother particle decaying into $N$ daughter particles with no intermediate resonances. Samples of phase space Monte Carlo events are widely used in high energy physics for calculating phase volumes. They are also used as the starting point to implement and describe decays containing one or more intermediate resonances and to simulate the response of detectors [1].

MCBooster is a header-only, C++11-compliant library for the generation of large samples of phase-space Monte Carlo events on massively parallel platforms. It was released on GitHub 111https://github.com/MultithreadCorner/MCBooster in the spring of 2016. The core libraries implement the Raubold-Lynch algorithm [1]; they are able to generate the full kinematics of decays with up to nine particles in the final state. MCBooster supports the generation of sequential decays as well as the parallel evaluation of arbitrary functions over the generated events. The output of MCBooster accords with popular and well-tested software packages such as GENBOD (W515 from CERNLIB) [2] and TGenPhaseSpace in the ROOT framework [3]. MCBooster is built on top of the Thrust library and runs on Linux systems. It deploys transparently on NVidia CUDA-enabled GPUs as well as multicore CPUs.

This contribution summarizes the main features of MCBooster. A basic description of the user interface is presented in section 3. Some examples of applications are provided in section 6. Measurements of performance on a variety of platforms and in a variety of environments are discussed in section 5.

2 Phase space Monte Carlo

Considering conservation of energy and momentum, and the relativistic mass-energy relation ($m^{2}=E^{2}-p^{2}$), the kinematics of the decay of a mother particle at rest to an N-particle final state can be described by $4N-4-3-N=3N-7$ independent parameters ($4N$ for the four-momenta of the $N$ daughters, 4 for conservation of energy and momentum, 3 for angles in the mother’s center-of-momentum about which the daughters’ momenta can be rotated without changing the kinematics, and $N$ for the mass-energy relations of the daughters).

The Raubold-Lynch algorithm is described in detail on [1]. Basically, the method models an N-body decay as a recursive series of two-body decays [1, 4]. The four-momenta of the final state particles are randomly generated and are kinematically consistent. All events are accepted, and a weight proportional to the local phase space density is calculated in order to ensure the correct event distribution. The resulting sample can be unweighted in parallel using one of the methods provided by the package.

3 Implementation highlights

MCBooster is implemented on top of the Thrust library as a header-only, C++11-compliant library and runs on Linux systems. It supports systems compatible with NVidia CUDA and OpenMP. In each call to the principal method, events are generated in large bunches (typically multiples of a million), depending on the memory available on the device.

Some of the main features of MCBooster are

• •

The complete set of four-vectors is available for the final states. MCBooster can output weighted and unweighted samples.

• •

Interfaces exist for intrusive and non-intrusive evaluation of arbitrary functions of the phase space coordinates.

• •

The MCBooster interface also supports the generation of sequential decays.

Three- and four-vectors are described by the classes Vector3R and Vector4R. Event generation is managed by the class PhaseSpace, which is configured via the constructor with the masses of the particles and the number of events to generate. The methods PhaseSpace::Generate(Vector4R mother) and PhaseSpace::Generate(MCBooster::Particles mothers) perform the actual event generation. Generated events are kept on the device. The containers have the lifetime of the corresponding MCBooster::PhaseSpace object. Allocated memory can be reused in subsequent runs. The user can export the generated events from the device to the host and store them in the dedicated container class Events.

MCBooster provides three interfaces to perform the evaluation of arbitrary functions over the generated events in parallel. To use such interfaces, it is necessary only to derive a functor from the corresponding interface and to implement the corresponding operator().

4 A worked example

To exemplify some of the basic functionality of MCBooster, the code snippets below show how to generate 10 million events and calculate some parameters corresponding to the decay chain $B^{0}\to K^{\ast}J/\psi$, with $K^{\ast}\to K^{+}\pi^{-}$ and $J/\psi\to\mu^{+}\mu^{-}$. The parameters are

• •

$M(K,\pi)$, the invariant mass of the $K\pi$ system;

• •

$\cos(\theta_{K})$, the helicity angle of the $K^{\ast}$;

• •

$\cos(\theta_{\mu})$, the helicity angle of the $J/\psi$;

• •

$\Delta\phi$ difference between the decay planes of the $K^{\ast}$ and $J/\psi$.

Given it’s very narrow width, the invariant mass of the $J/\psi$ is fixed at its nominal value, an excellent approximation given its very narrow width, one can use the two-body $K\pi$ invariant mass plus the angular variables to fully specify the kinematic state of each generated event.

4.1 Generation of the events

Typical C++ code to generate and export to the host memory the $B^{0}\to K^{+}\pi^{-}{J/\psi}$ events, with ${J/\psi}\to\mu^{+}\mu^{-}$ is shown in 1. The comments in the listing explain the details about each line. Basically, the class PhaseSpace is instantiated and configured to generate the first level of the decay chain $B^{0}\to K^{+}\pi^{-}{J/\psi}$. The $B^{0}$s are at rest. The $B^{0}\to K^{+}\pi^{-}{J/\psi}$ events are then copied to the host memory. In the next step, a second PhaseSpace object is instantiated and configured to generate the ${J/\psi}\to\mu^{+}\mu^{-}$ decay. The generation is performed invoking the PhaseSpace::Generate method and passing as argument the list of ${J/\psi}$ daughters taken from the $B^{0}\to K^{+}\pi^{-}{J/\psi}$ decay. The generated ${J/\psi}\to\mu^{+}\mu^{-}$ decays are copied to the host memory as well.

4.2 Evaluation of function over the events

Using MCBooster::Evaluate and MCBooster::EvaluateArray routines, it is possible to evaluate arbitrary function objects in parallel, taking events stored in the device memory as parameters. The results can be kept in the device or copied automatically to the host memory depending on the parameters passed to the algorithms. 3 shows how to make a dataset with four parameters, as discussed in section 4. The results of the calculation are stored in the host memory.

5 Performance

To evaluate the performance of MCBooster, the time spent to perform the various operations has been measured running on different NVidia GPUs and on a multicore CPU deploying different number of OpenMP threads.

The performance measurements made with four CUDA-enabled devices with different architectures are summarized in Table 1. Different parameters concerning each NVidia graphics card need to be considered to understand these numbers. The main parameters are listed below

• •

Quadro K2200 (Maxwell/5.0): 640 CUDA Cores @ 1.12 GHz. FP64/FP32 = 1:32

• •

GeForce GTX 970 (Maxwell/5.2): 1664 CUDA Cores @ 1.18 GHz. FP64/FP32 = 1:32

• •

GeForce Titan Z (Kepler/3.5): 2880 CUDA Cores @ 0.88 GHz. FP64/FP32 = 1:3

• •

Tesla K40c (Kepler/3.5): 2880 CUDA Cores @ 0.75 GHz. FP64/FP32 = 1:3

The performance variation as a function of problem size was measured using the Tesla K40c device. Figure 1 shows that the time taken to generate the underlying three-body phase-space decays is essentially independent of the number of events between 500 thousand and 50 million events, and then begins to rise very gently. The (small absolute) overhead setting up the problem completely dominates at lower statistics, and continues to dominate, even at the highest statistics considered. Figure 2 shows that the time to generate a sample grows roughly linearly as a function of the number of particles in the final state. This is the expected behavior as the Raubold-Lynch algorithm adds one recursive step for each additional final state particle.

The performance using the OpenMP back-end was measured using a 24 physical-core, 48 logical-core Xeon ES-2680 @2.5 GHz CPU. The results are summarized in the Table 2. As a further figure of merit to describe the performance using the OpenMP back-end, the time spent to generate a 10 million event sample with 9 particles in final state was also measured:

• •

MCBooster takes 0.74 seconds using 24 OpemMP threads.

• •

TGenPhaseSpace takes 22 seconds.

One sees that MCBooster provides an order of magnitude performance boost on such a system, even though it was optimized for execution on NVidia GPUs.

6 Summary

The basic design, performance and functionality of the header-only, C++11-compliant library MCBooster have been introduced. The basic interfaces are discussed through the working example presented . The performance measurements for running MCBooster using CUDA and OpemMP back-ends are discussed in section 5, and show that MCBooster can be up to 100 times faster than conventional software, depending on the graphics card or number of threads deployed. Since MCBooster is header only, no additional building process needs be done beyond the inclusion of the required headers. One example of integration and usage of MCBooster can be found the on poster [5].

MCBooster is open source. The project is hosted on GitHub at

https://github.com/MultithreadCorner/MCBooster.

The package includes a suite of examples and also a simple application to generate samples and save them in ROOT TTrees.

7 Acknowledgments

This work was performed with support from NSF Award PHY-1414736. NVidia provided K40 GPUs for our use through its University Partnership program. We are pleased to acknowledge valuable technical support from Christoph Hasse and Brad Hittle.