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

反应 $Q$ 值重建

反应 $Q$ 值基本定义

对于核反应

$$ A(a,b)B, $$

反应 $Q$ 值定义为末态与初态动能之差：

$$ Q = E_{KB}+E_{Kb}-E_{KA}-E_{Ka}. \tag{1} $$

在相对论形式下，粒子 $i$ 的总能量写为

$$ E_i = E_{Ki}+m_i, \tag{2} $$

其中 $E_{Ki}$ 为动能，$m_i$ 为静质量。在自然单位 $c=1$ 下，将式(2)代入式(1)，可得

$$ Q = m_A+m_a-m_B-m_b. \tag{3} $$

若粒子 $i$ 处于激发态，其质量写为

$$ m_i = m_{i0}+E_{x,i}, $$

其中 $m_{i0}$ 为基态质量，$E_{x,i}$ 为激发能。于是，当所有初末态粒子均处于基态时，基态 $Q$ 值为

$$ Q_{gs}=m_{A0}+m_{a0}-m_{B0}-m_{b0}. \tag{4} $$

若末态核 $B$ 处于激发态，则有

$$ Q = Q_{gs}-E_{x,B}. \tag{5} $$

式 (5) 使用的是退激前末态核的质量。如果核已发射未探测的 γ 射线，带电粒子动能之和缺少的是实验室系 γ 能量。后文会区分这一精确关系与按激发能识别分支的近似。

---

$^{14}\mathrm C$ 束流在 $(\mathrm{CH})_n$ 靶上的反应：

$$ {}^{1}\mathrm H({}^{14}\mathrm C,{}^{14}\mathrm C^{*}){}^{1}\mathrm H, $$

和

$$ {}^{12}\mathrm C({}^{14}\mathrm C,{}^{14}\mathrm C^{*}){}^{12}\mathrm C. $$

两类反应都采用相同的级联两体衰变图像。入射 $^{14}\mathrm C$ 束流首先与靶核作用，产生反冲粒子和激发态 $^{14}\mathrm C^{*}$；随后

$$ {}^{14}\mathrm C^{*}\rightarrow {}^{4}\mathrm{He}+{}^{10}\mathrm{Be}^{*}, $$

再进一步

$$ {}^{10}\mathrm{Be}^{*}\rightarrow {}^{10}\mathrm{Be}+\gamma. $$

对 H、C 两类过程，反冲核的初末态质量相同，因此基态反应 Q 值相同： $$Q_{\mathrm{gs}}=m_{14\mathrm C}-m_\alpha-m_{10\mathrm{Be}}\simeq-12.01\ \mathrm{MeV}.$$ 若 γ 未探测，真值带电动能求和严格满足 $$Q_{\mathrm{charged}}=T_\alpha+T_{\mathrm{Be}}+T_R-T_{\mathrm{beam}}=Q_{\mathrm{gs}}-E_{\gamma,\mathrm{lab}}.$$ 实验室系 $E_\gamma$ 随发射方向发生 Doppler 变化。因此约 3.368 MeV 的 Q 峰间隔可用来识别分支，但逐事件缺失能量并不严格等于 3.368 MeV。

---

$Q$ 值测量对于不变质量谱（Invariant Mass）的意义

母核激发能与衰变碎片的相对能、质量阈值有关。只测量 α 和 γ 退激后的 $^{10}\mathrm{Be}$ 时，$(P_\alpha+P_{\mathrm{Be}})^2$ 是测到的两粒子系统的不变质量，并非完整母核质量；完整四动量还包含未测量的 $P_\gamma$。

本例用 Q 谱区分 $^{10}\mathrm{Be}$ 基态和 3.368 MeV 分支，再选择相应的重建处理。有限分辨使分支有重叠，门选只能提高分支纯度，不能保证每个事件都分类正确。Q 值在本例是重要的分支判据，不是不变质量分析在所有情形下都必须先做的一步。

