---

Relativistic Phase Space

In relativistic kinematics, the momenta of the final-state particles cannot be chosen arbitrarily. They are constrained by the particle masses, the requirement of positive energy, and the conservation of total 4-momentum. The set of all final states that satisfy these conditions is called the phase space.

Phase space appears whenever one asks which final-state configurations are kinematically allowed. It is therefore a kinematic concept: it specifies the allowed region of momentum space before any detailed interaction dynamics are introduced.

1. Kinematic Constraints

For a final state with particles of 4-momenta $p_1,\dots,p_n$, the allowed region is determined by three basic conditions:

1. each final-state particle must have its physical mass;
2. each final-state particle must have positive energy;
3. the total 4-momentum must be conserved.

Explicitly, these conditions mean

$$ p_i^2 = m_i^2, \qquad E_i > 0, \qquad P = \sum_{i=1}^n p_i. $$

The first relation means that each final-state particle is treated as a real particle with mass (m_i). Its energy is therefore not independent of its momentum, but fixed by

$$ E_i = \sqrt{\vec p_i^{\,2} + m_i^2}. $$

In practice, this means that once the masses and 3-momenta are specified, the energies are already determined. Together with total energy-momentum conservation, this produces strong correlations among the final-state particles.

2. Phase Space and the Matrix Element

For a decay or reaction, the differential rate or cross section is generally written in the form

$$ d\Gamma \;\; \text{or} \;\; d\sigma \;\propto\; |\mathcal M|^2\, d\Phi_n. $$

The two factors in this expression play distinct roles. The phase-space element $d\Phi_n$ determines the kinematically allowed region of the final state, while the matrix element $|\mathcal M|^2$, as specified by the theoretical model of the interaction, determines how the decay or reaction rate is distributed within that region.

This distinction is essential. Phase space tells us which final-state configurations are possible from kinematics alone, whereas the matrix element tells us how the underlying interaction favors or suppresses different allowed configurations. As a result, different theoretical models can lead to different event distributions even within the same kinematically allowed region.

Lorentz-Invariant Phase-Space Element

The Lorentz-invariant (n)-body phase-space element is

$$ d\Phi_n(P;p_1,\dots,p_n) = \delta^{(4)}\!\left(P-\sum_{i=1}^{n} p_i\right) \prod_{i=1}^{n}\frac{d^3p_i}{(2\pi)^3\,2E_i}. $$

This expression summarizes the kinematic constraints in compact form.

• The four-dimensional delta function enforces exact conservation of total energy and momentum.
• The factor $d^3p_i/(2E_i)$ is the Lorentz-invariant measure for each final-state particle.
• Their product describes the full region of momentum space allowed by the kinematics.

For practical use, the key point is simple: phase space is the set of final-state momenta that satisfy all physical constraints simultaneously.

Multi-Body Phase Space

For a two-body final state, the kinematics are highly constrained once the parent 4-momentum and daughter masses are fixed.

For three-body and higher final states, a continuous region of allowed configurations remains. This is the multi-body phase space.

Multi-body phase space is important for studying:

• kinematic distributions;
• invariant-mass spectra;
• detector acceptance;
• reconstruction effects;
• background shapes.

The same event can be described in different frames, such as the laboratory frame, the parent rest frame, or the center-of-mass frame. Lorentz-invariant quantities, such as invariant masses, remain unchanged under these frame transformations.

---

Using TGenPhaseSpace in ROOT

TGenPhaseSpace is a ROOT class for generating kinematically allowed multi-body decay configurations.

Its role should be clearly limited:

TGenPhaseSpace is a phase-space generator, not a reaction or decay dynamics model.

It takes as input:

• the parent 4-vector;
• the daughter particle masses.

Then it generates daughter 4-vectors satisfying:

$$ P_{\text{parent}} = \sum_i p_i, \qquad p_i^2=m_i^2, \qquad E_i>0. $$

A typical workflow is:

1. define the parent 4-vector;
2. define the daughter masses;
3. call SetDecay(...);
4. call Generate() for each event;
5. obtain daughter 4-vectors with GetDecay(i).

Practical Meaning of the Event Weight

