核物理实验数据处理 Experimental Data Analysis in Nuclear Physics

---

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) = (2\pi)^4\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.

In this section the goal is to compare distribution shapes, not absolute reaction yields. The event weights of separately configured generators do not set the relative branching ratio between two decay mechanisms. When overlaying their spectra, normalize each to unit area; keep the stored event weights when filling each spectrum.

%jsroot on

Direct Three-Alpha Decay

We study $^{12}\mathrm C^*\to3\alpha$. The parent mass is fixed so that the available kinetic energy in its rest frame is $Q_{3\alpha}=9.345$ MeV, corresponding to $E_x\simeq16.620$ MeV. The direct sample uses constant-matrix-element three-body phase space; the sequential sample proceeds through a broad $^{8}\mathrm{Be}(2^+)$ state. These are two teaching models of the decay, not a calculation of their physical branching fractions.

Both samples are generated in the parent rest frame and then boosted along $+z$. The boost uses a momentum scale computed from a 0.675 MeV proton. Because the parent mass is fixed independently, this is a specified moving frame, not an exact reconstruction of the $p+{}^{11}\mathrm B$ entrance channel. It illustrates how energies and angles change while invariant masses do not.

The output tree events stores kinetic energies in MeV and angles in degrees: T_rest, theta_rest, phi_rest refer to the parent rest frame; T_move, theta_move, phi_move refer to the boosted frame. Each branch contains three elements, one per alpha particle. weight is used when filling distributions.

In the direct sample, the three labels are the order returned by GetDecay(i) and have no experimental identity. GetDecay(i) returns a pointer to generator-owned storage. Copy the four-vector before the next Generate() call.

%%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;
}

%%cpp -d
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.

    gRandom->SetSeed(4201); // Generate 使用 gRandom；固定种子复现相同事件

    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

We study $^{12}\mathrm C^*\to3\alpha$. The parent mass is fixed so that the available kinetic energy in its rest frame is $Q_{3\alpha}=9.345$ MeV, corresponding to $E_x\simeq16.620$ MeV. The direct sample uses constant-matrix-element three-body phase space; the sequential sample proceeds through a broad $^{8}\mathrm{Be}(2^+)$ state. These are two teaching models of the decay, not a calculation of their physical branching fractions.

Both samples are generated in the parent rest frame and then boosted along $+z$. The boost uses a momentum scale computed from a 0.675 MeV proton. Because the parent mass is fixed independently, this is a specified moving frame, not an exact reconstruction of the $p+{}^{11}\mathrm B$ entrance channel. It illustrates how energies and angles change while invariant masses do not.

The output tree events stores kinetic energies in MeV and angles in degrees: T_rest, theta_rest, phi_rest refer to the parent rest frame; T_move, theta_move, phi_move refer to the boosted frame. Each branch contains three elements, one per alpha particle. weight is used when filling distributions.

For each sequential event, sample the intermediate mass from a Breit–Wigner distribution restricted to $2m_\alpha<m_8<M_{12}-m_\alpha$, then generate the two successive two-body decays. The width parameter is the Breit–Wigner FWHM. Since the second SetDecay receives the already moving $^8$Be four-vector, its daughters are returned in the $^{12}$C rest frame; no extra boost is needed at that step.

The code randomly permutes the three final labels before storage. This avoids tying an array index to the primary alpha, but does not implement quantum-mechanical symmetrization or interference. Both steps assume isotropic decay and omit spin-dependent angular correlations. In ROOT's default two-body mode each normalized generation weight is 1 (up to roundoff); their product does not add a physical mass-dependent resonance probability. Here that mass distribution is precisely the truncated Breit–Wigner distribution just specified.

%%cpp -d
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
    //
    // Randomly permute the final labels before storing the event.

    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(4202); // 抽样中间态质量、打乱三个 alpha 的标签
    gRandom->SetSeed(4203); // 两级 TGenPhaseSpace 使用的随机数序列

    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_direct_3alpha();
gen_sequential_3alpha();

Wrote direct_3alpha.root with 500000 events.
Wrote sequential_3alpha.root with 500000 events.

Comparison of the Two Decay Modes

Run both generator calls before making these plots. tree->Draw("T_rest", "weight") fills all three array elements, so an inclusive spectrum has three alpha entries per event. In "T_rest[1]:T_rest[0]", ROOT places the expression before the colon on the vertical axis. Compare the rest-frame and moving-frame panels within each model; their event identities and weights are unchanged.

%%cpp -d
TCanvas *c = nullptr;

%%cpp -d
void plot_3alpha_correlations(const char *filename, const char *tag)
{
    TFile *fin = TFile::Open(filename);
    if (!fin || fin->IsZombie()) throw std::runtime_error("Cannot open input ROOT file");
    TTree *tree = fin->Get<TTree>("events");
    if (!tree) throw std::runtime_error("Missing events tree");

    gStyle->SetOptStat(0);
    gStyle->SetPadLeftMargin(0.15);
    gStyle->SetPadRightMargin(0.14);
    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();
}

%%cpp -d
void plot_all_alpha_spectra(const char *filename, const char *tag)
{
    TFile *fin = TFile::Open(filename);
    if (!fin || fin->IsZombie()) throw std::runtime_error("Cannot open input ROOT file");
    TTree *tree = fin->Get<TTree>("events");
    if (!tree) throw std::runtime_error("Missing events tree");

    gStyle->SetOptStat(0);
    gStyle->SetPadLeftMargin(0.15);
    gStyle->SetPadRightMargin(0.14);
    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);Sum of weights", 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);Sum of weights", 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

The direct-decay spectrum is continuous: the fixed available energy can be shared in many ways among three particles. The weighted shape is the phase-space prediction, not a detector response.

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

Sequential

The sequential sample has additional structure because two particles share the intermediate-resonance energy. Its broad width smears the feature. The vertical scales of the separate direct and sequential panels are not relative yields.

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

Two-dimensional energy-correlation

Direct

Energy conservation correlates the two displayed energies with the unplotted third energy. A boost changes these energy–energy plots even though it preserves all invariant masses.

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

Sequential

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

For an unweighted sample, accept a proposed event with probability $w/W$, where $W$ bounds every generated weight. A largest weight observed in a finite pilot sample, even multiplied by a safety factor, is not a proof of that bound.

For ROOT 6's default SetDecay(..., "") implementation, Generate() contains the GENBOD normalization factor; its normalized weight is bounded by 1. GetWtMax() returns that internal factor, not the upper bound to divide the returned weight by. The distinction follows from the ROOT implementation. This statement does not cover the "Fermi" option or additional user weights.

5. Use rejection sampling consistently for unweighted events

For the default generator used here, rejection sampling can be written as:

double w;
do {
    w = gen.Generate();
    if (w < 0.0 || w > 1.0 + 1e-12)
        throw std::runtime_error("Weight outside the assumed bound");
} while (gRandom->Rndm() > w);
TLorentzVector daughter = *gen.GetDecay(0);
// The accepted event now has unit statistical weight: do not weight it again.

Keep weighted generation for the examples above. Use this alternative only when an equal-weight event sample is needed.

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.