$Q$ 值重建方法

反应初态

若没有逐事件束流能量测量，可采用束流调谐或独立测量给出的平均值。已知电荷态 $q$ 时，磁刚度 $B\rho=p/q$ 给出动量，再由 $T=\sqrt{p^2+m^2}-m$ 换算动能。是否能忽略束流能散取决于具体束流和测量精度，不能仅按稳定束流或放射性束流分类。本例沿用 4.4 的 317.1 MeV 均值和 1.83 MeV FWHM。

对于放射性束流，则通常有$\Delta E/E\sim 1\%$。 这会直接影响 $Q$ 值分辨。

反应末态

反应中，没有测量$\gamma$信息，因此实验上可用于 $Q$ 值重建的带电末态粒子为

$$ {}\mathrm H/{}^{12}\mathrm C,\quad {}^{4}\mathrm{He},\quad {}^{10}\mathrm{Be}, $$

其可测量量为动能 $E_{Ki}$ 和动量 $\vec p_i$。

（1）利用束流信息直接重建

若束流逐事件能量已知，则可定义

$$ Q_{\mathrm{event}} = E_{k,\alpha}+E_{k,^{10}\mathrm{Be}}+E_{k,\mathrm{recoil}}-E_{k,\mathrm{beam}}^{\mathrm{event}} . $$

若束流逐事件能量不可得，而仅知道平均束流能量 $\bar E_{k,\mathrm{beam}}$，则可定义

$$ Q_{\mathrm{mean}} = E_{k,\alpha}+E_{k,^{10}\mathrm{Be}}+E_{k,\mathrm{recoil}}-\bar E_{k,\mathrm{beam}} . $$ 其中 $Q_{\mathrm{mean}}$ 更接近实验中仅由磁刚度 $B\rho$ 给出平均束流能量时的情况，其分辨通常受束流能散影响更明显。

（2）利用三体末态动量守恒重建束流能量（3-body）

若三种带电末态粒子 $\alpha$、$^{10}\mathrm{Be}$、recoil 均被测得，则可先由它们的动能与方向重建三维动量

$$ \vec p_i = \sqrt{E_{k,i}^2+2m_iE_{k,i}}\ \hat r_i . $$

末态总动量为

$$ \vec p_{\mathrm{out}}=\vec p_{\alpha}+\vec p_{^{10}\mathrm{Be}}+\vec p_{\mathrm{recoil}} . $$

没有漏掉其他末态粒子时，$\vec p_{\mathrm{out}}$ 等于入射动量；本例对未测 γ 分支忽略 $\vec p_\gamma$，所以这里是近似。由带电末态动量和估计束流动能：

$$ E'_{k,\mathrm{beam}} = \sqrt{m_{^{14}\mathrm C}^2+p_{\mathrm{out}}^2}-m_{^{14}\mathrm C}. $$

于是得到逐事件的三体 Q 值重建：

$$ Q_{3\mathrm{body}} = E_{k,\alpha}+E_{k,^{10}\mathrm{Be}}+E_{k,\mathrm{recoil}}-E'_{k,\mathrm{beam}} . $$

这种方法不用束流平均能量，但会传播三个粒子的能量和角度误差；未探测的 γ 动量也会引入偏差。因此它能否改善 Q 分辨，需要用同一组数据比较，不能只由公式保证。

（3）反冲粒子缺失时的二体重建（2-body）

若 recoil 未被测得，但 $\alpha$ 与 $^{10}\mathrm{Be}$ 被测得，同时已知平均束流能量 $\bar E_{k,\mathrm{beam}}$，则可先构造平均束流动量 $\vec p_{\mathrm{beam}}^{\,\mathrm{mean}}$，再由动量守恒得到 recoil 动量：

$$ \vec p_{\mathrm{recoil}}^{\,\mathrm{rec}} = \vec p_{\mathrm{beam}}^{\,\mathrm{mean}} - \vec p_{\alpha} - \vec p_{^{10}\mathrm{Be}} . $$

