pyFANG

Python implementation of the Focused Angular N-body event Generator (FANG)

pyFANG is a fast Monte Carlo phase space event generator for particle physics. It generates N-body decay events in Lorentz-Invariant Phase Space (LIPS) with optional angular constraints on final-state particles, allowing you to focus event generation into specific detector regions instead of generating over the full solid angle.

This is useful for simulations where you only care about particles hitting a detector at known positions, making your Monte Carlo sampling dramatically more efficient.

If you use pyFANG in your research, please cite the paper. See Citation below.

🚀 Quick Start & Tutorial

The easiest way to learn pyFANG is to walk through our interactive tutorial:

👉 View the Tutorial Notebook (tutorial.ipynb) This notebook provides a step-by-step guide to generating events, setting up detectors, and computing phase space volumes.

What Does pyFANG Do?

In particle physics, when a parent particle decays into N daughter particles, the daughters can fly off in any direction allowed by energy-momentum conservation. Simulating this process is called phase space generation.

Standard generators such as GENBOD (or ROOT’s TGenPhaseSpace) sample the full $4\pi$ center-of-mass solid angle uniformly. However, in most experiments, the fiducial acceptance of the lab-frame detector covers only a small fraction of that sphere. This mismatch leads to significant computational waste, as the vast majority of generated events fail the required fiducial cuts and are discarded. This inefficiency is further compounded in coincidence measurements involving two or more detectors; requiring multiple particles to simultaneously satisfy separate cuts drastically reduces the success rate, making it difficult to accumulate a statistically significant sample of valid events.

pyFANG solves this problem. You define detector objects that describe where your detectors are and what shape they have, then pyFANG generates only events where the particles hit those detectors. The phase space weights are correctly adjusted so that all physical observables remain accurate.

Key idea: Instead of generating 10 million events and keeping 1,000 that hit your detector, pyFANG generates 1,000 events that all hit your detector, each with a correct weight.

Features

Detector objects -- Define your detectors as CircleDetector, StripDetector, RingDetector, or PointDetector and pass them to the generator
Correct weights -- Phase space weights are properly calculated so observables are physically accurate
Fast -- Core algorithms are JIT-compiled with Numba for near-C++ performance
Batch mode -- Generate millions of unconstrained events at once with NumPy vectorization
Multiple solutions -- Automatically finds all kinematically allowed solutions for a given detector configuration
Pure Python -- No C/C++ compilation needed; just pip install and go
Parallel-ready -- Built-in support for multiprocessing across CPU cores

Installation

Requirements

Python 3.8+
NumPy >= 1.20
Numba >= 0.55
Matplotlib >= 3.3 (for visualization)

Install from PyPI

pip install FANG-py

Install from Source

git clone https://github.com/itay-space/FANG.git
cd FANG
pip install .

For development (editable install):

pip install -e .

Quick Start

from pyFANG import TFANG
import numpy as np

# Create a generator
gen = TFANG(seed=42)

# Define the decay:
#   Parent 4-momentum [px, py, pz, mass] -- here a particle at rest with mass 5 GeV
#   Decaying into 3 daughters, each with mass 0.5 GeV
parent = np.array([0.0, 0.0, 0.0, 5.0])
masses = np.array([0.5, 0.5, 0.5])

gen.SetDecay(parent, masses)

# Generate a single event
gen.Generate()

# Get the weight and the daughter 4-momenta
weight = gen.GetWeight(0)
daughters = gen.GetDecays(0)

print(f"Weight: {weight:.6f}")
for i, p in enumerate(daughters):
    energy = np.sqrt(p[0]**2 + p[1]**2 + p[2]**2 + p[3]**2)
    print(f"  Particle {i}: px={p[0]:.4f}, py={p[1]:.4f}, pz={p[2]:.4f}, m={p[3]:.4f}, E={energy:.4f}")

Usage Guide

Basic Phase Space Generation

Generate single unconstrained events one at a time:

from pyFANG import TFANG
import numpy as np

gen = TFANG(seed=42)

# Parent particle: moving along z with momentum 5 GeV/c, mass 12 GeV
parent = np.array([0.0, 0.0, 5.0, 12.0])

# 5-body decay, all daughters with mass 1.0 GeV
masses = np.array([1.0, 1.0, 1.0, 1.0, 1.0])

gen.SetDecay(parent, masses)

# Generate events in a loop
sum_weights = 0.0
n_events = 100000

for _ in range(n_events):
    gen.Generate()
    weight = gen.GetWeight(0)
    sum_weights += weight

    # Access individual daughter momenta
    daughters = gen.GetDecays(0)
    # daughters[i] is [px, py, pz, mass] for particle i

