Invariant Mass

Invariant mass is a Lorentz-invariant quantity that characterizes the total energy and momentum of a particle or system, remaining constant across all reference frames. It’s a cornerstone for analyzing particle collisions and decays.

The invariant mass for two particles with four-momenta $ P_1 = (E_1, \vec{p}_1) $ and $ P_2 = (E_2, \vec{p}_2) $ is defined as:

$$ M^2_{Inv} = (P_1 + P_2)^2 = (E_1 + E_2)^2 - |\vec{p}_1 + \vec{p}_2|^2 $$

Here, $ E_i $ is the energy, $ \vec{p}_i $ is the three-momentum, and units assume $ c = 1 $. For a system of $ n $ particles, this generalizes to:

$$ M(n)^2_{Inv} = (P_1 + P_2 + \dots + P_n)^2 $$

For a single particle, the invariant mass equals its rest mass, $ m_0 $, since in its rest frame, $ P = (m_0, 0, 0, 0) $, yielding:

$$ M^2_{Inv} = m_0^2 $$

This invariance makes it ideal for identifying particles or resonances, regardless of the observer’s frame.

Intermediate States and Resonances

In collision experiments, short-lived intermediate particles or resonances often decay into secondary particles. To confirm their existence, physicists measure the four-momenta of decay products and compute the invariant mass of specific combinations. A peak in the invariant mass distribution indicates a resonance. For example, in $ Z \to e^+ + e^- $, the invariant mass of the electron-positron pair peaks around 91 GeV, confirming the Z boson.

Two-Body Collision Reactions

Consider a reaction $ A + B \to \text{products} $. Two processes highlight distinct invariant mass behaviors:

1. Process (a): $ A + B \to c + d + e $
   This is a 2 → 3 scattering process, producing three particles directly. The invariant mass of a pair, e.g., $ c + d $, follows a statistical distribution governed by phase space, which defines kinematic constraints (energy-momentum conservation). This results in a broad, continuous spectrum:

   $$ M^2_{c+d} \sim \text{phase space distribution} $$

2. Process (b): $ A + B \to c + d, d \to e + f $
   Here, $ c $ and an intermediate particle $ d $ are produced, with $ d $ decaying into $ e + f $. The invariant mass of $ e + f $ follows a Breit-Wigner distribution, centered around the mass of $ d $, with a width determined by its lifetime:

   $$ M^2_{e+f} \sim \text{Breit-Wigner}(m_d, \Gamma_d) $$

   For example, if $ d $ is a rho meson ($ \rho \to \pi^+ + \pi^- $), the pion pair’s invariant mass distribution peaks around 770 MeV with a Breit-Wigner shape.

Missing Mass

Missing mass complements invariant mass in experiments where particles (e.g., neutrinos) escape detection. By measuring the four-momenta of detected particles and knowing the initial state (e.g., $ P_A + P_B $), the missing four-momentum $ P_{\text{miss}} $ is calculated as:

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

The invariant mass squared of this missing four-momentum gives the missing mass:

$$ M^2_{\text{miss}} = P_{\text{miss}}^2 $$

This quantity reveals properties of undetected particles, such as their mass, and is consistent with the total invariant mass of the system.

Example: Detecting Short-Lived Resonances

• Particles liken $\rho,\omega,\phi$ have an average lifetime of 10$^{-22}$-10$^{-23}$s

  • How do we know of their existence if they live so shortly?

• reaction pp$\to$ pp$\pi^+\pi^-$

  • We identify all four particles in the final state and measure their momentum
  • Let's focus on the pion pair, the total energy & momentum are: $E=E_{\pi^+}+E_{\pi^-} \space\space\space\space P=P_{\pi^+}+P_{\pi^-}$
  • The invariant mass is: $M=(E^2- P^2)^{1/2}$
  • The event distribution for the variable M from $\pi^+\pi^-$ will look like

  

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

// 全局变量
TCanvas *c1 = new TCanvas("c1", "Invariant Mass", 600, 400);
TH1F *h_total = new TH1F("h_total", "#pi^{+}#pi^{-} Invariant Mass;M_{inv} (GeV/c^{2});Events", 100, 0.0, 3.0);
TH1F *h_rho = new TH1F("h_rho", "Rho Contribution", 100, 0.0, 3.0);
TH1F *h_nonrho = new TH1F("h_nonrho", "Non-Rho Contribution", 100, 0.0, 3.0);

