基于 Geant4 的实验模拟¶
1 本节目标¶
本节对文献[Commun Phys 6, 220 (2023)]的实验进行 Geant4 模拟，以服务实验设计和实验分析。前面的 4.4 已经建立了反应学模型，4.5 与 4.6 分别给出了 Q 值重建和 invariant mass reconstruction 的分析思路。本节的任务是在这些工作的基础上，引入真实的靶、死层和探测器几何，使模拟结果不仅包含末态粒子的真值运动学信息，还包含粒子在材料中的能量损失、在探测器中的能量沉积以及实验读出形式，从而为后续的 PID、gate、Q 值重建和 invariant mass reconstruction 提供输入。
本实验关心的三体末态为 $$ {}^{10}\mathrm{Be}+\alpha+p . $$ 其中前向发射的 $^{10}\mathrm{Be}$ 和 $\alpha$ 主要由零度望远镜 $T_0$ 探测，反冲质子主要由侧向望远镜 $T_{1x}$、$T_{2x}$ 探测。因此，Geant4 模拟不能停留在 step 级信息上，而应最终整理为与实验分析兼容的事件级量。
2 反应机制与基本近似¶
本实验的核心反应道为 $$ {}^{1}\mathrm H(^{14}\mathrm C,\,^{14}\mathrm C^\ast){}^{1}\mathrm H,\qquad {}^{14}\mathrm C^\ast\rightarrow {}^{10}\mathrm{Be}^{(\ast)}+\alpha . $$ 第一步是非弹散射形成 $^{14}\mathrm C^\ast$，第二步是 $^{14}\mathrm C^\ast$ 的二体衰变。对于第一步反应，质心系散射角分布采用经验形式 $$ \frac{d\sigma}{d\Omega^\ast}\propto \exp\!\left(-\frac{\theta^\ast}{\alpha}\right),\qquad \alpha=53^\circ . $$ 对于第二步衰变，采用 $4\pi$ 各向同性二体衰变近似。$^{10}\mathrm{Be}$ 末态至少保留基态和 $3.37\ \mathrm{MeV}$ 激发态，必要时再加入更高激发分支。
这里需要强调，Geant4 内部并不具备本实验所需的专用反应机制。Geant4 可以处理粒子在材料中的输运，但不会自动产生“指定 $^{14}\mathrm C^\ast$ 激发态、指定分支比、指定角分布、再顺序衰变”的信号事件。因此，本节采用如下近似方法：
首先，在 PrimaryGeneratorAction 中由用户自定义信号反应发生器，生成该事件的真值四动量；然后，将生成的反冲质子、$\alpha$ 和 $^{10}\mathrm{Be}$ 作为 primary 交给 Geant4；最后，由 Geant4 处理这些粒子在靶、死层和探测器中的真实输运与能量沉积。
这一步最自然的实现方式，是直接复用 4.4 的 sequential two-body decay 代码。也就是说，在 PrimaryGeneratorAction.cc 中保留 4.4 已经建立的反应学生成逻辑，而不再额外重写一套新的衰变程序。这样既保持系列讲义的一致性，也避免重复实现已经验证过的运动学部分。
3 束流、靶和探测器几何¶
PPAC 用于确定束流粒子的靶前轨迹。在本节中，PPAC 不作为真实气体探测器详细模拟，而作为束流抽样输入，提供靶前 $x$、$y$ 位置与入射方向分布。位置分辨取 $$ \sigma_x=\sigma_y=0.65\ \mathrm{mm} . $$ 因此，Geant4 中的入射束流不应视为理想平行束，而应从 PPAC 重建结果中逐事件抽样。
实验使用多种靶条件，包括 $(\mathrm{CH}_2)_n$、$(\mathrm{CD}_2)_n$、C 靶和空靶。束流能量按高斯分布处理。若某种靶条件下束流均值为 $E_0$，FWHM 为 $\Delta E$，则高斯标准差应写为 $$ \sigma_E=\frac{\Delta E}{2.355} . $$ 在程序中，应使用 $\sigma_E$ 进行随机抽样，而不能直接把 FWHM 当作高斯标准差。
前向零度望远镜 $T_0$ 是本实验最关键的探测器。它由三层 BB7 DSSD、后续 SSD 和 $2\times2$ CsI(Tl) 组成。三层 BB7 的厚度均取 $1000\ \mu\mathrm m$，有效面积为 $64\times64\ \mathrm{mm}^2$，条带数为 $32+32$。这里尤其要注意，$T_0$ 第一层硅探测器 T0D1 的厚度必须严格保留为 $1000\ \mu\mathrm m$。这是因为低能粒子能否穿透第一层硅，将直接决定低能区域的探测效率和后续 PID 结构；若第一层厚度写错，后面的接受度、低能截止以及能谱结构都会发生系统偏差。
侧向望远镜 $T_1$、$T_2$ 用于探测反冲质子。其放置角度分别取约 $30.8^\circ$ 和 $68.5^\circ$，距靶距离约 $170\ \mathrm{mm}$。在 Geant4 中，应将这两个望远镜作为真实几何体摆放在相应方向上，而不是仅在分析阶段用一个理想角度窗口替代。
4 程序结构¶
建议程序按如下结构组织：
project/
├── CMakeLists.txt
├── include/
│ ├── DetectorConstruction.hh
│ ├── PrimaryGeneratorAction.hh
│ ├── PhysicsList.hh
│ ├── SensitiveDetector.hh
│ ├── EventAction.hh
│ └── RunAction.hh
├── src/
│ ├── main.cc
│ ├── DetectorConstruction.cc
│ ├── PrimaryGeneratorAction.cc
│ ├── PhysicsList.cc
│ ├── SensitiveDetector.cc
│ ├── EventAction.cc
│ └── RunAction.cc
└── data/
├── beam_ppac.root
└── mass.txt
其中，DetectorConstruction.cc 负责材料和几何；PrimaryGeneratorAction.cc 负责束流抽样、反应点抽样和反应学生成；PhysicsList.cc 负责输运过程；SensitiveDetector.cc 负责记录 step 级 hit；EventAction.cc 负责从 hit 整理出实验事件；RunAction.cc 负责整轮运行开始和结束时的文件处理。
5 DetectorConstruction.cc：材料与几何的实现¶
1. 材料定义¶
在 DetectorConstruction.cc 中，首先定义世界体、硅和 CsI 材料。常用材料可直接由 G4NistManager 获取。对于 $(\mathrm{CH}_2)_n$ 靶，则定义自定义材料：
auto nist = G4NistManager::Instance();
auto worldMat = nist->FindOrBuildMaterial("G4_Galactic");
auto Si = nist->FindOrBuildMaterial("G4_Si");
auto CsI = nist->FindOrBuildMaterial("G4_CESIUM_IODIDE");
G4double density = 0.93*g/cm3;
auto H = nist->FindOrBuildElement("H");
auto C = nist->FindOrBuildElement("C");
auto CH2 = new G4Material("Polyethylene_CH2", density, 2);
CH2->AddElement(C, 1);
CH2->AddElement(H, 2);
2. 靶的定义¶
靶应定义为一个薄片几何体，其厚度取实验实际靶厚。靶厚不仅影响 Geant4 中粒子的输运，也影响 PrimaryGeneratorAction 中反应点 $z$ 的抽样，因此这两个部分必须保持一致。
G4double targetXY = 30.0*mm;
G4double targetThickness = 100.0*um; // 改为实际靶厚
auto solidTarget =
new G4Box("Target",
targetXY/2,
targetXY/2,
targetThickness/2);
auto logicTarget =
new G4LogicalVolume(solidTarget, CH2, "TargetLV");
new G4PVPlacement(nullptr,
G4ThreeVector(0,0,0),
logicTarget,
"Target",
logicWorld,
false,
0);
3. $T_0$ 的定义¶
BB7 探测器可先定义为整块硅片，然后在 hit 整理阶段再由局域坐标映射为 strip 编号。这样做比把每一条 strip 都建成独立几何体更简洁，也更适合本节后续的数据整理方式。
G4double bb7XY = 64.0*mm;
G4double bb7Thickness = 1000.0*um;
auto solidBB7 =
new G4Box("BB7",
bb7XY/2,
bb7XY/2,
bb7Thickness/2);
auto logicBB7 =
new G4LogicalVolume(solidBB7, Si, "T0D_BB7_LV");
G4double zT0 = 170.0*mm;
for(int i=0; i<3; i++){
G4double z = zT0 + i*3.0*mm; // 层间距按实验结构修改
new G4PVPlacement(nullptr,
G4ThreeVector(0,0,z),
logicBB7,
"T0D_BB7",
logicWorld,
false,
i);
}
4. $T_1$、$T_2$ 的定义¶
侧向望远镜的中心方向由几何角度确定。对于 $T_1$，其中心位置可写为
G4double R = 170.0*mm;
G4double thetaT1 = 30.8*deg;
G4ThreeVector posT1(R*std::sin(thetaT1),
0,
R*std::cos(thetaT1));
auto rotT1 = new G4RotationMatrix();
rotT1->rotateY(thetaT1);
对于 $T_2$，将角度改为 68.5*deg 即可。DSSD、SSD 和 CsI 依次沿望远镜法向方向排布。
6 PrimaryGeneratorAction.cc：束流、反应点与反应学生成¶
1. 束流输入¶
束流信息从 beam_ppac.root 中读取。PrimaryGeneratorAction.hh 中应保存输入文件与输入树，并建立对应的 branch address。
#include "TFile.h"
#include "TTree.h"
class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction {
public:
PrimaryGeneratorAction();
virtual void GeneratePrimaries(G4Event*) override;
private:
TFile* fBeamFile;
TTree* fBeamTree;
double bx, by, btx, bty;
long nBeamEntries;
void SampleBeam(double& x, double& y,
double& tx, double& ty,
double& ekin);
};
构造函数中读取束流树：
fBeamFile = TFile::Open("beam_ppac.root");
fBeamTree = (TTree*)fBeamFile->Get("beam");
fBeamTree->SetBranchAddress("x", &bx);
fBeamTree->SetBranchAddress("y", &by);
fBeamTree->SetBranchAddress("tx", &btx);
fBeamTree->SetBranchAddress("ty", &bty);
nBeamEntries = fBeamTree->GetEntries();
逐事件抽样时，位置加入 PPAC 的位置分辨，能量按高斯分布抽样：
void PrimaryGeneratorAction::SampleBeam(double& x,
double& y,
double& tx,
double& ty,
double& ekin)
{
long id = G4RandFlat::shootInt(nBeamEntries);
fBeamTree->GetEntry(id);
x = bx + G4RandGauss::shoot(0.0, 0.65*mm);
y = by + G4RandGauss::shoot(0.0, 0.65*mm);
tx = btx;
ty = bty;
double mean = 317.1*MeV;
double fwhm = 1.98*MeV;
double sigma = fwhm/2.355;
ekin = G4RandGauss::shoot(mean, sigma);
}
其中，tx 和 ty 用小角近似理解为
$$
t_x=\frac{p_x}{p_z},\qquad t_y=\frac{p_y}{p_z},
$$
束流方向单位矢量可写为
G4ThreeVector dir(tx, ty, 1.0);
dir = dir.unit();
2. 反应点抽样¶
反应点不应固定在靶中心，而应在靶厚范围内均匀抽样。若反应点位于 $z_{\mathrm{vertex}}$，则入射束流在到达该位置前已经在靶中走过一段路径，其有效入射能量应相应降低。
设入射能量为 $E_{\mathrm{in}}$，入射前在靶中走过的路径长度为 $\ell_{\mathrm{in}}$，则反应点处的能量写为 $$ E_{\mathrm{reaction}} = E_{\mathrm{in}} - \left(\frac{dE}{dx}\right)\ell_{\mathrm{in}} . $$
为了避免重复实现已有过程，停能本领应尽量调用 Geant4 内部机制。可使用 G4EmCalculator 查询离子在靶材料中的总能量损失率：
G4EmCalculator emCal;
auto ionTable = G4IonTable::GetIonTable();
auto C14 = ionTable->GetIon(6,14,0.0);
double dedx =
emCal.ComputeTotalDEDX(Ein, C14, targetMaterial);
double path =
(zVertex + targetThickness/2.0)/std::cos(thetaBeam);
double Ereac = Ein - dedx*path;
3. 第一步反应的运动学¶
第一步反应为 $$ {}^{14}\mathrm C + p \rightarrow {}^{14}\mathrm C^\ast + p . $$
若入射 $^{14}\mathrm C$ 质量为 $m_a$，靶质子质量为 $m_A$，出射 $^{14}\mathrm C^\ast$ 质量为 $m_b=m_{14C}+E_x$，反冲质子质量为 $m_B$，则入射总能量为 $$ E_a=T_a+m_a , $$ 入射动量大小为 $$ p_a=\sqrt{E_a^2-m_a^2} . $$
程序中可写为
double EaTot = Ereac + m14C;
double pa = std::sqrt(EaTot*EaTot - m14C*m14C);
G4ThreeVector beamDir(tx, ty, 1.0);
beamDir = beamDir.unit();
设 Mandelstam 变量为 $s$，则质心系中末态动量大小为 $$ p^\ast = \frac{\sqrt{[s-(m_b+m_B)^2][s-(m_b-m_B)^2]}}{2\sqrt{s}} . $$
程序中为
double mC14star = m14C + Ex;
double mOutP = mp;
double pcm =
std::sqrt((s - std::pow(mC14star + mOutP,2)) *
(s - std::pow(mC14star - mOutP,2)))
/(2.0*sqrtS);
散射角采用经验分布抽样。实现时可以用接受拒绝法。若抽样函数权重取 $$ w(\theta^\ast)=\sin\theta^\ast\exp\!\left(-\frac{\theta^\ast}{\alpha}\right), $$ 则程序可写为
double weight = std::sin(theta)*std::exp(-theta/alpha);
if(y < weight) return theta;
随后对 $\phi^\ast$ 取均匀分布：
double theta = SampleThetaCM();
double phi = G4RandFlat::shoot(0.0, twopi);
CLHEP::Hep3Vector pStarVec(pcm*std::sin(theta)*std::cos(phi),
pcm*std::sin(theta)*std::sin(phi),
pcm*std::cos(theta));
4. 激发态与 sequential decay¶
本节中，$^{14}\mathrm C^\ast$ 的激发态与 $^{10}\mathrm{Be}$ 的末态分支仍按 4.4 的方式处理。若某核处于激发态，则其参与运动学的质量应写为
$$
m^\ast = m_{\rm gs}+E_x .
$$
在程序中，建议保留一组 Be10Branch 与 C14Level 结构体，先抽取 $^{10}\mathrm{Be}$ 分支，再抽取相应的 $^{14}\mathrm C^\ast$ 激发能，然后直接调用 4.4 中已经建立的 sequential two-body decay 逻辑，生成末态 $p$、$\alpha$ 和 $^{10}\mathrm{Be}$ 的四动量。
因此，PrimaryGeneratorAction.cc 中应至少包含以下函数：
void SampleBeam(double& x, double& y,
double& tx, double& ty,
double& ekin);
void SampleVertex(double& vx, double& vy, double& vz);
double SampleEx14C(...);
double SampleEx10Be(...);
void GenerateReactionKinematics(...);
最后，将生成出的末态粒子转换为 Geant4 primary，并挂到同一个 G4PrimaryVertex 上：
auto vertex = new G4PrimaryVertex(G4ThreeVector(vx, vy, vz), 0.0);
auto pDef = G4Proton::Definition();
auto alphaDef = G4Alpha::Definition();
auto beDef = G4IonTable::GetIonTable()->GetIon(4,10,0.0);
auto pPrim = new G4PrimaryParticle(pDef, px_p, py_p, pz_p);
auto aPrim = new G4PrimaryParticle(alphaDef, px_a, py_a, pz_a);
auto bePrim = new G4PrimaryParticle(beDef, px_be, py_be, pz_be);
vertex->SetPrimary(pPrim);
vertex->SetPrimary(aPrim);
vertex->SetPrimary(bePrim);
event->AddPrimaryVertex(vertex);
在这里，不再手工实现靶中能损、多重散射或死层修正。这些过程全部交由 Geant4 的输运模型处理。
7 PhysicsList.cc：输运过程的选取¶
本节的重点不是由 Geant4 自动生成信号反应，而是让 Geant4 正确处理末态粒子在材料中的输运。因此，应尽量使用 Geant4 内部已有、经过验证的电磁与输运机制，避免自行重复实现 ionisation、multiple scattering 或 stopping power。
推荐采用 G4EmStandardPhysics_option4 作为低能电磁过程的核心配置。若使用 modular physics list，可写为
RegisterPhysics(new G4EmStandardPhysics_option4());
RegisterPhysics(new G4DecayPhysics());
若使用 reference physics list，也可写为
auto physicsList = new FTFP_BERT;
physicsList->ReplacePhysics(new G4EmStandardPhysics_option4());
对于本实验中的质子、$\alpha$ 和 $^{10}\mathrm{Be}$，需要重点保留的物理过程包括：
质子的 ionisation、multiple scattering 和 elastic scattering；$\alpha$ 与重离子的 ionisation、multiple scattering 以及相应的 nuclear stopping。这样处理之后，粒子在靶、死层、硅和 CsI 中的能量沉积将由 Geant4 自动给出，从而能够自然反映几何厚度与材料结构对实验响应的影响。
8 SensitiveDetector.cc：step 级 hit 的记录¶
SensitiveDetector 负责记录 Geant4 的原始 hit 信息，而不直接构造实验事件。对于每个 hit，至少应保留以下量：
探测器编号 detID，层编号 layerID，能量沉积 edep，粒子编号 trackID，粒子种类 pdg，全局位置 posGlobal，局域位置 posLocal，条带编号 stripX、stripY，以及 step 前后动能 preKinE、postKinE。
实现时，先取得该 step 的总能量沉积，若 edep<=0 则直接返回；否则从 PreStepPoint 和 Touchable 中读取位置和 copy number，并将结果写入 hit：
G4bool SiSD::ProcessHits(G4Step* step,
G4TouchableHistory*)
{
double edep = step->GetTotalEnergyDeposit();
if(edep <= 0) return false;
auto pre = step->GetPreStepPoint();
auto post = step->GetPostStepPoint();
auto touchable = pre->GetTouchableHandle();
int copyNo = touchable->GetCopyNumber();
SiHit* hit = new SiHit;
hit->detID = copyNo;
hit->edep = edep;
hit->trackID = step->GetTrack()->GetTrackID();
hit->pdg = step->GetTrack()->GetDefinition()->GetPDGEncoding();
hit->preKinE = pre->GetKineticEnergy();
hit->postKinE = post->GetKineticEnergy();
hitsCollection->insert(hit);
return true;
}
如果希望后续按实验读出形式整理，则在这里还应保留局域坐标，以便后续映射为 strip 编号。
9 EventAction.cc：从 hit 到实验事件¶
EventAction 负责把 SensitiveDetector 记录下来的 step 级 hit 整理为实验事件级数据。这一步包括四个连续过程：同一通道内多 step 合并、连续位置离散为 strip 读出、能量分辨与 threshold 处理、按望远镜组织实验事件。
1. 通道合并¶
由于同一粒子在同一条带、同一层中可能产生多个 Geant4 step，因此不能把每个 step 都视为一个独立读出。应按探测器、层号和条带号进行聚合。可以定义
struct ChannelID {
int detID;
int layerID;
int stripX;
int stripY;
bool operator<(const ChannelID& other) const
{
if(detID != other.detID) return detID < other.detID;
if(layerID!= other.layerID) return layerID< other.layerID;
if(stripX != other.stripX) return stripX < other.stripX;
return stripY < other.stripY;
}
};
随后用 std::map<ChannelID,double> 或 std::map<ChannelID,ChannelSignal> 对能量沉积进行累加：
std::map<ChannelID, double> edepMap;
for(auto hit : hits){
ChannelID ch{hit->detID,
hit->layerID,
hit->stripX,
hit->stripY};
edepMap[ch] += hit->edep;
}
若需要同时保留位置、时间和多重性，则可定义
struct ChannelSignal {
double edep = 0.0;
double time = 1e99;
double x = 0.0;
double y = 0.0;
double z = 0.0;
int multiplicity = 0;
};
2. strip 映射¶
Geant4 中的 hit 位置是连续量，而实验读出是离散的 strip 编号。因此在 EventAction 中，需要将局域坐标映射为 strip 号，并在需要时再换算成 strip 中心位置。
对于 BB7，条带间距为 $$ p=\frac{64\ \mathrm{mm}}{32}=2\ \mathrm{mm} . $$ 若条带编号为 $n_x$、$n_y$，则条带中心位置为 $$ x_{\mathrm{strip}}=-32\ \mathrm{mm}+\left(n_x+\frac12\right)p , $$ $$ y_{\mathrm{strip}}=-32\ \mathrm{mm}+\left(n_y+\frac12\right)p . $$
程序可写为
double GetStripCenterX(int stripX)
{
G4double pitch = 64.0*mm / 32.0;
return -32.0*mm + (stripX + 0.5)*pitch;
}
double GetStripCenterY(int stripY)
{
G4double pitch = 64.0*mm / 32.0;
return -32.0*mm + (stripY + 0.5)*pitch;
}
3. 能量分辨与 threshold¶
实验读出中，DSSD、SSD 和 CsI 的能量测量都具有有限分辨；同时，电子学读出还存在触发阈值。因此，在通道合并之后，应对每个通道施加能量展宽和 threshold。
一个简单的处理方式为
double SmearEnergy(double E, int detType)
{
double sigma = 0.0;
if(detType == kDSSD){
sigma = 0.03*MeV;
}
else if(detType == kCsI){
sigma = 0.02*E;
}
return G4RandGauss::shoot(E, sigma);
}
然后在写入事件之前施加阈值：
double threshold = 0.2*MeV;
if(E_smeared < threshold) continue;
这一步应放在事件整理层，而不放在 SensitiveDetector 中。因为 step 级 hit 记录的任务是尽可能完整地保存输运结果，而实验读出形式的形成属于事件级处理。
4. 望远镜级事件组织¶
在完成通道合并、strip 离散化、展宽和 threshold 后，应将信号组织为望远镜级读出。对于 $T_0$，至少需要整理出三层 DSSD、后续 SSD 和 CsI 的能量，以及第一层条带编号和多重性；对于 $T_1$、$T_2$，至少需要整理出 DSSD、SSD 和 CsI 的能量、条带编号和多重性。
此外，还应在事件级形成基本符合标志，例如：
- 是否 hit $T_0$
- 是否 hit $T_1$
- 是否 hit $T_2$
- 是否形成 triple coincidence
这些量将在后续 ROOT 输出与离线分析中直接使用，因此应在 EventAction.cc 中完成整理。
10 main.cc、编译与运行¶
主程序中，将各个用户类依次挂到 RunManager 上：
int main(int argc, char** argv)
{
auto* runManager = new G4RunManager;
runManager->SetUserInitialization(new DetectorConstruction);
runManager->SetUserInitialization(new PhysicsList);
runManager->SetUserAction(new PrimaryGeneratorAction);
runManager->SetUserAction(new EventAction);
runManager->SetUserAction(new RunAction);
runManager->Initialize();
G4UImanager* UImanager = G4UImanager::GetUIpointer();
if (argc > 1) {
G4String command = "/control/execute ";
G4String fileName = argv[1];
UImanager->ApplyCommand(command + fileName);
}
delete runManager;
return 0;
}
CMakeLists.txt 可写为
cmake_minimum_required(VERSION 3.16)
project(g4sim)
find_package(Geant4 REQUIRED ui_all vis_all)
find_package(ROOT REQUIRED)
include(${Geant4_USE_FILE})
include_directories(${PROJECT_SOURCE_DIR}/include)
include_directories(${ROOT_INCLUDE_DIRS})
file(GLOB SOURCES src/*.cc)
file(GLOB HEADERS include/*.hh)
add_executable(g4sim ${SOURCES} ${HEADERS})
target_link_libraries(g4sim ${Geant4_LIBRARIES} ${ROOT_LIBRARIES})
编译过程为
mkdir build
cd build
cmake ..
make -j4
准备运行宏文件 run.mac：
/run/initialize
/run/beamOn 100000
随后执行
./g4sim run.mac
程序运行结束后，将生成供后续分析使用的事件 ROOT 文件。
11 RunAction.cc 与 ROOT 文件输出¶
Geant4 事件在 EventAction.cc 中已经被整理为实验事件级量。接下来的任务，是将这些量写入 ROOT 文件，供后续离线分析直接读取。这里不输出 step 级原始信息，而输出与实验分析兼容的事件树。
程序中可将 ROOT 文件的创建、树的建立和分支定义集中在一个单独的类中，例如 RootIO，再由 RunAction 在每一轮运行开始和结束时负责打开、写入和关闭文件。
RunAction.cc 的基本结构可写为
void RunAction::BeginOfRunAction(const G4Run*)
{
RootIO::Instance()->OpenFile("g4sim.root");
}
void RunAction::EndOfRunAction(const G4Run*)
{
RootIO::Instance()->Write();
RootIO::Instance()->Close();
}
这里输出文件名取为 g4sim.root。若需要区分不同靶条件、不同激发态方案或不同几何设置，可以在文件名中加入相应标签，例如 g4sim_ch2.root、g4sim_cd2.root 或 g4sim_ex.root。
12 ROOT 树的组织方式¶
ROOT 输出应尽量贴近实验分析所使用的变量，而不是直接暴露 Geant4 内部 step 级细节。因此，树中至少应包含以下三类信息：
第一类是真值信息，用于后续将模拟结果与真实输入反应学进行对照；
第二类是望远镜级实验读出信息，用于直接构造 $\Delta E-E$、gate、Q 值与 invariant mass reconstruction；
第三类是符合标志，用于快速筛选事件类型。
树名可定义为 evt。在 RootIO::OpenFile() 中可写为
void RootIO::OpenFile(const std::string& filename)
{
fFile = new TFile(filename.c_str(), "RECREATE");
fTree = new TTree("evt", "Geant4 simulated events");
}
1. 真值分支¶
真值分支用于保存反应点、激发态以及末态粒子的真值四动量。建议至少包含：
- 反应点：
vx,vy,vz - 激发能：
ex14c,ex10be - 反冲质子真值四动量：
p_px,p_py,p_pz,p_E - $\alpha$ 真值四动量：
a_px,a_py,a_pz,a_E - $^{10}\mathrm{Be}$ 真值四动量：
be_px,be_py,be_pz,be_E
分支定义可写为
fTree->Branch("vx", &vx, "vx/D");
fTree->Branch("vy", &vy, "vy/D");
fTree->Branch("vz", &vz, "vz/D");
fTree->Branch("ex14c", &ex14c, "ex14c/D");
fTree->Branch("ex10be", &ex10be, "ex10be/D");
fTree->Branch("p_px", &p_px, "p_px/D");
fTree->Branch("p_py", &p_py, "p_py/D");
fTree->Branch("p_pz", &p_pz, "p_pz/D");
fTree->Branch("p_E", &p_E, "p_E/D");
fTree->Branch("a_px", &a_px, "a_px/D");
fTree->Branch("a_py", &a_py, "a_py/D");
fTree->Branch("a_pz", &a_pz, "a_pz/D");
fTree->Branch("a_E", &a_E, "a_E/D");
fTree->Branch("be_px", &be_px, "be_px/D");
fTree->Branch("be_py", &be_py, "be_py/D");
fTree->Branch("be_pz", &be_pz, "be_pz/D");
fTree->Branch("be_E", &be_E, "be_E/D");
这些量主要用于理解模拟输入、检查不同激发态分支的重建效果，以及在需要时和实验分析结果进行对照。
2. $T_0$ 的实验读出分支¶
对于零度望远镜 $T_0$，建议保留三层 DSSD、SSD 和 CsI 的能量读出，以及第一层 DSSD 的条带编号和多重性。这样即可直接构造前向望远镜的 $\Delta E-E$ 图和位置读出。
建议包含：
t0_E_D1t0_E_D2t0_E_D3t0_E_SSDt0_E_CsIt0_StripX_D1t0_StripY_D1t0_mult_D1t0_mult_D2t0_mult_D3
代码可写为
fTree->Branch("t0_E_D1", &t0_E_D1, "t0_E_D1/D");
fTree->Branch("t0_E_D2", &t0_E_D2, "t0_E_D2/D");
fTree->Branch("t0_E_D3", &t0_E_D3, "t0_E_D3/D");
fTree->Branch("t0_E_SSD", &t0_E_SSD, "t0_E_SSD/D");
fTree->Branch("t0_E_CsI", &t0_E_CsI, "t0_E_CsI/D");
fTree->Branch("t0_StripX_D1", &t0_StripX_D1, "t0_StripX_D1/I");
fTree->Branch("t0_StripY_D1", &t0_StripY_D1, "t0_StripY_D1/I");
fTree->Branch("t0_mult_D1", &t0_mult_D1, "t0_mult_D1/I");
fTree->Branch("t0_mult_D2", &t0_mult_D2, "t0_mult_D2/I");
fTree->Branch("t0_mult_D3", &t0_mult_D3, "t0_mult_D3/I");
3. $T_1$、$T_2$ 的实验读出分支¶
侧向望远镜主要探测反冲质子，因此对于每个侧向望远镜，至少保留 DSSD、SSD、CsI 的能量以及 DSSD 条带位置和多重性。
对于 $T_1$：
fTree->Branch("t1_E_DSSD", &t1_E_DSSD, "t1_E_DSSD/D");
fTree->Branch("t1_E_SSD", &t1_E_SSD, "t1_E_SSD/D");
fTree->Branch("t1_E_CsI", &t1_E_CsI, "t1_E_CsI/D");
fTree->Branch("t1_StripX", &t1_StripX, "t1_StripX/I");
fTree->Branch("t1_StripY", &t1_StripY, "t1_StripY/I");
fTree->Branch("t1_mult", &t1_mult, "t1_mult/I");
对于 $T_2$：
fTree->Branch("t2_E_DSSD", &t2_E_DSSD, "t2_E_DSSD/D");
fTree->Branch("t2_E_SSD", &t2_E_SSD, "t2_E_SSD/D");
fTree->Branch("t2_E_CsI", &t2_E_CsI, "t2_E_CsI/D");
fTree->Branch("t2_StripX", &t2_StripX, "t2_StripX/I");
fTree->Branch("t2_StripY", &t2_StripY, "t2_StripY/I");
fTree->Branch("t2_mult", &t2_mult, "t2_mult/I");
4. 符合标志分支¶
为了在分析中快速筛选事件，应同时保存事件级的符合标志：
flag_hitT0flag_hitT1flag_hitT2flag_triple
其中，flag_triple 表示 $T_0$、$T_1$、$T_2$ 三个望远镜均有有效信号。
fTree->Branch("flag_hitT0", &flag_hitT0, "flag_hitT0/I");
fTree->Branch("flag_hitT1", &flag_hitT1, "flag_hitT1/I");
fTree->Branch("flag_hitT2", &flag_hitT2, "flag_hitT2/I");
fTree->Branch("flag_triple", &flag_triple, "flag_triple/I");
13 事件写入¶
在 EventAction.cc 中，当事件整理完成后，将对应变量写入 RootIO 的缓存成员，然后调用 Fill()。例如
void EventAction::EndOfEventAction(const G4Event*)
{
auto* io = RootIO::Instance();
io->vx = vx;
io->vy = vy;
io->vz = vz;
io->ex14c = ex14c;
io->ex10be = ex10be;
io->t0_E_D1 = t0_E_D1;
io->t0_E_D2 = t0_E_D2;
io->t0_E_D3 = t0_E_D3;
io->t0_E_SSD = t0_E_SSD;
io->t0_E_CsI = t0_E_CsI;
io->t1_E_DSSD = t1_E_DSSD;
io->t1_E_SSD = t1_E_SSD;
io->t1_E_CsI = t1_E_CsI;
io->t2_E_DSSD = t2_E_DSSD;
io->t2_E_SSD = t2_E_SSD;
io->t2_E_CsI = t2_E_CsI;
io->flag_hitT0 = flag_hitT0;
io->flag_hitT1 = flag_hitT1;
io->flag_hitT2 = flag_hitT2;
io->flag_triple = flag_triple;
io->Fill();
}
这里的关键点是：事件整理完成后再统一写入 ROOT，而不是在 SensitiveDetector 中逐个 step 写入。这样输出树中的每一行才对应一个实验事件，而不是一个 Geant4 step。
4.7.14 ROOT 文件的打开方式¶
程序运行结束后，将得到 ROOT 文件，例如
g4sim.root
打开 ROOT 文件的方法为
root -l g4sim.root
进入 ROOT 之后，首先读取文件和树：
TFile *f = TFile::Open("g4sim.root");
TTree *tr = (TTree*)f->Get("evt");
tr->Print();
这样即可查看树中包含的分支。
如果希望直接在 ROOT 命令行中检查事件数，可写为
tr->GetEntries();
14 模拟数据的后续分析¶
ROOT 文件生成之后，相应的数据分析按照 4.5 和 4.6 的相应部分进行。也就是说，对模拟数据施加和实际数据分析完全相同的条件、gate 和处理步骤，使模拟结果能够直接与实际数据进行对应比较。
具体而言，首先按照实际数据分析中的粒子鉴别方法，对 $T_0$、$T_1$ 和 $T_2$ 的探测器信号构造相同的 $\Delta E-E$ 变量，并施加与实际数据一致的 PID gate，选出 $\alpha$、$^{10}\mathrm{Be}$ 和反冲质子事件。随后，按照 4.5 中的处理流程，对选出的事件进行反应 Q 值重建，以区分不同的 $^{10}\mathrm{Be}$ 末态分支，并建立与实际数据完全一致的 Q 值条件。
在完成 Q 值选择之后，再按照 4.6 中 invariant mass reconstruction 的相应步骤，对 $\alpha+{}^{10}\mathrm{Be}$ 末态进行不变质量重建。这里所使用的事件选择、粒子组合方式、门选条件和重建顺序，都与实际数据分析保持一致。这样得到的模拟谱线可以直接与实验谱线进行比较，用于判断几何接受度、探测器响应、能量损失和门选条件对最终谱形的影响。
采用这种处理方式的目的是使 Geant4 模拟真正进入实验分析链，而不是停留在独立的数值模拟层面。通过施加和实际数据分析完全相同的条件和步骤，模拟结果中的每一个结构都可以在与实验相同的分析框架下进行理解，从而用于实验设计优化、探测器响应评估以及谱结构解释。