print(f"Average weight: {sum_weights / n_events:.4f}")

Adding Detectors

Constrain daughter particles to point toward specific detectors. This is the core feature of FANG. You create detector objects that describe the shape and direction of each detector, then add them to the generator:

from pyFANG import TFANG, CircleDetector
import numpy as np

gen = TFANG(seed=42)

parent = np.array([0.0, 0.0, 5.0, 12.0])
masses = np.array([1.0, 1.0, 1.0, 1.0, 1.0])
gen.SetDecay(parent, masses)

# Constrain particle 0: must go into a circular cone centered on (0, 0, 1)
# with solid angle 0.5 steradians
gen.AddDetector(CircleDetector(
    direction=np.array([0.0, 0.0, 1.0]),  # detector direction (will be normalized)
    solid_angle=0.5                         # solid angle in steradians
))

# Constrain particle 1: must go into a circular cone centered on (1, 0, 0)
# with solid angle 0.3 steradians
gen.AddDetector(CircleDetector(
    direction=np.array([1.0, 0.0, 0.0]),
    solid_angle=0.3
))

# Particles 2, 3, 4 remain unconstrained (full 4pi)

# Generate and iterate over solutions
gen.Generate()

for i in range(gen.GetNSolutions()):
    weight = gen.GetWeight(i)
    daughters = gen.GetDecays(i)
    # Process this solution...

Detectors are applied to daughter particles in order: the first AddDetector call constrains particle 0, the second constrains particle 1, and so on. Unconstrained particles sample the full 4pi solid angle.

Detector Types

pyFANG provides four detector shapes. Each is a Python class you instantiate and pass to gen.AddDetector():

`CircleDetector` -- Circular cone

The most common type. Particles land anywhere inside a cone of the given solid angle.

from pyFANG import CircleDetector

det = CircleDetector(
    direction=np.array([0.0, 0.0, 1.0]),  # center direction
    solid_angle=0.5                         # cone solid angle in steradians (0 to 4pi)
)
gen.AddDetector(det)

`PointDetector` -- Fixed direction

The particle direction is fixed exactly to the given vector. No solid angle needed. Useful for single-angle differential cross section calculations.

from pyFANG import PointDetector

det = PointDetector(direction=np.array([0.0, 0.0, 1.0]))
gen.AddDetector(det)

`StripDetector` -- Rectangular strip

A rectangular region in (theta, phi) space. azimuthal_coverage is the fraction of 2pi covered in the azimuthal direction (between 0 and 1).

from pyFANG import StripDetector

det = StripDetector(
    direction=np.array([0.0, 1.0, 0.0]),  # center direction
    solid_angle=1.2,                        # solid angle in steradians
    azimuthal_coverage=0.4                  # fraction of 2pi (0 < value <= 1)
)
gen.AddDetector(det)

`RingDetector` -- Annular ring

An annular region around the direction vector, excluding the center of the cone.

from pyFANG import RingDetector

det = RingDetector(
    direction=np.array([0.0, 0.0, 1.0]),  # center direction
    solid_angle=0.5,                        # solid angle in steradians
    opening=0.1                             # ring opening parameter (positive)
)
gen.AddDetector(det)

Mixing detector types

You can use different detector types for different particles in the same event:

from pyFANG import TFANG, CircleDetector, PointDetector, StripDetector
import numpy as np

gen = TFANG(seed=42)
gen.SetDecay(np.array([0.0, 0.0, 5.0, 12.0]), np.array([1.0]*5))

# Particle 0: circular cone detector
gen.AddDetector(CircleDetector(direction=np.array([0.0, 0.0, 1.0]), solid_angle=0.5))

# Particle 1: fixed direction (e.g., beam axis)
gen.AddDetector(PointDetector(direction=np.array([0.0, 0.0, 1.0])))

# Particle 2: strip detector (40% azimuthal coverage)
gen.AddDetector(StripDetector(direction=np.array([0.0, 1.0, 0.0]), solid_angle=1.2, azimuthal_coverage=0.4))

# Particles 3, 4: unconstrained

Batch Generation

GenerateBatch generates many events at once and returns NumPy arrays. It works with or without detectors:

gen = TFANG(seed=42)
parent = np.array([0.0, 0.0, 5.0, 12.0])
masses = np.array([1.0, 1.0, 1.0, 1.0, 1.0])
gen.SetDecay(parent, masses)

# Generate 1 million events at once
momenta, weights = gen.GenerateBatch(1_000_000)

# momenta.shape = (1000000, 5, 4)  -- [event, particle, (px,py,pz,m)]
# weights.shape  = (1000000,)