Each call to

weight = gen.Generate();

returns an unnormalized phase-space weight for the generated event.

In the examples below, histograms are filled using this returned weight:

hist->Fill(x, weight);

This is the standard weighted-event usage of TGenPhaseSpace.

Simulation

In a full dynamical calculation, either for direct or sequential decay, the event distribution should be written as $$ d\Gamma \propto |\mathcal M|^2 d\Phi_3 . $$ Here the three-body phase-space factor $d\Phi_3$ is generated numerically with TGenPhaseSpace, which produces final-state four-momenta satisfying the on-shell conditions and total four-momentum conservation.

In the simplified code below, the direct case assumes $|\mathcal M|^2=1$, corresponding to pure three-body phase space. The sequential case uses a practical resonance model: the $^{8}\mathrm{Be}(2^+)$ mass is sampled from a Breit-Wigner distribution, followed by two successive two-body phase-space decays. This sequential treatment is therefore a practical event-generation prescription, not a full coherent three-body matrix-element calculation.

Direct Three-Alpha Decay

As a first step, we consider the direct three-alpha decay channel of the resonant reaction

$$ p + {}^{11}\mathrm{B} \rightarrow {}^{12}\mathrm{C}^{*}(16.62~\mathrm{MeV},2^-) \rightarrow \alpha + \alpha + \alpha . $$

The event generation is performed in the rest frame of the excited carbon nucleus ${}^{12}\mathrm{C}^{*}$. The $z$-axis is chosen along the incident proton-beam direction. The decay

$$ {}^{12}\mathrm{C}^{*} \rightarrow \alpha + \alpha + \alpha $$

is generated directly using three-body phase space.

After the three alpha particles are generated in the ${}^{12}\mathrm{C}^{*}$ rest frame, each event is boosted to a moving frame. This moving frame corresponds to the recoil motion of the ${}^{12}\mathrm{C}^{*}$ system for a proton beam kinetic energy

$$ E_p = 0.675~\mathrm{MeV} $$

in the $p + {}^{11}\mathrm{B}$ entrance channel, with the ${}^{11}\mathrm{B}$ target at rest.

For each alpha particle, the kinetic energy

$$ T_i = E_i - m_\alpha $$

and the angular variables

$$ (\theta_i,\phi_i) $$

are recorded both in the ${}^{12}\mathrm{C}^{*}$ rest frame and in the moving frame. The three alpha particles are stored in the order returned by the phase-space generator. No event-by-event energy sorting or angular sorting is applied.

The event weight returned by Generate() is also stored.

In direct three-body decay, the three alpha particles are ordered as returned by TGenPhaseSpace. This order has no physical significance; it is merely a convenience for storing the generated values.

%%cpp -d

#include "TFile.h"
#include "TTree.h"
#include "TGenPhaseSpace.h"
#include "TLorentzVector.h"
#include "TVector3.h"
#include "TMath.h"
#include "TRandom3.h"
#include "TCanvas.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TStyle.h"

#include <iostream>
#include <cmath>

double PhiDeg(const TLorentzVector &p)
{
    double phi = p.Phi() * TMath::RadToDeg();
    if (phi < 0.0) phi += 360.0;
    return phi;
}

