Invariant Mass, Missing Mass, and Dalitz Analysis

Invariant Mass

Invariant mass is a Lorentz-invariant quantity that characterizes the total energy and momentum of a particle or a system of particles. Because it has the same value in all inertial frames, it is one of the most useful observables in particle and nuclear physics.

For two particles with four-momenta $P_1=(E_1,\vec p_1)$ and $P_2=(E_2,\vec p_2)$, the invariant mass is defined as

$$ M_{\mathrm{inv}}^2=(P_1+P_2)^2 =(E_1+E_2)^2-\left|\vec p_1+\vec p_2\right|^2. $$

For a system of $n$ particles, this generalizes to

$$ M_{\mathrm{inv}}^2(n)=\left(P_1+P_2+\cdots+P_n\right)^2. $$

For a single particle, the invariant mass is simply its rest mass. For a system of final-state particles, however, the invariant mass allows us to reconstruct a short-lived parent state from its decay products.

---

Intermediate States and Resonances

In many reactions, the observed particles are produced through intermediate states. If such a state is short-lived, it is not directly seen as a stable particle in the detector. Instead, its presence is inferred from the invariant mass of its decay products.

Suppose an intermediate state $R$ decays as

$$ R \to 1 + 2. $$

Then the parent four-momentum is reconstructed as

$$ P_R=P_1+P_2, $$

and its invariant mass is

$$ M_R=\sqrt{P_R^2}. $$

If a large number of events are accumulated, a resonance appears as an enhancement or peak in the invariant-mass spectrum. This is one of the basic methods used in experimental particle physics to identify unstable states.

---

Two-Body Collision Reactions

Consider a reaction of the type

$$ A+B\to \text{products}. $$

A final-state particle pair may be produced in two different ways.

In the non-resonant case, the particles are produced directly, for example,

$$ A+B\to c+d+e. $$

Then the invariant mass of a selected pair, such as $c+d$, varies smoothly according to phase space.

In the resonant case, an intermediate state is produced first and then decays:

$$ A+B\to c+R, \qquad R\to e+f. $$

Now the invariant mass of the pair $e+f$ clusters around the mass of the resonance $R$.

This is the physical reason invariant-mass spectra are so useful: they distinguish smooth phase-space production from the production of an actual intermediate state.

---

Missing Mass

Missing mass is used when not all final-state particles are detected.

If the initial-state four-momentum is known and only part of the final state is measured, the missing four-momentum is defined as

$$ P_{\mathrm{miss}}=(P_A+P_B)-P_{\mathrm{detected}}. $$

The missing mass is then

$$ M_{\mathrm{miss}}^2=P_{\mathrm{miss}}^2. $$

So invariant mass and missing mass are closely related:

• invariant mass reconstructs a state from detected daughters;
• missing mass reconstructs an undetected particle or system from the known initial state and the measured remainder.

Both are direct applications of four-momentum conservation.

---

Example: Detecting Short-Lived Resonances

A standard example is the reaction

$$ pp \to pp\pi^+\pi^-. $$

All four final-state particles are measured. One then computes the invariant mass of the pion pair,

$$ M_{\pi^+\pi^-}=\sqrt{(P_{\pi^+}+P_{\pi^-})^2}. $$

If part of the event sample proceeds through the chain

$$ pp\to pp\rho, \qquad \rho\to\pi^+\pi^-, $$

then the $\pi^+\pi^-$ invariant-mass distribution contains a resonant enhancement near the $\rho$ mass, sitting on top of a broader non-resonant continuum from direct $pp\to pp\pi^+\pi^-$ production.

This example illustrates the essential logic of invariant-mass analysis: a resonance is identified through the kinematics of its decay products.

Code example

The following example keeps the original analysis idea, but makes the two contributions explicit:

• a non-resonant four-body channel $pp\to pp\pi^+\pi^-$;
• a resonant channel $pp\to pp\rho$, followed by $\rho\to\pi^+\pi^-$.