由此可求 recoil 的重建动能

$$ E_{k,\mathrm{recoil}}^{\,\mathrm{rec}} = \sqrt{m_{\mathrm{recoil}}^2+\left(p_{\mathrm{recoil}}^{\mathrm{rec}}\right)^2} - m_{\mathrm{recoil}} . $$

于是得到

$$ Q_{2\mathrm{body}} = E_{k,\alpha}+E_{k,^{10}\mathrm{Be}}+E_{k,\mathrm{recoil}}^{\,\mathrm{rec}} - \bar E_{k,\mathrm{beam}} . $$

该方法适用于 recoil 探测效率较低或实验几何不覆盖 recoil 的情况。

---

由 Q 值反推出 $^{10}\mathrm{Be}$ 激发能

可定义分支识别量 $$E_{x,\mathrm{est}}(^{10}\mathrm{Be})=Q_{\mathrm{gs}}-Q_{\mathrm{rec}}.$$ 基态事件集中在零附近，第一激发态事件集中在约 3.368 MeV 附近。由于 γ 的 Doppler 变化、漏掉的 γ 动量及测量误差，这不是逐事件精确的激发能测量。

下面保留 Qevt、Qmean、Q3、Q2 四种计算，使用同一输入样本比较。ekBeam 和 exBe10 仅作真值对照；图中分开的 H/C 曲线使用 processID 判断反冲核质量，属于已知真值类别的性能检查。处理实验中的未知类别时，要分别提出 H/C 假设或采用 4.6 的关联选择，不能读取这个真值标签。

代码中 Kine 只是把公式放在一起：P 计算动量大小，pv 根据角度构造三维动量，lv 加入总能量；Qevt/Qmean/Q3/Q2 分别实现上述四个公式，ExInv 计算未补 γ 的两碎片不变质量。SetBranchAddress 将当前事件数据读入数组，GetEntry(i) 后才计算第 $i$ 个事件。

输入是 4.4 生成的 C14_CHn_He4Be10x.root。文件里的质量参数是 GeV，统一乘 1000 变成 MeV；树里的动能已经是 MeV，不再转换。

%jsroot on

%%cpp -d
#include <cmath>
#include <algorithm>
#include "TFile.h"
#include "TTree.h"
#include "TParameter.h"
#include "TMath.h"
#include "TVector3.h"
#include "TLorentzVector.h"
#include "TCanvas.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TLegend.h"
#include "TLine.h"
#include "TStyle.h"

// 以下两个画布在后面的显示单元中继续使用。
TCanvas *c1 = nullptr;
TCanvas *c2 = nullptr;

struct Kine {
    double mH1, mC12, mHe4, mBe10, mC14, ekBeamMean;

    double mRecoil(int pid) const { return pid == 1 ? mH1 : mC12; }
    double Qgs(int pid) const { return mC14 - mHe4 - mBe10; }
    double MomentumMagnitude(double ek, double m) const { return std::sqrt(ek*ek + 2*m*ek); }

    double ParticleMass(int pid, int i, double ex) const {
        double mass[] = {mHe4, mBe10, mRecoil(pid), mBe10 + ex};
        return mass[i];
    }

    TVector3 MomentumVector(int pid, int i, double ek[], double th[], double ph[], double ex) const {
        TVector3 v;
        v.SetMagThetaPhi(MomentumMagnitude(ek[i], ParticleMass(pid,i,ex)), th[i]*TMath::DegToRad(), ph[i]*TMath::DegToRad());
        return v;
    }

    TLorentzVector FourMomentum(int pid, int i, double ek[], double th[], double ph[], double ex) const {
        return TLorentzVector(MomentumVector(pid,i,ek,th,ph,ex), ek[i] + ParticleMass(pid,i,ex));
    }