void gen_direct_3alpha(Long64_t N = 500000)
{
    // ------------------------------------------------------------
    // Units:
    //   internal calculation: GeV
    //   stored kinetic energy: MeV
    //   stored angles: degree
    //
    // Reaction:
    //   p + 11B -> 12C*(16.62 MeV, 2-) -> alpha + alpha + alpha
    //
    // The three alpha particles are stored in the order returned by
    // TGenPhaseSpace. No event-by-event reordering is applied.
    // ------------------------------------------------------------

    const double mp     = 0.9382720813;  // GeV
    const double mAlpha = 3.727379378;   // GeV

    // 12C*(16.62 MeV), Q value relative to 3-alpha threshold
    const double Q3a   = 0.009345;       // GeV
    const double mC12s = 3.0 * mAlpha + Q3a;

    // Moving frame corresponding to proton beam energy Ep = 0.675 MeV
    const double Tp    = 0.000675;       // GeV
    const double pBeam = std::sqrt(Tp * (Tp + 2.0 * mp));
    const double EMove = std::sqrt(mC12s * mC12s + pBeam * pBeam);

    TVector3 betaMove(0.0, 0.0, pBeam / EMove);

    // 12C* rest frame
    TLorentzVector C12_rest(0.0, 0.0, 0.0, mC12s);

    TGenPhaseSpace gen;
    double masses[3] = {mAlpha, mAlpha, mAlpha};

    if (!gen.SetDecay(C12_rest, 3, masses)) {
        std::cerr << "SetDecay failed for direct decay." << std::endl;
        return;
    }

    TFile *fout = new TFile("direct_3alpha.root", "RECREATE");
    TTree *tree = new TTree("events", "Direct 3alpha decay");

    double T_rest[3], theta_rest[3], phi_rest[3];
    double T_move[3], theta_move[3], phi_move[3];
    double weight;

    tree->Branch("T_rest",     T_rest,     "T_rest[3]/D");
    tree->Branch("theta_rest", theta_rest, "theta_rest[3]/D");
    tree->Branch("phi_rest",   phi_rest,   "phi_rest[3]/D");

    tree->Branch("T_move",     T_move,     "T_move[3]/D");
    tree->Branch("theta_move", theta_move, "theta_move[3]/D");
    tree->Branch("phi_move",   phi_move,   "phi_move[3]/D");

    tree->Branch("weight",     &weight,    "weight/D");

    for (Long64_t iev = 0; iev < N; ++iev) {

        weight = gen.Generate();

        for (int i = 0; i < 3; ++i) {

            TLorentzVector aRest = *gen.GetDecay(i);

            // In 12C* rest frame
            T_rest[i]     = (aRest.E() - mAlpha) * 1000.0;
            theta_rest[i] = aRest.Theta() * TMath::RadToDeg();
            phi_rest[i]   = PhiDeg(aRest);

            // In moving frame
            TLorentzVector aMove = aRest;
            aMove.Boost(betaMove);

            T_move[i]     = (aMove.E() - mAlpha) * 1000.0;
            theta_move[i] = aMove.Theta() * TMath::RadToDeg();
            phi_move[i]   = PhiDeg(aMove);
        }

        tree->Fill();
    }

    tree->Write();
    fout->Close();

    std::cout << "Wrote direct_3alpha.root with " << N << " events." << std::endl;
}

Sequential Three-Alpha Decay

As an alternative decay mechanism, we consider the sequential three-alpha decay channel of the same resonant reaction,

$$ p + {}^{11}\mathrm{B} \rightarrow {}^{12}\mathrm{C}^{*}(16.62~\mathrm{MeV},2^-) \rightarrow \alpha_0 + {}^{8}\mathrm{Be}(2^+) \rightarrow \alpha_0 + \alpha_1 + \alpha_2 . $$

The event generation is again performed starting from the rest frame of the excited carbon nucleus ${}^{12}\mathrm{C}^{*}$, with the $z$-axis chosen along the incident proton-beam direction.

The sequential decay is generated in two steps. First,

$$ {}^{12}\mathrm{C}^{*} \rightarrow \alpha_0 + {}^{8}\mathrm{Be}(2^+) $$

is generated in the ${}^{12}\mathrm{C}^{*}$ rest frame. The mass of the intermediate ${}^{8}\mathrm{Be}(2^+)$ state is sampled event by event from a Breit-Wigner distribution, with values restricted to the kinematically allowed region.

Second,

$$ {}^{8}\mathrm{Be}(2^+) \rightarrow \alpha_1 + \alpha_2 $$

is generated. The resulting final-state alpha particles are expressed in the original ${}^{12}\mathrm{C}^{*}$ rest frame. The full three-alpha final state is then boosted to the same moving frame used in the direct-decay case. This moving frame corresponds to a proton beam kinetic energy

$$ E_p = 0.675~\mathrm{MeV} $$

for the $p + {}^{11}\mathrm{B}$ entrance channel with the target at rest.

For each of the three alpha particles, the kinetic energy

$$ T_i = E_i - m_\alpha $$