#include <TROOT.h>
#include <TMath.h>
#include <TRandom3.h>
#include <TH1F.h>
#include <TCanvas.h>
#include <TLorentzVector.h>
#include <TGenPhaseSpace.h>
#include <TLegend.h>
#include <iostream>

void invariant_mass_demo()
{
    TRandom3 rand(0);
    gROOT->GetListOfCanvases()->Delete();

    const double m_pion    = 0.13957;   // GeV
    const double m_proton  = 0.93827;   // GeV
    const double rho_mass  = 0.775;     // GeV
    const double rho_width = 0.149;     // GeV
    const double sqrt_s    = 7.0;       // GeV
    const int n_events     = 100000;
    const double frac_rho  = 0.10;

    TH1F *h_total = new TH1F(
        "h_total",
        "#pi^{+}#pi^{-} Invariant Mass;M_{#pi^{+}#pi^{-}} (GeV/c^{2});Events",
        100, 0.25, 1.50
    );

    TH1F *h_rho = new TH1F(
        "h_rho",
        "Resonant Contribution",
        100, 0.25, 1.50
    );

    TH1F *h_nonrho = new TH1F(
        "h_nonrho",
        "Non-resonant Contribution",
        100, 0.25, 1.50
    );

    double E_beam  = sqrt_s / 2.0;
    double pz_beam = TMath::Sqrt(E_beam * E_beam - m_proton * m_proton);

    TLorentzVector p4_p1(0, 0,  pz_beam, E_beam);
    TLorentzVector p4_p2(0, 0, -pz_beam, E_beam);
    TLorentzVector p4_initial = p4_p1 + p4_p2;

    TGenPhaseSpace phaseSpace_pprho;
    TGenPhaseSpace phaseSpace_rho;
    TGenPhaseSpace phaseSpace_full;

    double masses_pprho[3] = {m_proton, m_proton, rho_mass};
    double masses_rho[2]   = {m_pion, m_pion};
    double masses_full[4]  = {m_proton, m_proton, m_pion, m_pion};

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

        bool is_rho = (rand.Uniform() < frac_rho);

        if (is_rho) {
            double M_rho = 0.0;
            bool ok = false;

            for (int trial = 0; trial < 100; ++trial) {
                M_rho = rand.BreitWigner(rho_mass, rho_width);
                if (M_rho < 2.0 * m_pion) continue;
                if (M_rho > sqrt_s - 2.0 * m_proton) continue;

                masses_pprho[2] = M_rho;
                if (phaseSpace_pprho.SetDecay(p4_initial, 3, masses_pprho)) {
                    ok = true;
                    break;
                }
            }

            if (!ok) continue;

            double weight1 = phaseSpace_pprho.Generate();
            if (weight1 <= 0.0) continue;

            TLorentzVector p4_rho = *phaseSpace_pprho.GetDecay(2);

            if (!phaseSpace_rho.SetDecay(p4_rho, 2, masses_rho)) continue;
            double weight2 = phaseSpace_rho.Generate();
            if (weight2 <= 0.0) continue;

            TLorentzVector p4_pion1 = *phaseSpace_rho.GetDecay(0);
            TLorentzVector p4_pion2 = *phaseSpace_rho.GetDecay(1);

            double inv_mass = (p4_pion1 + p4_pion2).M();
            double weight = weight1 * weight2;

            h_rho->Fill(inv_mass, weight);
            h_total->Fill(inv_mass, weight);
        }
        else {
            if (!phaseSpace_full.SetDecay(p4_initial, 4, masses_full)) continue;

            double weight = phaseSpace_full.Generate();
            if (weight <= 0.0) continue;

            TLorentzVector p4_pion1 = *phaseSpace_full.GetDecay(2);
            TLorentzVector p4_pion2 = *phaseSpace_full.GetDecay(3);

            double inv_mass = (p4_pion1 + p4_pion2).M();

            h_nonrho->Fill(inv_mass, weight);
            h_total->Fill(inv_mass, weight);
        }
    }

    TCanvas *c1 = new TCanvas("c1", "Invariant Mass", 700, 500);

    h_total->SetLineColor(kBlack);
    h_total->SetFillColor(kGray + 1);
    h_total->SetFillStyle(1001);
    h_total->Draw("HIST");

    h_rho->SetLineColor(kBlue);
    h_rho->SetFillColor(kBlue);
    h_rho->SetFillStyle(3004);
    h_rho->Draw("HIST SAME");

    h_nonrho->SetLineColor(kRed);
    h_nonrho->SetFillColor(kRed);
    h_nonrho->SetFillStyle(3005);
    h_nonrho->Draw("HIST SAME");

    TLegend *leg = new TLegend(0.58, 0.68, 0.88, 0.88);
    leg->AddEntry(h_total,  "Total", "f");
    leg->AddEntry(h_rho,    "#rho #rightarrow #pi^{+}#pi^{-}", "f");
    leg->AddEntry(h_nonrho, "Non-resonant", "f");
    leg->Draw();
    c1->Draw();
}