    double Qevt(double ek[], double ekBeam) const { return ek[0]+ek[1]+ek[2] - ekBeam; }
    double Qmean(double ek[]) const { return ek[0]+ek[1]+ek[2] - ekBeamMean; }

    double Q3(int pid, double ek[], double th[], double ph[], double ex) const {
        TVector3 p = MomentumVector(pid,0,ek,th,ph,ex) + MomentumVector(pid,1,ek,th,ph,ex) + MomentumVector(pid,2,ek,th,ph,ex);
        return ek[0]+ek[1]+ek[2] - (std::sqrt(mC14*mC14 + p.Mag2()) - mC14);
    }

    double Q2(int pid, double ek[], double th[], double ph[], double ex) const {
        TVector3 pR = TVector3(0,0,MomentumMagnitude(ekBeamMean,mC14)) - MomentumVector(pid,0,ek,th,ph,ex) - MomentumVector(pid,1,ek,th,ph,ex);
        double mR = mRecoil(pid);
        return ek[0]+ek[1] + (std::sqrt(mR*mR + pR.Mag2()) - mR) - ekBeamMean;
    }

    double ExInv(int pid, double ek[], double th[], double ph[], double ex) const {
        return (FourMomentum(pid,0,ek,th,ph,ex) + FourMomentum(pid,1,ek,th,ph,ex)).M() - mC14;
    }
};

int colH = kBlue+1, colC = kRed+1;
void vline(double x, double y1, double y2, int col, int sty=2) {
    TLine *l = new TLine(x,y1,x,y2);
    l->SetLineColor(col); l->SetLineStyle(sty); l->Draw();
}

读取质量与束流参数，并连接 Tree 的各个 branch。质量从 GeV 转为 MeV 后，所有重建量统一用 MeV。

TFile *f = TFile::Open("C14_CHn_He4Be10x.root");
if (!f || f->IsZombie()) throw std::runtime_error("Run Section 4.4 first");
TTree *t = f->Get<TTree>("tree");
if (!t) throw std::runtime_error("Missing tree");

Kine K;
const double u = 1000.0; // 文件中质量为 GeV，转换为 MeV；动能分支已经是 MeV
K.mH1        = ((TParameter<Double_t>*)f->Get("massH1"))->GetVal() * u;
K.mC12       = ((TParameter<Double_t>*)f->Get("massC12"))->GetVal() * u;
K.mHe4       = ((TParameter<Double_t>*)f->Get("massHe4"))->GetVal() * u;
K.mBe10      = ((TParameter<Double_t>*)f->Get("massBe10"))->GetVal() * u;
K.mC14       = ((TParameter<Double_t>*)f->Get("massC14"))->GetVal() * u;
K.ekBeamMean = ((TParameter<Double_t>*)f->Get("ekBeamMean"))->GetVal();

Int_t pid;
Double_t ekBeam, exBe10, w;
Double_t ek[4], th[4], ph[4];
t->SetBranchAddress("processID",    &pid);
t->SetBranchAddress("ekBeam",       &ekBeam);
t->SetBranchAddress("exBe10",       &exBe10);
t->SetBranchAddress("weight_total", &w);
t->SetBranchAddress("ek",           ek);
t->SetBranchAddress("theta",        th);
t->SetBranchAddress("phi",          ph);

为各个重建方法分别准备直方图；第二个下标 0、1 对应 H 靶和 C 靶。

