$ {}^{14}\mathrm C^*\!\to{}^{4}\mathrm{He}+{}^{10}\mathrm{Be}^*{} $ 反应运动学模拟
Introduction
本节以 $^{14}\mathrm C$ 束流在含氢、碳的靶上发生非弹激发，随后衰变为 $\alpha+{}^{10}\mathrm{Be}^{(*)}$ 的反应为例。实验背景见 Han et al., Communications Physics 6, 220 (2023)。实验通过测量末态粒子，确定衰变分支并重建 $^{14}\mathrm C$ 的激发能。
这里保留两个反应过程： $$^{14}\mathrm C+{}^1\mathrm H\to{}^{14}\mathrm C^*+p,$$ $$^{14}\mathrm C+{}^{12}\mathrm C\to{}^{14}\mathrm C^*+{}^{12}\mathrm C,$$ 随后均有 $$^{14}\mathrm C^*\to\alpha+{}^{10}\mathrm{Be}^{(*)}.$$ 如果生成 $^{10}\mathrm{Be}$ 的 3.368 MeV 激发态，还生成其 $\gamma$ 退激：$^{10}\mathrm{Be}^*\to{}^{10}\mathrm{Be}+\gamma$。
模拟分两部分：先用 TGenPhaseSpace 生成满足四动量守恒的末态，再对带电粒子的能量、方向施加简化分辨。两类靶反应写入同一个 ROOT 文件，供 4.5 比较 Q 值重建方法、4.6 选择反应道并重建不变质量。processID 保存靶核真值，用于检查分类结果，不能当作实验中已经知道的靶核标签。
这是一组运动学与重建练习用的 toy MC。它不模拟探测器几何接受度、靶中能损、角分布的反应动力学或绝对产额；这些内容在 4.7 引入输运时再处理。
四动量表示与两步衰变
在自然单位 $c=1$ 下，单粒子四动量写作 $$ P=(E,\vec p),\qquad P^2=E^2-\vec p^{\,2}=m^2. $$
所有运动学关系都由四动量守恒给出。若初态总四动量为 $$ P_{\mathrm{tot}}=P_{\mathrm{beam}}+P_{\mathrm{target}}, $$ 则对过程 1，有 $$ P_{\mathrm{tot}}=P_{p}+P_{^{14}\mathrm C^*}, \qquad P_{^{14}\mathrm C^*}=P_{\alpha}+P_{^{10}\mathrm Be^*}. $$ 对过程 2，有 $$ P_{\mathrm{tot}}=P_{^{12}\mathrm C\mathrm{(recoil)}}+P_{^{14}\mathrm C^*}, \qquad P_{^{14}\mathrm C^*}=P_{\alpha}+P_{^{10}\mathrm Be^*}. $$ 若 $^{10}\mathrm Be^*$ 处于激发态，还要继续满足 $$ P_{^{10}\mathrm Be^*}=P_{^{10}\mathrm Be}+P_{\gamma}. $$ 程序实现时，先构造
compound = beam + target;
然后按过程类型连续调用 TGenPhaseSpace::SetDecay() 与 Generate() 来完成逐步衰变生成。
激发态质量的处理
若某核处于激发态，其参与运动学的质量应写成 $m^* = m_{\rm gs} + E_x,$ 其中 $m_{\rm gs}$ 为基态质量，$E_x$ 为激发能。本节中采用 $$ m_{{}^{14}\mathrm C^*}=m_{{}^{14}\mathrm C}+E_x({}^{14}\mathrm C),\qquad m_{{}^{10}\mathrm{Be}^*}=m_{{}^{10}\mathrm{Be}}+E_x({}^{10}\mathrm{Be}), $$ $^{10}\mathrm{Be}$ 激发能为 $0$ 和 $3.368$ MeV；$^{14}\mathrm C$ 则按不同分支抽取若干激发态，并用 Breit–Wigner 分布对各态的激发能做抽样。
衰变物理过程
下面列出的能级、宽度及分支权重构成教学输入表，不是对该文献能级表和实测分支比的逐项复现。br 是抽样的相对权重，width 是 Breit–Wigner 的 FWHM。物理阈值由每次 SetDecay() 检查；不允许的候选事件不进入输出树。
%jsroot on
flowchart LR
A["¹⁴C* → α + ¹⁰Be*"]
A --> B["¹⁰Be(gs)<br/>Ex = 0 MeV<br/>50%"]
A --> C["¹⁰Be*(3.368 MeV)<br/>50%"]
B --> B1["Ex(¹⁴C) = 14.9 MeV<br/>Γ = 0.10 MeV<br/>BR = 0.25"]
B --> B2["Ex(¹⁴C) = 15.6 MeV<br/>Γ = 0.18 MeV<br/>BR = 0.50"]
B --> B3["Ex(¹⁴C) = 16.4 MeV<br/>Γ = 0.14 MeV<br/>BR = 0.25"]
C --> C1["Ex(¹⁴C) = 17.3 MeV<br/>Γ = 0.12 MeV<br/>BR = 0.50"]
C --> C2["Ex(¹⁴C) = 18.5 MeV<br/>Γ = 0.08 MeV<br/>BR = 0.50"]
C --> G["¹⁰Be*(3.368 MeV) → ¹⁰Be(gs) + γ"]
G --> G1["<br/>Eγ ≈ 3.368 MeV（ ¹⁰Be* 静止系）"]
探测器能量分辨
用 Gaussian 表示单粒子动能读数的简化响应： $$\sigma_E^2=\omega F T+\sigma_{\mathrm{noise}}^2.$$ 取 $F=0.12$、$\omega=3.6\times10^{-3}$ keV、$\sigma_{\mathrm{noise}}=50$ keV。式中 $T$ 也用 keV；代码内部动能用 GeV，因此先乘 $10^6$，计算和抽样后再换回 GeV。这里的 $T$ 被视作一个能量读数，不逐层模拟能量沉积；这些数值是本练习的响应假设，不代表完整望远镜的实测分辨。
eRes 保留正值重抽样的简化处理，因此其输出是正值截断 Gaussian。远高于噪声时影响很小，接近阈值时则不能把它当作无偏的读出模型。
角分辨模型
分别在极角与方位角上加均匀扰动：
$$\delta\theta,\delta\phi\sim U(-\Delta/2,\Delta/2),\qquad\Delta=1^\circ.$$
这是坐标角上的教学展宽，标准差为 $\Delta/\sqrt{12}$，不是 $1^\circ$ 的 Gaussian 分辨，也不是方向球面上的各向同性展宽。angleRes 返回单位方向；随后由测量动能和质量重新计算动量大小。真实位置探测器的角分辨还取决于反应点、探测器距离和位置分辨。
输入参数
束流平均动能为 317.1 MeV（整颗 $^{14}\mathrm C$，不是 MeV/u），本练习取束流 FWHM 为 1.83 MeV，换算为 $\sigma=1.83/2.355$ MeV。模拟分别尝试生成 $10^6$ 个 H 靶和 C 靶事件；运动学阈值和下述选择可能使实际保存数减少。
将 mass.txt 与 macro 放在同一目录。每行依次为 Z A mass，第三列是核的静止能量，单位 keV，不是质量过剩。例如 $(Z,A)=(1,1)$ 的值约为 938272 keV。getMass 返回 GeV；写入树的动能、激发能使用 MeV，角度使用度。
程序流程
每个事件沿同一条流程处理：
- 抽取束流动能，构造束流与静止靶核的四动量。
- 按输入表选择 $^{10}\mathrm{Be}$ 分支和 $^{14}\mathrm C^*$ 激发能。
- 依次生成反冲核、$^{14}\mathrm C^*$ 及其衰变碎片；激发态 $^{10}\mathrm{Be}$ 再做 γ 退激。
- 检查三个带电末态的真值动能均不低于 0.1 MeV，随后施加读出展宽。这里保留这一生成样本选择，不把它称为真实电子学触发阈值。
- 写入同一个
tree，保存最终带电粒子的动能、方向及用于检验的真值。
每次两体 Generate() 的默认归一化权重约为 1；保留 weight_total 分支是为了保持与前面加权分析的接口一致。H/C 的相对事件数和能级分支概率来自教学输入，不从这些权重推导实验产额。
当前保存样本已要求反冲核通过真值动能下限，因此 4.6 的“两体分析”指不使用反冲核测量量，并不等于未经反冲核选择的完整两体触发样本。
代码实现
先看数据在程序中的含义，再读完整实现：
| 变量或函数 | 作用 |
|---|---|
C14Level、Be10Branch |
保存输入的激发能、宽度和相对抽样权重 |
getMass(Z,A) |
从 mass.txt 读取核质量，并把 keV 换成 GeV |
SelectIndex、SampleBreitWigner |
分别抽取离散分支和连续激发能 |
FillChargedParticle |
将四动量转为动能和角度，并施加读出展宽 |
ResetEvent |
每个事件开始时清空缓存，避免留下上一事件的值 |
Simulation |
串联上述过程，逐事件 Fill()，最后保存 ROOT 文件 |
ek[0]、ek[1]、ek[2] 分别是 α、最终 $^{10}\mathrm{Be}$ 和反冲核的模拟读数；theta、phi 使用相同下标。ek[3] 及对应角度只在 hasGamma==1 时保存 γ 退激前的 $^{10}\mathrm{Be}^*$ 真值。exC14、exBe10、ekBeam、eGamma 也是真值或生成输入，不是全部可供实验重建使用的测量量。
完整实现依照上述顺序分段注释。运行页面顶部的 c14_simulation.C 会在当前目录生成 C14_CHn_He4Be10x.root；同名旧输出会被覆盖，需保留时请先另存。
%%cpp -d
#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include "TFile.h"
#include "TTree.h"
#include "TParameter.h"
#include "TRandom.h"
#include "TMath.h"
#include "TLorentzVector.h"
#include "TVector3.h"
#include "TGenPhaseSpace.h"
using namespace std;
// 数据结构
struct C14Level {
double ex; // MeV
double width; // MeV
double br; // 分支比
};
struct Be10Branch {
double ex; // MeV
double br; // 分支比
vector<C14Level> levels; // 该 10Be 支路下对应的 14C* 能级
};
// 质量读取
%%cpp -d
double getMass(int zz, int aa, const char* inputFile = "mass.txt")
{
const double keV2GeV = 1.0e-6;
ifstream inFile(inputFile);
if (!inFile.is_open()) {
cerr << "Cannot open mass file: " << inputFile << endl;
return -1.0;
}
int Z, A;
double m;
while (inFile >> Z >> A >> m) {
if (Z == zz && A == aa) return m * keV2GeV;
}
cerr << "Mass not found for Z=" << zz << ", A=" << aa << endl;
return -1.0;
}
%%cpp -d
int SelectIndex(const vector<double>& ratios)
{
double sum = 0.0;
for (size_t i = 0; i < ratios.size(); ++i) sum += ratios[i];
if (sum <= 0.0) return -1;
double r = gRandom->Uniform(0.0, sum);
double cumulative = 0.0;
for (size_t i = 0; i < ratios.size(); ++i) {
cumulative += ratios[i];
if (r < cumulative) return (int)i;
}
return (int)ratios.size() - 1;
}
%%cpp -d
double SampleBreitWigner(double mean, double width,
double xmin, double xmax,
int maxRetry = 100)
{
for (int i = 0; i < maxRetry; ++i) {
double x = gRandom->BreitWigner(mean, width);
if (x >= xmin && x < xmax) return x;
}
return -1.0;
}
%%cpp -d
double eRes(double eraw, bool addRes = true) // input: GeV
{
if (!addRes) return eraw;
const double keVtoGeV = 1.0e-6;
const double GeVtokeV = 1.0e6;
const double sig_noise = 50.0; // keV
const double fano_factor = 0.12;
const double omega = 3.6e-3; // keV
const double energy_keV = eraw * GeVtokeV;
const double sig_int = TMath::Sqrt(omega * fano_factor * energy_keV);
const double sig_tot = TMath::Sqrt(sig_noise * sig_noise + sig_int * sig_int);
double smeared_keV = -1.0;
do {
smeared_keV = gRandom->Gaus(energy_keV, sig_tot);
} while (smeared_keV <= 0.0);
return smeared_keV * keVtoGeV;
}
%%cpp -d
TVector3 angleRes(const TVector3& pvraw, bool addRes = true)
{
TVector3 dir = pvraw.Unit();
if (!addRes) return dir;
const double deg = TMath::Pi() / 180.0;
const double delta = 1.0 * deg;
double theta = dir.Theta();
double phi = dir.Phi();
theta += gRandom->Uniform(-0.5 * delta, 0.5 * delta);
phi += gRandom->Uniform(-0.5 * delta, 0.5 * delta);
if (theta < 0.0) theta = -theta;
if (theta > TMath::Pi()) theta = 2.0 * TMath::Pi() - theta;
TVector3 out;
out.SetMagThetaPhi(1.0, theta, phi);
return out.Unit();
}
%%cpp -d
void FillChargedParticle(const TLorentzVector& lv, double massGeV,
double& ekMeV, double& thetaDeg, double& phiDeg,
bool smearE = true, bool smearA = true)
{
const double MeV = 1.0e3;
const double RadToDeg = 180.0 / TMath::Pi();
double Traw = lv.E() - massGeV;
TVector3 dir = angleRes(lv.Vect(), smearA);
ekMeV = eRes(Traw, smearE) * MeV;
thetaDeg = dir.Theta() * RadToDeg;
phiDeg = dir.Phi() * RadToDeg;
if (phiDeg < 0.0) phiDeg += 360.0;
}
%%cpp -d
void ResetEvent(double ek[4], double theta[4], double phi[4],
double weight[3], int& hasGamma,
double& eGamma, double& weight_total,
double& exBe10, double& exC14)
{
for (int i = 0; i < 4; ++i) {
ek[i] = -999.0;
theta[i] = -999.0;
phi[i] = -999.0;
}
for (int i = 0; i < 3; ++i) weight[i] = 1.0;
hasGamma = 0;
eGamma = 0.0;
weight_total = 1.0;
exBe10 = -999.0;
exC14 = -999.0;
}
%%cpp -d
vector<Be10Branch> BuildDecayTable()
{
vector<Be10Branch> table;
Be10Branch b0;
b0.ex = 0.0;
b0.br = 0.50;
b0.levels.push_back({14.9, 0.10, 0.25});
b0.levels.push_back({15.6, 0.18, 0.50});
b0.levels.push_back({16.4, 0.14, 0.25});
table.push_back(b0);
Be10Branch b1;
b1.ex = 3.368;
b1.br = 0.50;
b1.levels.push_back({17.3, 0.12, 0.50});
b1.levels.push_back({18.5, 0.08, 0.50});
table.push_back(b1);
return table;
}
先设置质量、束流和衰变分支。固定分支表对应的 be10Ratios 在事件循环外构造一次；抽样次序和随机数种子保持固定。
gRandom->SetSeed(4401); // 固定种子，便于复现教学结果
const double GeV = 1.0e-3; // MeV -> GeV
const double MeV = 1.0e3; // GeV -> MeV
const int PROCESS_H = 1;
const int PROCESS_C = 2;
const double exC14Min = 0.0;
const double exC14Max = 30.0;
const double thresholdGeV = 100.0e-6; // 100 keV
// Ground-state masses (GeV)
const double mH1 = getMass(1, 1);
const double mHe4 = getMass(2, 4);
const double mBe10 = getMass(4, 10);
const double mC12 = getMass(6, 12);
const double mC14 = getMass(6, 14);
const double mGamma = 0.0;
if (mH1 < 0 || mHe4 < 0 || mBe10 < 0 || mC12 < 0 || mC14 < 0) {
cout << "Mass loading failed." << endl;
return;
}
// 衰变分支表
vector<Be10Branch> decayTable = BuildDecayTable();
vector<double> be10Ratios;
for (const Be10Branch &branch : decayTable) be10Ratios.push_back(branch.br);
// Beam input
const double ekBeamMean = 317.1; // MeV
const double ekBeamSigma = 1.83 / 2.355; // MeV
// Statistics
const Long64_t nEventEach = 1000000;
// Output
建立输出 Tree。数组的四个位置依次保存 α、末态 10Be、反冲核，以及退激前的 10Be*。
TFile* fout = new TFile("C14_CHn_He4Be10x.root", "recreate");
TTree* simTree = new TTree("tree", "toy MC for 14C on (CH)n target");
Int_t processID; // 1: H 靶, 2: C 靶
Int_t hasGamma; // 0: 无 gamma, 1: 有 gamma
Double_t ekBeam; // MeV
Double_t ek[4]; // MeV
Double_t theta[4]; // degree
Double_t phi[4]; // degree
Double_t exBe10; // MeV
Double_t exC14; // MeV
Double_t eGamma; // MeV, lab frame
Double_t weight[3];
Double_t weight_total;
simTree->Branch("processID", &processID, "processID/I");
simTree->Branch("hasGamma", &hasGamma, "hasGamma/I");
simTree->Branch("ekBeam", &ekBeam, "ekBeam/D");
simTree->Branch("ek", ek, "ek[4]/D");
simTree->Branch("theta", theta, "theta[4]/D");
simTree->Branch("phi", phi, "phi[4]/D");
simTree->Branch("exBe10", &exBe10, "exBe10/D");
simTree->Branch("exC14", &exC14, "exC14/D");
simTree->Branch("eGamma", &eGamma, "eGamma/D");
simTree->Branch("weight", weight, "weight[3]/D");
simTree->Branch("weight_total", &weight_total, "weight_total/D");
fout->cd();
TParameter<Double_t>("massH1", mH1).Write();
TParameter<Double_t>("massC12", mC12).Write();
TParameter<Double_t>("massHe4", mHe4).Write();
TParameter<Double_t>("massBe10", mBe10).Write();
TParameter<Double_t>("massC14", mC14).Write();
TParameter<Double_t>("ekBeamMean", ekBeamMean).Write();
TParameter<Long64_t>("nEventEach", nEventEach).Write();
// 两个过程写入同一个 TTree
//
// 统一约定：
// ek[0],theta[0],phi[0] : alpha
// ek[1],theta[1],phi[1] : final 10Be
// ek[2],theta[2],phi[2] : recoil (H 靶为 p, C 靶为 12C)
// ek[3],theta[3],phi[3] : gamma 退激前的 10Be*（仅 hasGamma=1 时有意义）
逐事件生成反应和衰变，再加入已定义的能量、角度响应。每段注释对应一个物理步骤，通过条件的事件写入 simTree。
for (int proc = PROCESS_H; proc <= PROCESS_C; ++proc) {
Long64_t entriesBefore = simTree->GetEntries();
const double mTarget = (proc == PROCESS_H) ? mH1 : mC12;
const double mRecoil = (proc == PROCESS_H) ? mH1 : mC12;
TLorentzVector target(0.0, 0.0, 0.0, mTarget);
for (Long64_t n = 0; n < nEventEach; ++n) {
processID = proc;
ResetEvent(ek, theta, phi, weight, hasGamma, eGamma, weight_total, exBe10, exC14);
// 1. Beam
double Tbeam_GeV = gRandom->Gaus(ekBeamMean * GeV, ekBeamSigma * GeV);
if (Tbeam_GeV <= 0.0) continue;
double Ebeam_GeV = Tbeam_GeV + mC14;
double pbeam_GeV = TMath::Sqrt(Ebeam_GeV * Ebeam_GeV - mC14 * mC14);
TLorentzVector beam(0.0, 0.0, pbeam_GeV, Ebeam_GeV);
TLorentzVector compound = beam + target;
// 2. 选 10Be 支路
int iBe10 = SelectIndex(be10Ratios);
if (iBe10 < 0) continue;
exBe10 = decayTable[iBe10].ex;
// 3. 在选中的 10Be 支路下，再选 14C* 能级
vector<double> c14Ratios;
for (size_t i = 0; i < decayTable[iBe10].levels.size(); ++i) {
c14Ratios.push_back(decayTable[iBe10].levels[i].br);
}
int iC14 = SelectIndex(c14Ratios);
if (iC14 < 0) continue;
C14Level level = decayTable[iBe10].levels[iC14];
exC14 = SampleBreitWigner(level.ex, level.width, exC14Min, exC14Max, 100);
if (exC14 < 0.0) continue;
// 4. compound -> recoil + 14C*
TGenPhaseSpace decay1;
double masses1[2] = {mRecoil, mC14 + exC14 * GeV};
if (!decay1.SetDecay(compound, 2, masses1)) continue;
weight[0] = decay1.Generate();
TLorentzVector* recoil = decay1.GetDecay(0);
TLorentzVector* c14x = decay1.GetDecay(1);
// 5. 14C* -> alpha + 10Be*
TGenPhaseSpace decay2;
double masses2[2] = {mHe4, mBe10 + exBe10 * GeV};
if (!decay2.SetDecay(*c14x, 2, masses2)) continue;
weight[1] = decay2.Generate();
TLorentzVector* alpha = decay2.GetDecay(0);
TLorentzVector* be10x = decay2.GetDecay(1);
// 6. 若 10Be 为激发态，则总是继续做 gamma 退激
TLorentzVector be10_final = *be10x;
if (exBe10 > 1.0e-6) {
TGenPhaseSpace decay3;
double masses3[2] = {mBe10, mGamma};
if (!decay3.SetDecay(*be10x, 2, masses3)) continue;
weight[2] = decay3.Generate();
TLorentzVector* be10 = decay3.GetDecay(0);
TLorentzVector* gamma = decay3.GetDecay(1);
be10_final = *be10;
eGamma = gamma->E() * MeV; // 实验室系 gamma 能量
hasGamma = 1;
}
weight_total = weight[0] * weight[1] * weight[2];
// 7. 阈值判断
double Traw_alpha = alpha->E() - mHe4;
double Traw_recoil = recoil->E() - mRecoil;
double Traw_be10 = 0.0;
if (hasGamma) {
Traw_be10 = be10_final.E() - mBe10;
} else {
Traw_be10 = be10x->E() - (mBe10 + exBe10 * GeV);
}
if (Traw_alpha < thresholdGeV) continue;
if (Traw_be10 < thresholdGeV) continue;
if (Traw_recoil < thresholdGeV) continue;
// 8. 统一写树
FillChargedParticle(*alpha, mHe4, ek[0], theta[0], phi[0]);
if (hasGamma) {
FillChargedParticle(be10_final, mBe10, ek[1], theta[1], phi[1]);
FillChargedParticle(*be10x, mBe10 + exBe10 * GeV,
ek[3], theta[3], phi[3],
false, false);
} else {
FillChargedParticle(*be10x, mBe10 + exBe10 * GeV,
ek[1], theta[1], phi[1]);
}
FillChargedParticle(*recoil, mRecoil, ek[2], theta[2], phi[2]);
ekBeam = Tbeam_GeV * MeV;
simTree->Fill();
}
cout << "processID=" << proc << ": attempted=" << nEventEach
<< ", saved=" << simTree->GetEntries()-entriesBefore << endl;
}
simTree->Write();
fout->Close();
cout << "Simulation finished: C14_CHn_He4Be10x.root" << endl;
processID=1: attempted=1000000, saved=991356 processID=2: attempted=1000000, saved=994314 Simulation finished: C14_CHn_He4Be10x.root
结果画图
TFile *f = new TFile("C14_CHn_He4Be10x.root");
if (f->IsZombie()) throw std::runtime_error("Run the generation cells above first");
TTree *tree = f->Get<TTree>("tree");
tree->Print(); // 一行对应一个通过生成条件的事件
gStyle->SetOptStat(0);
gStyle->SetPalette(kBird);
gStyle->SetPadLeftMargin(0.15);
gStyle->SetPadRightMargin(0.16);
****************************************************************************** *Tree :tree : toy MC for 14C on (CH)n target * *Entries : 1985670 : Total = 333749050 bytes File Size = 207346880 * * : : Tree compression factor = 1.61 * ****************************************************************************** *Br 0 :processID : processID/I * *Entries : 1985670 : Total Size= 7947705 bytes File Size = 43993 * *Baskets : 49 : Basket Size= 838656 bytes Compression= 180.63 * *............................................................................* *Br 1 :hasGamma : hasGamma/I * *Entries : 1985670 : Total Size= 7947652 bytes File Size = 850723 * *Baskets : 49 : Basket Size= 838656 bytes Compression= 9.34 * *............................................................................* *Br 2 :ekBeam : ekBeam/D * *Entries : 1985670 : Total Size= 15893582 bytes File Size = 13613910 * *Baskets : 85 : Basket Size= 1676288 bytes Compression= 1.17 * *............................................................................* *Br 3 :ek : ek[4]/D * *Entries : 1985670 : Total Size= 63568625 bytes File Size = 54490212 * *Baskets : 302 : Basket Size= 6704128 bytes Compression= 1.17 * *............................................................................* *Br 4 :theta : theta[4]/D * *Entries : 1985670 : Total Size= 63569543 bytes File Size = 54610910 * *Baskets : 302 : Basket Size= 6705152 bytes Compression= 1.16 * *............................................................................* *Br 5 :phi : phi[4]/D * *Entries : 1985670 : Total Size= 63568931 bytes File Size = 54554490 * *Baskets : 302 : Basket Size= 6704640 bytes Compression= 1.17 * *............................................................................* *Br 6 :exBe10 : exBe10/D * *Entries : 1985670 : Total Size= 15893582 bytes File Size = 1065457 * *Baskets : 85 : Basket Size= 1676288 bytes Compression= 14.92 * *............................................................................* *Br 7 :exC14 : exC14/D * *Entries : 1985670 : Total Size= 15893493 bytes File Size = 14530834 * *Baskets : 85 : Basket Size= 1676288 bytes Compression= 1.09 * *............................................................................* *Br 8 :eGamma : eGamma/D * *Entries : 1985670 : Total Size= 15893582 bytes File Size = 8667480 * *Baskets : 85 : Basket Size= 1676288 bytes Compression= 1.83 * *............................................................................* *Br 9 :weight : weight[3]/D * *Entries : 1985670 : Total Size= 47677793 bytes File Size = 3097251 * *Baskets : 230 : Basket Size= 5028864 bytes Compression= 15.39 * *............................................................................* *Br 10 :weight_total : weight_total/D * *Entries : 1985670 : Total Size= 15894116 bytes File Size = 1806347 * *Baskets : 85 : Basket Size= 1676800 bytes Compression= 8.80 * *............................................................................*
Ek vs theta
分别比较 H/C 两类反应的动能—极角关联，判断前向碎片与反冲核的角度覆盖需求。这些图还未施加探测器几何接受度。
// c1: H/C target kinematics
// 2 rows x 3 cols, each pad ~300x300
TCanvas *c1 = new TCanvas("c1","H/C target kinematics",900,600);
c1->Divide(3,2);
// ---- H target ----
c1->cd(1);
gPad->SetRightMargin(0.12);
tree->Draw("ek[0]:theta[0]>>h1_H_a(300,0,120,300,0,140)",
"(processID==1)*weight_total","colz");
((TH2*)gROOT->FindObject("h1_H_a"))->SetTitle("H: #alpha;#theta_{#alpha} [deg];E_{k,#alpha} [MeV]");
c1->cd(2);
gPad->SetRightMargin(0.12);
tree->Draw("ek[1]:theta[1]>>h1_H_b(300,0,120,300,0,260)",
"(processID==1)*weight_total","colz");
((TH2*)gROOT->FindObject("h1_H_b"))->SetTitle("H: ^{10}Be;#theta_{^{10}Be} [deg];E_{k,^{10}Be} [MeV]");
c1->cd(3);
gPad->SetRightMargin(0.12);
tree->Draw("ek[2]:theta[2]>>h1_H_r(300,0,120,300,0,80)",
"(processID==1)*weight_total","colz");
((TH2*)gROOT->FindObject("h1_H_r"))->SetTitle("H: recoil p;#theta_{recoil} [deg];E_{k,recoil} [MeV]");
// ---- C target ----
c1->cd(4);
gPad->SetRightMargin(0.12);
tree->Draw("ek[0]:theta[0]>>h1_C_a(300,0,120,300,0,140)",
"(processID==2)*weight_total","colz");
((TH2*)gROOT->FindObject("h1_C_a"))->SetTitle("C: #alpha;#theta_{#alpha} [deg];E_{k,#alpha} [MeV]");
c1->cd(5);
gPad->SetRightMargin(0.12);
tree->Draw("ek[1]:theta[1]>>h1_C_b(300,0,120,300,0,260)",
"(processID==2)*weight_total","colz");
((TH2*)gROOT->FindObject("h1_C_b"))->SetTitle("C: ^{10}Be;#theta_{^{10}Be} [deg];E_{k,^{10}Be} [MeV]");
c1->cd(6);
gPad->SetRightMargin(0.12);
tree->Draw("ek[2]:theta[2]>>h1_C_r(300,0,120,300,0,320)",
"(processID==2)*weight_total","colz");
((TH2*)gROOT->FindObject("h1_C_r"))->SetTitle("C: recoil ^{12}C;#theta_{recoil} [deg];E_{k,recoil} [MeV]");
c1->Draw();
Energy-energy correlations
同一事件中的粒子能量不能独立任意变化；关联带体现总能量和动量守恒。绘图使用 weight_total，标签标明当前靶核与粒子组合。
TCanvas *c2 = new TCanvas("c2","Energy-energy correlations",600,600);
c2->Divide(2,2);
// ---- H target ----
c2->cd(1);
gPad->SetRightMargin(0.12);
tree->Draw("ek[1]:ek[0]>>h2_H_ab(300,0,140,300,0,260)",
"(processID==1)*weight_total","colz");
((TH2*)gROOT->FindObject("h2_H_ab"))->SetTitle("H: E_{k}(^{10}Be) vs E_{k}(#alpha);E_{k,#alpha} [MeV];E_{k,^{10}Be} [MeV]");
c2->cd(2);
gPad->SetRightMargin(0.12);
tree->Draw("ek[2]:ek[0]>>h2_H_ar(300,0,140,300,0,80)",
"(processID==1)*weight_total","colz");
((TH2*)gROOT->FindObject("h2_H_ar"))->SetTitle("H: E_{k}(recoil) vs E_{k}(#alpha);E_{k,#alpha} [MeV];E_{k,recoil} [MeV]");
// ---- C target ----
c2->cd(3);
gPad->SetRightMargin(0.12);
tree->Draw("ek[1]:ek[0]>>h2_C_ab(300,0,140,300,0,260)",
"(processID==2)*weight_total","colz");
((TH2*)gROOT->FindObject("h2_C_ab"))->SetTitle("C: E_{k}(^{10}Be) vs E_{k}(#alpha);E_{k,#alpha} [MeV];E_{k,^{10}Be} [MeV]");
c2->cd(4);
gPad->SetRightMargin(0.12);
tree->Draw("ek[2]:ek[0]>>h2_C_ar(300,0,140,300,0,320)",
"(processID==2)*weight_total","colz");
((TH2*)gROOT->FindObject("h2_C_ar"))->SetTitle("C: E_{k}(recoil) vs E_{k}(#alpha);E_{k,#alpha} [MeV];E_{k,recoil} [MeV]");
c2->Draw();
Angle-angle correlations
图中分别给出两个碎片、反冲核与 α 的极角关联，不是直接绘制两粒子的空间夹角。先看坐标轴的粒子下标，再对照 ek/theta/phi 的统一定义；计算空间夹角时还需使用方位角。
TCanvas *c3 = new TCanvas("c3","Angle-angle correlations",600,600);
c3->Divide(2,2);
// ---- H target ----
c3->cd(1);
gPad->SetRightMargin(0.12);
tree->Draw("theta[1]:theta[0]>>h3_H_ab(300,0,120,300,0,120)",
"(processID==1)*weight_total","colz");
((TH2*)gROOT->FindObject("h3_H_ab"))->SetTitle("H: #theta(^{10}Be) vs #theta(#alpha);#theta_{#alpha} [deg];#theta_{^{10}Be} [deg]");
c3->cd(2);
gPad->SetRightMargin(0.12);
tree->Draw("theta[2]:theta[0]>>h3_H_ar(300,0,120,300,0,120)",
"(processID==1)*weight_total","colz");
((TH2*)gROOT->FindObject("h3_H_ar"))->SetTitle("H: #theta(recoil) vs #theta(#alpha);#theta_{#alpha} [deg];#theta_{recoil} [deg]");
// ---- C target ----
c3->cd(3);
gPad->SetRightMargin(0.12);
tree->Draw("theta[1]:theta[0]>>h3_C_ab(300,0,120,300,0,120)",
"(processID==2)*weight_total","colz");
((TH2*)gROOT->FindObject("h3_C_ab"))->SetTitle("C: #theta(^{10}Be) vs #theta(#alpha);#theta_{#alpha} [deg];#theta_{^{10}Be} [deg]");
c3->cd(4);
gPad->SetRightMargin(0.12);
tree->Draw("theta[2]:theta[0]>>h3_C_ar(300,0,120,300,0,120)",
"(processID==2)*weight_total","colz");
((TH2*)gROOT->FindObject("h3_C_ar"))->SetTitle("C: #theta(recoil) vs #theta(#alpha);#theta_{#alpha} [deg];#theta_{recoil} [deg]");
c3->Draw();
Input structure check
这组图检查生成输入：激发能分支、目标事件数及其关联。它们来自模拟真值，作用是确认生成器按设定工作，不是用实验数据测得的能级强度。
TCanvas *c4 = new TCanvas("c4","Input structure check",900,900);
c4->Divide(2,2);
c4->cd(1);
tree->Draw("exC14>>h8_exC14_H(400,13,20)",
"(processID==1)*weight_total","hist");
((TH1*)gROOT->FindObject("h8_exC14_H"))->SetTitle("H target: E_{x}(^{14}C);E_{x}(^{14}C) [MeV];weighted counts");
c4->cd(2);
tree->Draw("exC14>>h8_exC14_C(400,13,20)",
"(processID==2)*weight_total","hist");
((TH1*)gROOT->FindObject("h8_exC14_C"))->SetTitle("C target: E_{x}(^{14}C);E_{x}(^{14}C) [MeV];weighted counts");
c4->cd(3);
tree->Draw("exBe10>>h8_exBe10(120,0,5)",
"weight_total","hist");
((TH1*)gROOT->FindObject("h8_exBe10"))->SetTitle("^{10}Be excitation;E_{x}(^{10}Be) [MeV];weighted counts");
c4->cd(4);
tree->Draw("eGamma>>h8_g(300,0,10)",
"(hasGamma==1)*weight_total","hist");
((TH1*)gROOT->FindObject("h8_g"))->SetTitle("#gamma energy in lab;E_{#gamma}^{lab} [MeV];weighted counts");
c4->Draw();
Appendix：TTree::Draw 中的 selection
在
tree->Draw("varexp", "selection")
中，selection 不只是“筛选条件”，更准确地说，它是每个事件的填图权重表达式。
selection 同时可以实现：
- 选事件
- 加权填图
常见写法：
tree->Draw("x>>h", "cut");
表示只做筛选；满足条件的事件按权重 1 填图。
tree->Draw("x>>h", "(cut)*weight");
表示先筛选，再按 weight 加权填图。
一句话总结：
selection 应理解为逐事件的权重表达式：0 表示不填，1 表示普通填图，其他数值表示带权填图。