invariant_mass_demo();

Interpretation

The non-resonant component produces a broad distribution determined mainly by available phase space.
The resonant component produces a localized enhancement in the $\pi^+\pi^-$ invariant-mass spectrum.
The sum of the two gives the physically more realistic picture: a resonance peak on top of a continuum background.

---

Example: Missing-Mass Reconstruction of a $\rho$ Meson

Consider the reaction

$$ pp \to pp\pi^+\pi^- . $$ Suppose only the two outgoing protons are detected, while the pion pair is not measured directly. The missing four-momentum is then defined as

$$ P_X = P_{\mathrm{initial}} - P_{p_1'} - P_{p_2'} , $$

and the corresponding missing mass is

$$ M_{\mathrm{miss}}(pp)=\sqrt{P_X^2}. $$

Here $X$ denotes the unobserved system recoiling against the two detected protons.

If the reaction proceeds through direct non-resonant production,

$$ pp \to pp\pi^+\pi^-, $$

then $X$ is simply the $\pi^+\pi^-$ system, and the missing-mass spectrum forms a smooth continuum.

If instead the reaction proceeds through an intermediate resonance,

$$ pp \to pp\rho,\qquad \rho\to\pi^+\pi^-, $$

then the same missing-mass variable reconstructs the recoiling system $X=\rho$, and the spectrum shows an enhancement near the $\rho$-meson mass.

This example illustrates the basic idea of missing-mass analysis: even when some final-state particles are not detected, an intermediate state can still be identified from the measured particles through four-momentum conservation.

Code example

---

#include <TMath.h>
#include <TRandom3.h>
#include <TH1F.h>
#include <TCanvas.h>
#include <TLine.h>
#include <TLorentzVector.h>
#include <TGenPhaseSpace.h>
#include <TStyle.h>
#include <iostream>

