反应 $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} $$
因此，当某个激发态的退激辐射没有被实验直接探测到时，其激发能将表现为重建 $Q$ 值相对于 $Q_{gs}$ 的偏移；反过来，也可以通过 $Q$ 值重建来反推出该激发能。
$^{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. $$
若该退激 $\gamma$ 未被探测，则缺失的激发能不会直接出现在带电粒子的动能求和中，而会表现为 $Q$ 值相对于 $Q_{gs}$ 的偏移。因此，本节 $Q$ 值分析的基本思想就是：利用带电末态粒子的运动学信息重建 $Q$ 值，并由其偏移反推出 $^{10}\mathrm{Be}^{*}$ 的激发能。
$Q$ 值测量对于不变质量谱（Invariant Mass）的意义¶
在实验中，母核 $^{14}\mathrm{C}^*$ 的激发能通常通过 $\alpha$ 与 $^{10}\mathrm{Be}$ 的不变质量谱（Invariant Mass Spectrum） 来重建： $$ M_{inv} = \sqrt{(\sum E_i)^2 - (\sum \vec{p}_i)^2} $$ $$ E_x(^{14}\mathrm{C}) = M_{inv} - m_{14\mathrm{C}, gs} $$ 但是，这要求准确知道末态粒子的状态。如果 $^{10}\mathrm{Be}$ 处于激发态（例如 3.368 MeV）并释放了未被探测的 $\gamma$ 光子，而我们在计算不变质量时依然错误地假设它处于基态且未丢失能量，那么重建出的 $^{14}\mathrm{C}^*$ 激发能将会产生明显偏移，导致物理假峰或能级错位。
因此，$Q$ 值重建是不变质量谱分析的严格前置步骤：
- 反应道甄别：通过重建的 $Q$ 值（如 $Q \approx Q_{gs}$ 或是 $Q \approx Q_{gs} - 3.368\text{ MeV}$），我们可以逐事件（Event-by-event）判定该事件中的 $^{10}\mathrm{Be}$ 到底处于基态还是激发态。
- 修正激发能：一旦通过 $Q$ 值谱确认了 $^{10}\mathrm{Be}^*$ 的分支，就能在后续的不变质量计算中，对那些伴随 $\gamma$ 丢失的事件补偿相应的激发能（或剔除对应本底），从而提取出纯净且位置准确的母核 $^{14}\mathrm{C}^*$ 激发态能级。
$Q$ 值重建方法¶
反应初态¶
束流动能 $E_{K,^{14}\mathrm C}$ 在实验中通常不能逐事件直接测量，其平均值 $\bar E_{K,^{14}\mathrm C}$ 由加速器的偶极磁铁 $B\rho$ 给出。对于稳定束流，能散可以近似忽略，$\Delta E_{K,^{14}\mathrm C}\to 0;$
对于放射性束流，则通常有$\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}}|$，从而反推出束流动能
$$ 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}$ 激发能¶
对于本节关注的 sequential decay 过程，若将母核衰变到 $\alpha + {}^{10}\mathrm{Be}^{*}$ 看作中间步，则在理想情况下有
$$ Q = Q_{\mathrm{gs}} - E_x(^{10}\mathrm{Be}), $$
因此可由重建的 Q 值反推出
$$ E_x(^{10}\mathrm{Be}) = Q_{\mathrm{gs}} - Q . $$
下面代码中exBe10 给出真实输入的 $^{10}\mathrm{Be}$ 激发能，因此可以将重建得到的
$$ E_x^{\mathrm{reco}}(^{10}\mathrm{Be}) = Q_{\mathrm{gs}} - Q $$
与 exBe10 逐一比较，以检验不同 Q 重建方法的分辨能力和适用范围。
%%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"
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 P(double ek, double m) const { return std::sqrt(ek*ek + 2*m*ek); }
double m(int pid, int i, double ex) const {
double mass[] = {mHe4, mBe10, mRecoil(pid), mBe10 + ex};
return mass[i];
}
TVector3 pv(int pid, int i, double ek[], double th[], double ph[], double ex) const {
TVector3 v;
v.SetMagThetaPhi(P(ek[i], m(pid,i,ex)), th[i]*TMath::DegToRad(), ph[i]*TMath::DegToRad());
return v;
}
TLorentzVector lv(int pid, int i, double ek[], double th[], double ph[], double ex) const {
return TLorentzVector(pv(pid,i,ek,th,ph,ex), ek[i] + m(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 = pv(pid,0,ek,th,ph,ex) + pv(pid,1,ek,th,ph,ex) + pv(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,P(ekBeamMean,mC14)) - pv(pid,0,ek,th,ph,ex) - pv(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 (lv(pid,0,ek,th,ph,ex) + lv(pid,1,ek,th,ph,ex)).M() - mC14;
}
};
void DrawQAnalysis()
{
TFile *f = TFile::Open("C14_CHn_He4Be10x.root");
TTree *t = (TTree*)f->Get("tree");
Kine K;
const double u = 1000.0; // GeV -> MeV，若原文件已是 MeV 改成 1.0
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]=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);
// 填充
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);
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, K.ExInv(pid, ek, th, ph, exBe10), w);
hInvQ2[j]->Fill(q2, K.ExInv(pid, ek, th, ph, exBe10), w);
}
// 画图
gStyle->SetOptStat(0);
int colH = kBlue+1, colC = kRed+1;
auto 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();
};
// ================ c1: Q 值谱对比 (2x2, logy) ================
TCanvas *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) ================
TCanvas *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);
}
}
DrawQAnalysis();
c1->Draw();
不同 Q 值重建方法的对比¶
展示如何克服束流能散带来的分辨率恶化，以及如何利用运动学差异区分靶本底。
- [左上] Qevt（理想 Q 值） 使用模拟的真实逐事件束流能量计算。H 和 C 靶均呈现尖锐的峰，代表探测系统在剔除束流能散影响后的极限物理分辨率。
- [右上] Qmean（平均束流 Q 值） 由于实验中用固定的平均能量代替真实的入射束流能量，受放射性束流固有的能散影响，H 和 C 靶的 Q 值分辨均出现严重的退化。
- [左下] Q3（三体动量守恒重建 Q 值）
利用出射三粒子动量和反推束流能量，逐事件修正能散。
- H 靶 (蓝线)：动量守恒有效抵消了能散的影响，基态与激发态峰位清晰锐利。
- C 靶 (红线)：出现显著的运动学展宽。由于反冲核 $^{12}\text{C}$ 质量较大且在厚靶中的能损和歧离较严重，探测系统微小的能量/角度误差被放大为显著的动量矢量偏差，导致分辨率大幅退化。
- [右下] Q2（二体缺失质量重建 Q 值） 模拟反冲核未被探测到的妥协方案，需引入平均束流动量推算缺失参数。由于未完全消除束流能散，分辨率介于 Qmean 和 Q3 之间。
c2->Draw();
重建激发能与不变质量¶
展示 Q 值测量对末态碎片和母核激发能重建的作用。
[第一行：左 vs 右] $E_x(^{10}\text{Be})$ 激发能一维谱 这两张图展示了通过测量的动量反推末态碎片 $^{10}\text{Be}$ 激发能的效果。
- 左 (H 靶)：红线（利用三体动量 Q3 重建）尖锐，说明对氢靶的重建精度高。
- 右 (C 靶)：因为碳靶的反冲核很重，存在严重的运动学展宽，导致这里重建出来的激发能变成了一个宽的分布。
[第二行：左 vs 右] Invariant Mass vs Q2 二维关联图 这两张图展示了提取母核 $^{14}\text{C}^*$ 激发能（Y轴，不变质量）时，为什么必须依赖 Q 值（X轴）。
- 左 (H 靶信号)：物理信号清晰地聚集在特定的垂直带上（即真实的 Q 值处）。当事例伴随 $\gamma$ 射线逃逸时，数据会落在左侧紫线处，此时 Y 轴计算出的母核激发能会整体偏低。但激发态与基态有足够的区分度。
- 右 (C 靶本底)：由于 C 靶的 Q 值发生严重展宽，其事例在图中弥散成很宽的区域。基态与激发态重叠。对比左图可知，只要在 X 轴的垂直物理线附近切一个狭窄区域，能够区分基态和激发态的贡献。
[第三行：左 vs 右] Invariant Mass vs Q3 二维关联图