void invariant_mass_demo() {
    // 随机数生成器
    TRandom3 rand(0);

    // 物理常数
    const double m_pion = 0.13957;    // 介子质量 (GeV/c^2)
    const double m_proton = 0.93827;  // 质子质量 (GeV/c^2)
    const double rho_mass = 0.775;    // ρ 介子质量 (GeV/c^2)
    const double rho_width = 0.149;   // ρ 介子宽度 (GeV/c^2)
    const double sqrt_s = 7.0;        // 质心系能量 (GeV)
    const int n_events = 100000;      // 事件数
    const double frac_rho = 0.1;      // ρ 过程的比例

    // 初始态：pp 碰撞
    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;

    // 相空间设置
    // 1. pp -> pp rho
    TGenPhaseSpace phaseSpace_pprho;
    double masses_pprho[3] = {m_proton, m_proton, rho_mass}; // 默认 ρ 质量，稍后调整

    // 2. rho -> pi+ pi-
    TGenPhaseSpace phaseSpace_rho;
    double masses_rho[2] = {m_pion, m_pion};

    // 3. pp -> pp pi+ pi- (四体衰变)
    TGenPhaseSpace phaseSpace_full;
    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;
        TLorentzVector p4_pion1, p4_pion2;

        if (is_rho) {
            // 共振过程：pp -> pp rho, rho -> pi+ pi-
            double M_rho = rand.BreitWigner(rho_mass, rho_width);
            masses_pprho[2] = M_rho; // 更新 ρ 质量
            phaseSpace_pprho.SetDecay(p4_initial, 3, masses_pprho);
            double weight1 = phaseSpace_pprho.Generate();

            // 获取 ρ 的四动量
            TLorentzVector p4_rho = *phaseSpace_pprho.GetDecay(2);

            // ρ 衰变成 pi+ pi-
            phaseSpace_rho.SetDecay(p4_rho, 2, masses_rho);
            double weight2 = phaseSpace_rho.Generate();
            p4_pion1 = *phaseSpace_rho.GetDecay(0);
            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 {
            // 非共振过程：pp -> pp pi+ pi-
            phaseSpace_full.SetDecay(p4_initial, 4, masses_full);
            double weight = phaseSpace_full.Generate();
            p4_pion1 = *phaseSpace_full.GetDecay(2); // pi+
            p4_pion2 = *phaseSpace_full.GetDecay(3); // pi-

            // 计算不变质量
            double inv_mass = (p4_pion1 + p4_pion2).M();
            h_nonrho->Fill(inv_mass, weight);
            h_total->Fill(inv_mass, weight);
        }
    }
}

invariant_mass_demo();

// 绘图
c1->Clear();
h_total->SetLineColor(kBlack);
h_total->SetFillColor(kGray+1);
h_total->SetFillStyle(1001);
h_total->Draw();

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

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

// 图例
TLegend *leg = new TLegend(0.6, 0.7, 0.9, 0.9);
leg->AddEntry(h_total, "Total", "f");
leg->AddEntry(h_rho, "#rho #rightarrow #pi^{+}#pi^{-}", "f");
leg->AddEntry(h_nonrho, "Non-resonant", "f");
leg->Draw();

h_total->GetXaxis()->SetTitleSize(0.05);
h_total->GetYaxis()->SetTitleSize(0.05);

c1->Draw();

Three-Body Decays and the Dalitz Plot

A three-body decay occurs when a particle of mass $ M $ decays into three particles with rest masses $ m_1 $, $ m_2 $, $ m_3 $:

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

Kinematics and Degrees of Freedom

Each decay product has a four-momentum $ P_i = (E_i, \vec{p}_i) $, totaling 12 parameters (4 per particle). Constraints reduce the free parameters:

• Momentum Conservation: In the parent’s rest frame, $ \vec{p}_1 + \vec{p}_2 + \vec{p}_3 = 0 $, eliminating 3 parameters (9 remain).
• Energy Conservation: $ E_1 + E_2 + E_3 = M $, removing 1 parameter (8 remain).
• On-Shell Conditions: Each particle satisfies $ P_i^2 = m_i^2 $, reducing 3 more parameters (5 remain).
• Rest Frame Geometry: The decay products lie in a plane. Fixing the plane’s orientation (2 angles) leaves 3 parameters. Choosing one particle’s direction in the plane (1 angle) reduces to 2 degrees of freedom.

Thus, a three-body decay is fully described by two kinematic variables, shaped by phase space, which dictates the allowed energy and momentum configurations.

Dalitz Plot Representation

The Dalitz plot visualizes three-body decays by plotting two kinematic variables, revealing decay dynamics and resonances. Common variables include the squared invariant masses of particle pairs:

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

These are related by:

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

Since $ P_1 + P_2 + P_3 = P_M $ and $ P_M^2 = M^2 $ in the parent’s rest frame, only two of $ m_{12}^2 $, $ m_{13}^2 $, $ m_{23}^2 $ are independent. The Dalitz plot typically plots two (e.g., $ m_{12}^2 $ vs. $ m_{23}^2 $), with phase space defining the allowed region, a bounded area determined by energy-momentum conservation. Resonances appear as bands when a pair’s invariant mass (e.g., $ m_{23} $) aligns with an intermediate state’s mass.

Alternatively, Dalitz’s original variables use kinetic energies $ T_i = E_i - m_i $:

$$ x = \sqrt{3} \frac{T_1 - T_2}{Q}, \quad y = \frac{2 T_3 - T_1 - T_2}{Q} $$