void missing_mass_rho_demo()
{
    gStyle->SetOptStat(0);
    gROOT->GetListOfCanvases()->Delete();

    const double m_pion    = 0.13957;   // GeV
    const double m_proton  = 0.93827;   // GeV
    const double rho_mass  = 0.775;     // GeV
    const double rho_width = 0.149;     // GeV
    const double sqrt_s    = 3.0;       // GeV
    const int n_events     = 100000;
    const double frac_rho  = 0.35;      // fraction of resonant events

    double E_beam  = sqrt_s / 2.0;
    double pz_beam = TMath::Sqrt(E_beam * E_beam - m_proton * m_proton);

    TLorentzVector p4_beam1(0, 0,  pz_beam, E_beam);
    TLorentzVector p4_beam2(0, 0, -pz_beam, E_beam);
    TLorentzVector p4_initial = p4_beam1 + p4_beam2;

    TH1F *h_mmiss = new TH1F(
        "h_mmiss",
        "pp #rightarrow ppX, X #rightarrow #pi^{+}#pi^{-};M_{miss}(pp) [GeV/c^{2}];Events",
        100, 0.25, 1.50
    );

    TRandom3 rand(0);

    // Non-resonant channel: pp -> pp pi+ pi-
    TGenPhaseSpace phaseSpace_nonres;
    double masses_nonres[4] = {m_proton, m_proton, m_pion, m_pion};

    if (!phaseSpace_nonres.SetDecay(p4_initial, 4, masses_nonres)) {
        std::cerr << "Non-resonant channel is kinematically forbidden." << std::endl;
        return;
    }

    // Resonant channel: pp -> pp rho, followed by rho -> pi+ pi-
    TGenPhaseSpace phaseSpace_res;
    TGenPhaseSpace phaseSpace_rho;
    double masses_res[3] = {m_proton, m_proton, rho_mass};
    double masses_rho[2] = {m_pion, m_pion};

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

        bool is_rho = (rand.Uniform() < frac_rho);

        if (is_rho) {
            double M_rho = 0.0;
            bool ok = false;

            for (int trial = 0; trial < 100; ++trial) {
                M_rho = rand.BreitWigner(rho_mass, rho_width);

                if (M_rho < 2.0 * m_pion) continue;
                if (M_rho > sqrt_s - 2.0 * m_proton) continue;

                masses_res[2] = M_rho;
                if (phaseSpace_res.SetDecay(p4_initial, 3, masses_res)) {
                    ok = true;
                    break;
                }
            }

            if (!ok) continue;

            double weight1 = phaseSpace_res.Generate();
            if (weight1 <= 0.0) continue;

            TLorentzVector p4_p1  = *phaseSpace_res.GetDecay(0);
            TLorentzVector p4_p2  = *phaseSpace_res.GetDecay(1);
            TLorentzVector p4_rho = *phaseSpace_res.GetDecay(2);

            // Decay rho -> pi+ pi- (not observed, but included for completeness)
            if (!phaseSpace_rho.SetDecay(p4_rho, 2, masses_rho)) continue;
            double weight2 = phaseSpace_rho.Generate();
            if (weight2 <= 0.0) continue;

            TLorentzVector p4_miss = p4_initial - p4_p1 - p4_p2;
            h_mmiss->Fill(p4_miss.M(), weight1 * weight2);
        }
        else {
            double weight = phaseSpace_nonres.Generate();
            if (weight <= 0.0) continue;

            TLorentzVector p4_p1 = *phaseSpace_nonres.GetDecay(0);
            TLorentzVector p4_p2 = *phaseSpace_nonres.GetDecay(1);

            TLorentzVector p4_miss = p4_initial - p4_p1 - p4_p2;
            h_mmiss->Fill(p4_miss.M(), weight);
        }
    }

    TCanvas *c6 = new TCanvas("c6", "Missing Mass of the Recoil System", 700, 500);
    h_mmiss->SetLineWidth(2);
    h_mmiss->Draw("HIST");

    TLine *line = new TLine(rho_mass, 0.0, rho_mass, 0.9 * h_mmiss->GetMaximum());
    line->SetLineStyle(2);
    line->Draw();

    c6->Draw();
}

missing_mass_rho_demo();

Three-Body Decays and the Dalitz Plot

A three-body decay,

$$ M\to 1+2+3, $$

contains more kinematic information than a two-body decay. For this reason, a single invariant-mass variable is usually not enough to describe the full decay structure.

The standard tool is the Dalitz plot, which gives a two-dimensional representation of the internal kinematics of a three-body final state.

---

Kinematics and Degrees of Freedom

Each of the three final-state particles has a four-momentum $P_i=(E_i,\vec p_i)$, so the final state begins with 12 kinematic quantities.

These are reduced by constraints:

• 3 from momentum conservation,
• 1 from energy conservation,
• 3 on-shell conditions $P_i^2=m_i^2$,
• 3 describing the overall orientation of the decay in space.

This leaves 2 independent internal variables.

That is why the kinematics of a three-body decay can be represented completely by a two-dimensional plot.

---

Dalitz Plot Representation

Define the pairwise invariant masses squared:

$$ m_{12}^2=(P_1+P_2)^2,\qquad m_{13}^2=(P_1+P_3)^2,\qquad m_{23}^2=(P_2+P_3)^2. $$

These satisfy

$$ m_{12}^2+m_{13}^2+m_{23}^2 = M^2+m_1^2+m_2^2+m_3^2. $$