and the angular variables

$$ (\theta_i,\phi_i) $$

are recorded both in the ${}^{12}\mathrm{C}^{*}$ rest frame and in the moving frame.

The generator internally distinguishes the primary alpha particle from the two alpha particles originating from the $^8Be$ decay. However, this information is not experimentally observable, because the final state contains three identical alpha particles. Therefore, before the event is written to the output tree, the three alpha particles are randomly permuted.

The total event weight is taken as the product of the two phase-space weights,

$$ w = w_1 w_2 , $$

where $w_1$ is the weight from the first decay and $w_2$ is the weight from the second decay.

void gen_sequential_3alpha(Long64_t N = 500000)
{
    // ------------------------------------------------------------
    // Units:
    //   internal calculation: GeV
    //   stored kinetic energy: MeV
    //   stored angles: degree
    //
    // Reaction:
    //   p + 11B -> 12C*(16.62 MeV, 2-)
    //             -> alpha0 + 8Be(2+)
    //             -> alpha0 + alpha1 + alpha2
    //
    // alpha[0] = primary alpha from 12C* decay
    // alpha[1], alpha[2] = two alphas from 8Be decay
    //
    // No event-by-event reordering is applied.
    // ------------------------------------------------------------

    const double mp     = 0.9382720813;  // GeV
    const double mAlpha = 3.727379378;   // GeV

    // 12C*(16.62 MeV), Q value relative to 3-alpha threshold
    const double Q3a   = 0.009345;       // GeV
    const double mC12s = 3.0 * mAlpha + Q3a;

    // Moving frame corresponding to proton beam energy Ep = 0.675 MeV
    const double Tp    = 0.000675;       // GeV
    const double pBeam = std::sqrt(Tp * (Tp + 2.0 * mp));
    const double EMove = std::sqrt(mC12s * mC12s + pBeam * pBeam);

    TVector3 betaMove(0.0, 0.0, pBeam / EMove);

    // 8Be(2+) parameters
    const double E8gs_above_2a = 0.00009184;  // GeV
    const double Ex8_2plus     = 0.003030;    // GeV
    const double Gamma8_2plus  = 0.001513;    // GeV

    const double m8_mean = 2.0 * mAlpha + E8gs_above_2a + Ex8_2plus;
    const double m8_min  = 2.0 * mAlpha;
    const double m8_max  = mC12s - mAlpha;

    TLorentzVector C12_rest(0.0, 0.0, 0.0, mC12s);

    TGenPhaseSpace gen1;
    TGenPhaseSpace gen2;
    TRandom3 rng(0);

    TFile *fout = new TFile("sequential_3alpha.root", "RECREATE");
    TTree *tree = new TTree("events", "Sequential 3alpha decay");

    double T_rest[3], theta_rest[3], phi_rest[3];
    double T_move[3], theta_move[3], phi_move[3];
    double weight;

    tree->Branch("T_rest",     T_rest,     "T_rest[3]/D");
    tree->Branch("theta_rest", theta_rest, "theta_rest[3]/D");
    tree->Branch("phi_rest",   phi_rest,   "phi_rest[3]/D");

    tree->Branch("T_move",     T_move,     "T_move[3]/D");
    tree->Branch("theta_move", theta_move, "theta_move[3]/D");
    tree->Branch("phi_move",   phi_move,   "phi_move[3]/D");

    tree->Branch("weight",     &weight,    "weight/D");

    for (Long64_t iev = 0; iev < N; ++iev) {

        // Sample 8Be(2+) mass with a truncated Breit-Wigner distribution
        double m8 = 0.0;

        while (true) {
            double u = rng.Uniform();
            m8 = m8_mean + 0.5 * Gamma8_2plus * std::tan(TMath::Pi() * (u - 0.5));

            if (m8 > m8_min && m8 < m8_max) break;
        }

        // First decay:
        // 12C* -> alpha0 + 8Be
        double masses1[2] = {mAlpha, m8};

        if (!gen1.SetDecay(C12_rest, 2, masses1)) {
            --iev;
            continue;
        }

        double w1 = gen1.Generate();

        TLorentzVector alpha[3];

        alpha[0] = *gen1.GetDecay(0);
        TLorentzVector be8 = *gen1.GetDecay(1);

        // Second decay:
        // 8Be -> alpha1 + alpha2
        double masses2[2] = {mAlpha, mAlpha};

        if (!gen2.SetDecay(be8, 2, masses2)) {
            --iev;
            continue;
        }

        double w2 = gen2.Generate();

        alpha[1] = *gen2.GetDecay(0);
        alpha[2] = *gen2.GetDecay(1);

        weight = w1 * w2;
        // --------------------------------------------------------
        // Random permutation to erase generator labels
        // --------------------------------------------------------
        for (int k = 2; k > 0; --k) {
            int j = int(rng.Uniform() * (k + 1));
            TLorentzVector tmp = alpha[k];
            alpha[k] = alpha[j];
            alpha[j] = tmp;
        }

        for (int i = 0; i < 3; ++i) {

            TLorentzVector aRest = alpha[i];

            // In 12C* rest frame
            T_rest[i]     = (aRest.E() - mAlpha) * 1000.0;
            theta_rest[i] = aRest.Theta() * TMath::RadToDeg();
            phi_rest[i]   = PhiDeg(aRest);

            // In moving frame
            TLorentzVector aMove = aRest;
            aMove.Boost(betaMove);

            T_move[i]     = (aMove.E() - mAlpha) * 1000.0;
            theta_move[i] = aMove.Theta() * TMath::RadToDeg();
            phi_move[i]   = PhiDeg(aMove);
        }

        tree->Fill();
    }

    tree->Write();
    fout->Close();

    std::cout << "Wrote sequential_3alpha.root with " << N << " events." << std::endl;
}