// 直方图: [0]=H, [1]=C
TH1F *hQ[4][2], *hEx[3][2];
TH2F *hInvQ[2], *hInvQ2[2];
const char *qn[] = {"Qevt","Qmean","Q3","Q2"};
for (int i = 0; i < 4; ++i) {
    hQ[i][0] = new TH1F(Form("%s_H",qn[i]), Form("%s;Q [MeV];counts",qn[i]), 500,-35,5);
    hQ[i][1] = new TH1F(Form("%s_C",qn[i]), Form("%s;Q [MeV];counts",qn[i]), 500,-35,5);
}
const char *en[] = {"ExTrue","ExQ3","ExQ2"};
for (int i = 0; i < 3; ++i) {
    hEx[i][0] = new TH1F(Form("%s_H",en[i]), "H: E_{x}(^{10}Be);E_{x} [MeV];counts", 300,-5,10);
    hEx[i][1] = new TH1F(Form("%s_C",en[i]), "C: E_{x}(^{10}Be);E_{x} [MeV];counts", 300,-5,10);
}
hInvQ[0]  = new TH2F("InvQ_H",  "H: inv.mass vs Q3;Q3 [MeV];E_{x}(^{14}C) [MeV]", 200,-20,5, 200,8,22);
hInvQ[1]  = new TH2F("InvQ_C",  "C: inv.mass vs Q3;Q3 [MeV];E_{x}(^{14}C) [MeV]", 200,-20,5, 200,8,22);
hInvQ2[0] = new TH2F("InvQ2_H", "H: inv.mass vs Q2;Q2 [MeV];E_{x}(^{14}C) [MeV]", 200,-20,5, 200,8,22);
hInvQ2[1] = new TH2F("InvQ2_C", "C: inv.mass vs Q2;Q2 [MeV];E_{x}(^{14}C) [MeV]", 200,-20,5, 200,8,22);

每个事件只计算一次 Q3、Q2 和不变质量，再用相同权重填入各幅对比图。

// 填充
for (Long64_t i = 0; i < t->GetEntriesFast(); ++i) {
    t->GetEntry(i);
    int j = (pid == 1) ? 0 : 1;

    double q3 = K.Q3(pid, ek, th, ph, exBe10);
    double q2 = K.Q2(pid, ek, th, ph, exBe10);
    double exInv = K.ExInv(pid, ek, th, ph, exBe10);
    hQ[0][j]->Fill(K.Qevt(ek, ekBeam), w);
    hQ[1][j]->Fill(K.Qmean(ek), w);
    hQ[2][j]->Fill(q3, w);
    hQ[3][j]->Fill(q2, w);

    hEx[0][j]->Fill(exBe10, w);
    hEx[1][j]->Fill(K.Qgs(pid) - q3, w);
    hEx[2][j]->Fill(K.Qgs(pid) - q2, w);

    hInvQ[j]->Fill(q3, exInv, w);
    hInvQ2[j]->Fill(q2, exInv, w);
}

std::cout << "Qgs = " << K.Qgs(1) << " MeV; events = " << t->GetEntries() << std::endl;

Qgs = -12.0125 MeV; events = 1985670

// 画图
gStyle->SetOptStat(0);

// ================ c1: Q 值谱对比 (2x2, logy) ================
c1 = new TCanvas("c1","Q",900,900);
c1->Divide(2,2);
for (int i = 0; i < 4; ++i) {
    c1->cd(i+1); gPad->SetLogy();
    hQ[i][0]->SetLineColor(colH);
    hQ[i][1]->SetLineColor(colC);
    double ymax = 2*std::max(hQ[i][0]->GetMaximum(), hQ[i][1]->GetMaximum());
    hQ[i][0]->SetMinimum(0.5); hQ[i][0]->SetMaximum(ymax);
    hQ[i][0]->Draw("hist"); hQ[i][1]->Draw("hist same");
    vline(K.Qgs(1), 0.5, ymax, colH); vline(K.Qgs(1)-3.368, 0.5, ymax, colH, 3);
    vline(K.Qgs(2), 0.5, ymax, colC); vline(K.Qgs(2)-3.368, 0.5, ymax, colC, 3);
    if (i == 0) {
        TLegend *leg = new TLegend(0.15,0.75,0.45,0.88);
        leg->AddEntry(hQ[i][0], "H target", "l");
        leg->AddEntry(hQ[i][1], "C target", "l");
        leg->Draw();
    }
}

// ================ c2: 激发能验证 + 不变质量关联 (3x2) ================