So only two of them are independent. A Dalitz plot usually displays one pair, such as

$$ m_{12}^2 \ \text{vs.}\ m_{13}^2. $$

Each event becomes one point in this plane.

The kinematic boundary defines the allowed region, while the density inside that region reflects the underlying decay dynamics.

---

Phase Space and Dynamics

At this point it is important to distinguish two separate ideas.

Pure phase space

Pure phase space determines:

• which kinematic region is allowed;
• the baseline population of events inside that region.

Dynamics

The actual decay dynamics determine:

• whether some parts of the Dalitz plot are enhanced or suppressed;
• whether intermediate resonances are present.

An intermediate resonance in one two-body subsystem appears as a band in the Dalitz plot, because many events accumulate at the same invariant mass in that subsystem.

This separation between kinematic boundary and dynamical population is the key to reading a Dalitz plot correctly.

---

Example: $K^-p\to\Lambda\pi^+\pi^-$

A standard example is

$$ K^-p\to\Lambda\pi^+\pi^-. $$

For this reaction, define

$$ m_{12}\equiv m(\Lambda\pi^-),\qquad m_{13}\equiv m(\Lambda\pi^+). $$

Then the Dalitz plot is formed by

$$ m_{12}^2 \ \text{vs.}\ m_{13}^2. $$

This reaction is especially useful because the two subsystems

$$ \Lambda\pi^- \quad \text{and} \quad \Lambda\pi^+ $$

may contain the resonances

$$ \Sigma^-(1385)\to\Lambda\pi^-, \qquad \Sigma^+(1385)\to\Lambda\pi^+. $$

So the Dalitz plot can show two resonance bands, one vertical and one horizontal.

---

Step 1: Pure phase-space Dalitz plot

The first step is to generate the three-body final state using only phase space. This shows the allowed kinematic region without any additional dynamical enhancement.

#include <TGenPhaseSpace.h>
#include <TLorentzVector.h>
#include <TH2F.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <iostream>

void dalitz_plot_phase_space()
{
    gStyle->SetOptStat(0);
    gROOT->GetListOfCanvases()->Delete();

    const double m_lambda = 1.115683;
    const double m_pion   = 0.139570;
    const double sqrt_s   = 1.95;
    const int n_events    = 300000;

    TLorentzVector p_initial(0, 0, 0, sqrt_s);

    double masses[3] = {m_lambda, m_pion, m_pion};
    TGenPhaseSpace phase_space;

    if (!phase_space.SetDecay(p_initial, 3, masses)) {
        std::cerr << "Error: phase-space setup failed." << std::endl;
        return;
    }

    TH2F *h_dalitz_ps = new TH2F(
        "h_dalitz_ps",
        "Phase-space Dalitz Plot; (m_{#Lambda#pi^{-}})^{2} (GeV^{2}/c^{4}); (m_{#Lambda#pi^{+}})^{2} (GeV^{2}/c^{4})",
        300, 1.5, 3.5,
        300, 1.5, 3.5
    );

    for (int i = 0; i < n_events; ++i) {
        double weight = phase_space.Generate();
        if (weight <= 0.0) continue;

        TLorentzVector p_lambda   = *phase_space.GetDecay(0);
        TLorentzVector p_pi_plus  = *phase_space.GetDecay(1);
        TLorentzVector p_pi_minus = *phase_space.GetDecay(2);

        double m12_sq = (p_lambda + p_pi_minus).M2();
        double m13_sq = (p_lambda + p_pi_plus ).M2();

        h_dalitz_ps->Fill(m12_sq, m13_sq, weight);
    }

    TCanvas *c1 = new TCanvas("c1", "Dalitz Plot: Phase Space", 800, 650);
    c1->SetLogz();
    h_dalitz_ps->Draw("COLZ");
    c1->Draw();
}

dalitz_plot_phase_space();

Interpretation of Step 1

This plot shows only the kinematic boundary and the phase-space population.
It should be read as the reference picture: what the Dalitz plot looks like before any specific resonance structure is added.

---