where $ Q = M - m_1 - m_2 - m_3 $ is the decay’s available energy. These map the same phase space but emphasize energy differences.

Phase Space and Dynamics

Phase space uniformly populates the Dalitz plot for non-resonant decays, producing a flat distribution within kinematic boundaries. Resonances or interactions (e.g., particle 2 and 3 forming a state of mass $ m_{23} \approx m_R $) cause enhancements, visible as bands or peaks. The plot’s density reflects the decay matrix element, modulated by phase space constraints, making it a powerful tool for identifying intermediate states.

Example:

• As an example, let's study the reaction:
  • $K^{-}p\to\Lambda\pi^{+}\pi^{-}(\Lambda\to\pi^{-}p)$
• We can measure two invariant masses
  • $m_{12}\equiv m(\Lambda\pi^{-})_{1}\quad m_{13}\equiv m(\Lambda\pi^{+})$
• The so-called“Dalitz plot”shows the relation between $({m}_{13})^{2}\mathrm{~and~(m}_{12})^{2}$

• The $\Sigma^{\pm}$ resonance appears as two bands in the Dalitz plot around 1.4 GeV
  • $\Sigma^{\pm}(1387)\to \Lambda \pi^{\pm}$

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

// Particle masses (GeV/c^2)
const double m_lambda = 1.115683;
const double m_pion = 0.139570;
const double m_kaon = 0.493677;
const double m_proton = 0.938272;
const double m_sigma = 1.387;
const double gamma_sigma = 0.036;

// Number of events
const int n_events = 500000;

// Center-of-mass energy (GeV)
const double sqrt_s = 1.95;

// Create canvas and draw
TCanvas *c2 = new TCanvas("c2", "Dalitz Plot", 800, 600);
c2->SetLogz();  // Log scale for visibility

// Create histogram for Dalitz plot
TH2F *h_dalitz = new TH2F("h_dalitz", "Dalitz Plot for K^{-} p #rightarrow #Lambda #pi^{+} #pi^{-}",
                            400, 1.5, 3.5, 400, 1.5, 3.5);

double breit_wigner(double m, double m0, double gamma) {
    // Relativistic Breit-Wigner probability density
    double m2 = m * m;
    double m02 = m0 * m0;
    double gamma2 = gamma * gamma;
    double denom = (m2 - m02) * (m2 - m02) + m02 * gamma2;
    return (m0 * gamma / TMath::Pi()) / denom;
}

void dalitz_plot() {
    TRandom3 rand(42);

    // Define total CM four-momentum
    TLorentzVector p_initial(0, 0, 0, sqrt_s);

    // Setup phase space generator
    TGenPhaseSpace phase_space;
    double masses[3] = {m_lambda, m_pion, m_pion};
    if (!phase_space.SetDecay(p_initial, 3, masses)) {
        std::cerr << "Error: Phase space setup failed!" << std::endl;
        return;
    }
    int n_accepted = 0;
    for (int i = 0; i < n_events; ++i) {
        double weight = phase_space.Generate();
        if (weight <= 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);

        TLorentzVector p12 = *p_lambda + *p_pi_minus;
        TLorentzVector p13 = *p_lambda + *p_pi_plus;
        double m12 = p12.M();
        double m13 = p13.M();

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

        double bw_plus = breit_wigner(m13, m_sigma, gamma_sigma);
        double bw_minus = breit_wigner(m12, m_sigma, gamma_sigma);
        double weight_res = 0.5 * (bw_plus + bw_minus);
        double max_weight = breit_wigner(m_sigma, m_sigma, gamma_sigma);
        weight_res = TMath::Min(weight_res / max_weight + 0.15, 1.0);

        if (rand.Uniform() < weight_res) {
            h_dalitz->Fill(m12_sq, m13_sq, weight);
            n_accepted++;
        }
    }


}

dalitz_plot();

c2->Clear();
h_dalitz->Draw("COLZ");

h_dalitz->SetXTitle("(m_{#Lambda #pi^{-}})^{2} (GeV^{2}/c^{4})");
h_dalitz->SetYTitle("(m_{#Lambda #pi^{+}})^{2} (GeV^{2}/c^{4})");
h_dalitz->SetZTitle("Event Density");

// Add legend
TLegend *leg = new TLegend(0.7, 0.8, 0.9, 0.9);
leg->AddEntry(line_minus, "#Sigma^{-}(1387)", "l");
leg->AddEntry(line_plus, "#Sigma^{+}(1387)", "l");
leg->Draw();

// Add resonance lines    
double m_sigma_sq = m_sigma * m_sigma;
TLine *line_minus = new TLine(m_sigma_sq, 1.0, m_sigma_sq, 4.0);
line_minus->SetLineColor(kRed);
line_minus->SetLineStyle(2);
    line_minus->Draw();

TLine *line_plus = new TLine(1.0, m_sigma_sq, 4.0, m_sigma_sq);
line_plus->SetLineColor(kBlue);
line_plus->SetLineStyle(2);
line_plus->Draw();
c2->Draw();