c1->Draw();

不同 Q 值重建方法的对比

四个面板使用同一组事件：

• Qevt：使用真实逐事件束流动能，作为去掉束流能散这一不确定来源后的对照；仍含读出分辨和未探测 γ 的影响。
• Qmean：换成平均束流能量，展示忽略逐事件能量变化带来的展宽。
• Q3：由三个带电粒子的动量和估计束流动能，比较改善了多少，以及哪些测量误差仍在。
• Q2：只使用两个前向碎片，并给反冲核指定质量假设；平均束流能量的影响仍保留。

H/C 两类曲线的基态 Q 值相同，但动量误差的传播不同，特别是反冲粒子质量、动能和方向的差别会改变重建分辨。此样本没有模拟厚靶能损或能量歧离，不能用这些效应解释图中的差别。图中虚线仅标示 $Q_{\mathrm{gs}}$ 和 $Q_{\mathrm{gs}}-3.368$ MeV 的分支参考位置。

c2 = new TCanvas("c2","Ex",900,1200);
c2->Divide(2,3);

// --- 第一行: E_x(10Be) 重建 (logy) ---
for (int j = 0; j < 2; ++j) {
    c2->cd(j+1); gPad->SetLogy();
    hEx[0][j]->SetLineColor(kBlack);
    hEx[1][j]->SetLineColor(kRed+1);
    hEx[2][j]->SetLineColor(kBlue+1); hEx[2][j]->SetLineStyle(2);
    double ymax = 2.0*std::max({hEx[0][j]->GetMaximum(), hEx[1][j]->GetMaximum(), hEx[2][j]->GetMaximum()});
    hEx[0][j]->SetMinimum(0.5); hEx[0][j]->SetMaximum(ymax);
    hEx[0][j]->Draw("hist"); hEx[1][j]->Draw("hist same"); hEx[2][j]->Draw("hist same");
    TLegend *leg = new TLegend(0.55,0.7,0.88,0.88);
    leg->AddEntry(hEx[0][j], "true", "l");
    leg->AddEntry(hEx[1][j], "from Q3", "l");
    leg->AddEntry(hEx[2][j], "from Q2", "l");
    leg->Draw();
}

// --- 第二行: inv.mass vs Q2 ---
for (int j = 0; j < 2; ++j) {
    c2->cd(j+3); gPad->SetRightMargin(0.12);
    hInvQ2[j]->Draw("colz");
    int p = j+1;
    vline(K.Qgs(p), 8, 22, kWhite); vline(K.Qgs(p)-3.368, 8, 22, kMagenta+1);
}
// --- 第三行: inv.mass vs Q3 ---
for (int j = 0; j < 2; ++j) {
    c2->cd(j+5); gPad->SetRightMargin(0.12);
    hInvQ[j]->Draw("colz");
    int p = j+1;
    vline(K.Qgs(p), 8, 22, kWhite); vline(K.Qgs(p)-3.368, 8, 22, kMagenta+1);
}
c2->Draw();

---

重建激发能与不变质量

第一行将 $Q_{\mathrm{gs}}-Q_3$、$Q_{\mathrm{gs}}-Q_2$ 与生成分支 exBe10 比较，观察峰位、宽度和两分支的重叠。真值只有输入的两个离散分支；重建分布的宽度不等于该核能级的自然宽度。

第二、三行分别给出未补 γ 的两碎片不变质量与 Q2、Q3 的关联。横轴用于区分衰变分支，纵轴反映漏掉 γ 后的重建结果。先比较同一靶核在两种方法下的分支分离，再比较 H/C 的差异。图中 H/C 已按模拟真值分别绘制；是否能在混合样本中选出 H，需要下一节的关联分析来验证。

这些二维图同时提醒我们：Q 和不变质量使用了相同的粒子能量与方向，二者不是独立观测量，Q 门本身可能改变选中谱的形状。