Step 2: Add a simple $\Sigma^\pm(1385)$ enhancement

The second step is to keep the same three-body phase-space generation, but add a simple phenomenological enhancement in the $\Lambda\pi^-$ and $\Lambda\pi^+$ subsystems. This makes the resonance bands visible.

#include <TGenPhaseSpace.h>
#include <TLorentzVector.h>
#include <TH2F.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TMath.h>
#include <TLine.h>
#include <iostream>

double bw_shape(double m, double m0, double gamma)
{
    double dm2 = m*m - m0*m0;
    return 1.0 / (dm2*dm2 + m0*m0*gamma*gamma);
}

void dalitz_plot_with_sigma_bands()
{
    gStyle->SetOptStat(0);
    gROOT->GetListOfCanvases()->Delete();

    const double m_lambda    = 1.115683;
    const double m_pion      = 0.139570;
    const double m_sigma     = 1.387;
    const double gamma_sigma = 0.036;
    const double sqrt_s      = 1.95;
    const int n_events       = 300000;

    TLorentzVector p_initial(0, 0, 0, sqrt_s);

    double masses[3] = {m_lambda, m_pion, m_pion};
    TGenPhaseSpace phase_space;

    if (!phase_space.SetDecay(p_initial, 3, masses)) {
        std::cerr << "Error: phase-space setup failed." << std::endl;
        return;
    }

    TH2F *h_dalitz_dyn = new TH2F(
        "h_dalitz_dyn",
        "Dalitz Plot with #Sigma^{#pm}(1385) enhancement; (m_{#Lambda#pi^{-}})^{2} (GeV^{2}/c^{4}); (m_{#Lambda#pi^{+}})^{2} (GeV^{2}/c^{4})",
        300, 1.5, 3.5,
        300, 1.5, 3.5
    );

    double bw_norm = bw_shape(m_sigma, m_sigma, gamma_sigma);

    for (int i = 0; i < n_events; ++i) {
        double w_ps = phase_space.Generate();
        if (w_ps <= 0.0) continue;

        TLorentzVector p_lambda   = *phase_space.GetDecay(0);
        TLorentzVector p_pi_plus  = *phase_space.GetDecay(1);
        TLorentzVector p_pi_minus = *phase_space.GetDecay(2);

        double m12 = (p_lambda + p_pi_minus).M();
        double m13 = (p_lambda + p_pi_plus ).M();

        double m12_sq = m12 * m12;
        double m13_sq = m13 * m13;

        double w_dyn = 1.0
                     + 2.0 * bw_shape(m12, m_sigma, gamma_sigma) / bw_norm
                     + 2.0 * bw_shape(m13, m_sigma, gamma_sigma) / bw_norm;

        h_dalitz_dyn->Fill(m12_sq, m13_sq, w_ps * w_dyn);
    }

    TCanvas *c2 = new TCanvas("c2", "Dalitz Plot: Sigma bands", 800, 650);
    c2->SetLogz();
    h_dalitz_dyn->Draw("COLZ");

    double m_sigma_sq = m_sigma * m_sigma;

    TLine *line_minus = new TLine(m_sigma_sq, 1.5, m_sigma_sq, 3.5);
    line_minus->SetLineColor(kRed);
    line_minus->SetLineStyle(2);
    line_minus->Draw();

    TLine *line_plus = new TLine(1.5, m_sigma_sq, 3.5, m_sigma_sq);
    line_plus->SetLineColor(kBlue);
    line_plus->SetLineStyle(2);
    line_plus->Draw();
    c2->Draw();
}

dalitz_plot_with_sigma_bands();

Interpretation of Step 2

Compared with the pure phase-space plot, this Dalitz plot shows two enhanced bands:

• a vertical band from $\Sigma^-(1385)\to\Lambda\pi^-$,
• a horizontal band from $\Sigma^+(1385)\to\Lambda\pi^+$.

This makes the physical meaning of a Dalitz plot very clear:

• the boundary is determined by kinematics;
• the population inside the boundary is determined by dynamics;
• resonant substructures appear as visible bands.