print(f"Phase space estimate: {weights.mean():.2f} +/- {weights.std() / np.sqrt(len(weights)):.2f}")

It also works with detectors:

gen.AddDetector(CircleDetector(direction=np.array([0.0, 0.0, 1.0]), solid_angle=0.5))
momenta, weights = gen.GenerateBatch(200_000)

Note: Without detectors, GenerateBatch uses a fast vectorized path. With detectors, constrained generation can yield 0 or multiple solutions per call, so the number of returned events may differ from nEvents.

Calculating Phase Space Volumes

pyFANG can estimate the full and partial (detector-constrained) phase space volumes via Monte Carlo integration:

from pyFANG import TFANG, CircleDetector
import numpy as np

gen = TFANG(seed=42)
parent = np.array([0.0, 0.0, 5.0, 12.0])
masses = np.array([1.0, 1.0, 1.0, 1.0, 1.0])
gen.SetDecay(parent, masses)

# Full phase space volume
volume, error = gen.GetPhaseSpace(nEvents=1_000_000)
print(f"Full phase space: {volume:.1f} +/- {error:.1f}")

# Add a detector, then compute partial phase space
gen.AddDetector(CircleDetector(direction=np.array([0.0, 0.0, 1.0]), solid_angle=0.5))
partial_vol, partial_err = gen.GetPartialPhaseSpace(nEvents=1_000_000)
print(f"Partial phase space: {partial_vol:.4f} +/- {partial_err:.4f}")

Parallel Execution

For large-scale simulations, use Python's multiprocessing to distribute event generation across CPU cores:

import multiprocessing as mp
from pyFANG import TFANG
import numpy as np

def generate_worker(args):
    n_events, seed = args
    gen = TFANG(seed=seed)
    gen.SetDecay(np.array([0.0, 0.0, 5.0, 12.0]), np.array([1.0]*5))
    momenta, weights = gen.GenerateBatch(n_events)
    return weights

n_workers = 8
events_per_worker = 1_000_000

with mp.Pool(n_workers) as pool:
    results = pool.map(
        generate_worker,
        [(events_per_worker, seed) for seed in range(n_workers)]
    )

all_weights = np.concatenate(results)
print(f"Generated {len(all_weights):,} events across {n_workers} workers")
print(f"Phase space: {all_weights.mean():.2f} +/- {all_weights.std()/np.sqrt(len(all_weights)):.2f}")

A ready-to-use parallel script is included at tests/runFANG_parallel.py:

python tests/runFANG_parallel.py --workers 8

API Reference

`TFANG(rng=None, seed=None)`

Create a new generator instance.

Parameter	Type	Description
`rng`	`np.random.Generator`, optional	External RNG to use. If provided, pyFANG will not own it.
`seed`	`int`, optional	Seed for internal RNG. Defaults to 0 if `rng` is not provided.

Detector Classes

Class	Parameters	Description
`CircleDetector(direction, solid_angle)`	`direction`: 3D vector, `solid_angle`: float (steradians)	Circular cone detector
`PointDetector(direction)`	`direction`: 3D vector	Fixed-direction point detector
`StripDetector(direction, solid_angle, azimuthal_coverage)`	`direction`: 3D vector, `solid_angle`: float, `azimuthal_coverage`: float (0-1)	Rectangular strip detector
`RingDetector(direction, solid_angle, opening)`	`direction`: 3D vector, `solid_angle`: float, `opening`: float (positive)	Annular ring detector

Configuration Methods

Method	Description
`SetDecay(S, masses) -> bool`	Configure the decay. `S` is the parent 4-momentum `[px, py, pz, mass]`, `masses` is an array of daughter masses (must have at least 2 elements). Returns `True` on success.
`AddDetector(detector)`	Add a detector constraint for the next particle. Pass a `CircleDetector`, `PointDetector`, `StripDetector`, or `RingDetector` instance.
`ClearDetectors()`	Remove all detector constraints.
`SetSeed(seed)`	Set a new RNG seed.

Generation Methods

Method	Description
`Generate() -> int`	Generate one event. Returns the number of kinematic solutions found.
`GenerateWithRandoms(randoms) -> int`	Generate one event using a pre-generated array of random numbers. Useful for reproducibility.
`GenerateBatch(nEvents, seed=None) -> (momenta, weights)`	Generate many events at once, with or without detectors. Returns momenta `(N, nBody, 4)` and weights `(N,)`.

Result Access Methods

Method	Description
`GetNSolutions() -> int`	Number of solutions from the last `Generate()` call.
`GetWeight(iSolution=0) -> float`	Phase space weight for solution `iSolution`.
`GetDecay(iSolution, iParticle) -> np.ndarray`	4-momentum of particle `iParticle` in solution `iSolution`. Can also be called as `GetDecay(iParticle)` to use solution 0.
`GetDecays(iSolution) -> list[np.ndarray]`	All daughter 4-momenta for solution `iSolution`.