gen_sequential_3alpha();

Wrote sequential_3alpha.root with 500000 events.

Comparison of the Two Decay Modes

TCanvas *c = nullptr;

void plot_3alpha_correlations(const char *filename, const char *tag)
{
    TFile *fin = TFile::Open(filename);
    TTree *tree = (TTree*)fin->Get("events");

    gStyle->SetOptStat(0);
    
    if (c) {
        delete c;
        c = nullptr;
    }    
   
    c = new TCanvas("c", "c", 900, 450);

    c->Divide(2, 1);

    c->cd(1);
    gPad->SetRightMargin(0.14);

    TH2D *h_rest_01 = new TH2D(
        Form("h_rest_01_%s", tag),
        Form("%s: rest frame [0]-[1];T_{0} (MeV);T_{1} (MeV)", tag), 260, 0.0, 9, 260, 0.0, 9);

    tree->Draw(
        Form("T_rest[1]:T_rest[0]>>h_rest_01_%s", tag), "weight", "colz");


    c->cd(2);
    gPad->SetRightMargin(0.14);

    TH2D *h_move_01 = new TH2D(
        Form("h_move_01_%s", tag),
        Form("%s: moving frame [0]-[1];T_{0} (MeV);T_{1} (MeV)", tag), 260, 0.0, 9, 260, 0.0, 9);

    tree->Draw(
        Form("T_move[1]:T_move[0]>>h_move_01_%s", tag), "weight", "colz");


    c->Draw();
}

void plot_all_alpha_spectra(const char *filename, const char *tag)
{
    TFile *fin = TFile::Open(filename);
    TTree *tree = (TTree*)fin->Get("events");

    gStyle->SetOptStat(0);

    if (c) {
        delete c;
        c = nullptr;
    }    
   
    c = new TCanvas("c", "c", 900, 450);

    c->Divide(2, 1);

    c->cd(1);

    TH1D *h_rest = new TH1D(
        Form("h_all_rest_%s", tag),
        Form("%s: all alpha, rest frame;T_{#alpha} (MeV);Counts", tag), 260, 0.0, 9);

    tree->Draw(Form("T_rest>>h_all_rest_%s", tag),  "weight", "hist");
    h_rest->Draw("hist");

    c->cd(2);

    TH1D *h_move = new TH1D(
        Form("h_all_move_%s", tag),
        Form("%s: all alpha, moving frame;T_{#alpha} (MeV);Counts", tag), 260, 0.0, 9);

    tree->Draw(Form("T_move>>h_all_move_%s", tag),  "weight", "hist");
    h_move->Draw("hist");

    c->Draw();
}