Phase Space Calculation

Method	Description
`GetPhaseSpace(nEvents=1000000) -> (volume, error)`	Estimate the full (unconstrained) phase space volume.
`GetPartialPhaseSpace(nEvents=1000000) -> (volume, error)`	Estimate the partial (detector-constrained) phase space volume. Requires detectors to be set.

Query Methods

Method	Description
`GetNBody() -> int`	Number of daughter particles.
`GetNDetectors() -> int`	Number of detectors set.
`IsConstrained() -> bool`	Whether any detectors are set.
`GetNRandomsPerEvent() -> int`	Number of random numbers needed per event.

4-Momentum Format

All 4-momenta in pyFANG use the format:

[px, py, pz, mass]

where px, py, pz are the three-momentum components and mass is the invariant mass. This applies to:

The parent 4-momentum in SetDecay
All daughter 4-momenta returned by GetDecay, GetDecays, and GenerateBatch

Energy is not stored explicitly. Compute it as:

energy = np.sqrt(px**2 + py**2 + pz**2 + mass**2)

Examples

The tests/ directory contains full working examples:

tests/runFANG.py -- Validates pyFANG against known results:
1. Full phase space volume for a 5-body decay (compared to published Table I)
2. Partial phase space with detector constraints (constrained vs. unconstrained-with-cuts)
3. Elastic electron-proton scattering cross section (compared to the Rosenbluth formula)
tests/runFANG_parallel.py -- Parallel event generation using multiprocessing

Run the validation:

python tests/runFANG.py

Citation

If you use pyFANG in your research, please cite the following paper:

Horin, I., Kreisel, A. & Alon, O. Focused angular N-body event generator (FANG). J. High Energ. Phys. 2025, 137 (2025).

DOI: 10.1007/JHEP12(2025)137

arXiv: 2509.11105

BibTeX:

@article{Horin2025FANG,
    author  = {Horin, Itay and Kreisel, Arik and Alon, or},
    title   = {Focused angular N-body event generator (FANG)},
    journal = {Journal of High Energy Physics},
    volume  = {2025},
    pages   = {137},
    year    = {2025},
    doi     = {10.1007/JHEP12(2025)137},
    eprint  = {2509.11105},
    archivePrefix = {arXiv},
    primaryClass  = {hep-ph}
}

Project Structure

FANG/
├── pyFANG/
│   ├── __init__.py         # Package exports
│   └── core.py             # Core FANG implementation (Numba JIT)
├── tests/
│   ├── runFANG.py          # Validation and demonstration
│   └── runFANG_parallel.py # Parallel execution example
├── pyproject.toml          # Package metadata
├── setup.py                # Setup script
├── requirements.txt        # Dependencies
├── LICENSE                 # MIT License
└── README.md               # This file

License

MIT License. See LICENSE for details.

Authors

Itay Horin
Arik Kreisel

Name		Name	Last commit message	Last commit date
Latest commit History 11 Commits
.github/workflows		.github/workflows
pyFANG		pyFANG
tests		tests
.gitignore		.gitignore
LICENSE		LICENSE
QUICKSTART.md		QUICKSTART.md
README.md		README.md
pyproject.toml		pyproject.toml
requirements.txt		requirements.txt
tutorial.ipynb		tutorial.ipynb

Folders and files

Latest commit

History

Repository files navigation

pyFANG

🚀 Quick Start & Tutorial

Table of Contents

What Does pyFANG Do?

Features

Installation

Requirements

Install from PyPI

Install from Source

Quick Start

Usage Guide

Basic Phase Space Generation

Adding Detectors

Detector Types

CircleDetector -- Circular cone

PointDetector -- Fixed direction

StripDetector -- Rectangular strip

RingDetector -- Annular ring

Mixing detector types

Batch Generation

Calculating Phase Space Volumes

Parallel Execution

API Reference

TFANG(rng=None, seed=None)

Detector Classes

Configuration Methods

Generation Methods

Result Access Methods

Phase Space Calculation

Query Methods

4-Momentum Format

Examples

Citation

Project Structure

License

Authors

About

Resources

License

Uh oh!

Stars

Watchers

Forks

Releases 1

Packages 0

Uh oh!

Contributors

Uh oh!

Languages

`CircleDetector` -- Circular cone

`PointDetector` -- Fixed direction

`StripDetector` -- Rectangular strip

`RingDetector` -- Annular ring

`TFANG(rng=None, seed=None)`

Packages