Inclusive alpha spectrum

Direct

plot_all_alpha_spectra("direct_3alpha.root", "Direct");

Sequential

plot_all_alpha_spectra("sequential_3alpha.root", "Sequential");

Two-dimensional energy-correlation

Direct

plot_3alpha_correlations("direct_3alpha.root", "Direct");

Sequeintial

plot_3alpha_correlations("sequential_3alpha.root", "Sequential");

---

Appendix: Practical Notes on Using TGenPhaseSpace

1. Check that the decay is kinematically allowed

Before generating events, first verify that the decay is allowed by kinematics. In practice, this means that the parent four-momentum must contain enough invariant mass to produce all daughter particles. In ROOT, this check is performed by SetDecay(...), which returns kTRUE if the decay is permitted and kFALSE otherwise. (ROOT)

TGenPhaseSpace phaseSpaceGen;

if (!phaseSpaceGen.SetDecay(parent, numDaughters, daughterMasses)) {
    std::cerr << "Decay is kinematically forbidden." << std::endl;
    return;
}

2. Understand what Generate() returns

Once SetDecay(...) succeeds, each call to Generate() produces one random kinematically allowed final state. ROOT explicitly states that the return value is the un-normalized weight of the current event. The daughter four-momenta are then obtained with GetDecay(i). (ROOT)

double weight = phaseSpaceGen.Generate();
TLorentzVector* daughterMomentum = phaseSpaceGen.GetDecay(i);

For many simple phase-space studies, this is all that is needed: generate an event, read the daughter four-vectors, and fill the desired observables.

3. Weighted events are the default

The most straightforward use of TGenPhaseSpace is to keep the generated events weighted and use the returned event weight when filling histograms. This is also how the official ROOT tutorial example is written. (ROOT)

hist->Fill(x, weight);

For teaching examples, quick checks of kinematic distributions, and qualitative comparisons, this is usually the simplest and most robust approach.

4. Additional care is needed for unweighted samples

If an unweighted event sample is required, a rejection-sampling step must be added. In that case, one needs a value $W_{\max}$ such that

$$ W \le W_{\max} $$

for every generated event. ROOT provides GetWtMax(), and the class stores an internal maximum-weight quantity fWtMax; in the implementation, Generate() uses that internal factor when constructing the event weight. (ROOT)

For a custom unweighting procedure, a practical approach is to estimate a working $W_{\max}$ from a large pilot run and then include a small safety margin. The key requirement is that the chosen $W_{\max}$ must remain a true upper bound for all accepted events.

TGenPhaseSpace phaseSpaceGen;
if (!phaseSpaceGen.SetDecay(parent, numDaughters, daughterMasses)) return;

double Wmax = 0.0;
const int Nsamples = 1000000;

for (int i = 0; i < Nsamples; ++i) {
    double weight = phaseSpaceGen.Generate();
    if (weight > Wmax) Wmax = weight;
}

Wmax *= 1.05;  // small safety margin

5. Use rejection sampling consistently for unweighted events

Once a suitable $W_{\max}$ has been chosen, unweighted events can be generated by accepting each trial event with probability $W/W_{\max}$. A simple implementation is:

while (true) {
    double weight = phaseSpaceGen.Generate();
    double u = gRandom->Rndm();

    if (u < weight / Wmax) {
        TLorentzVector* daughterMomentum = phaseSpaceGen.GetDecay(i);
        break;
    }
}

This produces accepted events with equal statistical weight, provided that $W_{\max}$ is indeed an upper bound. If a later event exceeds the chosen $W_{\max}$, then the sample is no longer strictly correct and the procedure should be restarted with a larger bound.

6. Recommended practical usage

For most exploratory studies, phase-space checks, and pedagogical examples, the safest workflow is:

1. call SetDecay(...) and verify that it succeeds;
2. call Generate() to obtain one event at a time;
3. access daughter four-momenta through GetDecay(i);
4. keep the events weighted unless an unweighted sample is specifically required.