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

课程主页 · 上一节：4.6

基于 Geant4 的实验模拟

1 本节目标

4.4 给出了反应点处的粒子四动量，4.5–4.6 用简化读数练习重建。本节沿同一反应继续讨论：粒子经过靶和探测器后，怎样由能量沉积形成可分析的事件。

保持以下分工：反应发生器生成末态；Geant4 处理材料输运；读出模型把沉积转换成通道信号；离线分析完成 PID、粒子关联、能损修正、Q 与不变质量重建。真值只作检查，不替代测量输入。

下面各小节给出模块接口和关键代码片段，不是一份已齐备的 Geant4 工程。诸如 SiHit、RootIO 是应用程序需要实现的类，不是 Geant4 自动提供的功能；beam_ppac.root 和完整几何参数也需实验输入。

2 反应机制与基本近似

反应仍为 $^{14}\mathrm C+p\to{}^{14}\mathrm C^*+p$，随后 $^{14}\mathrm C^*\to\alpha+{}^{10}\mathrm{Be}^{(*)}$。第一步是非弹散射，第二步是激发核衰变。可先沿用 4.4 的各向同性两体抽样，检查分析链；研究接受度对角分布的依赖时，再给第一步指定分布。例如 $$\frac{d\sigma}{d\Omega^*}\propto\exp(-\theta^*/\alpha).$$ 下面以 $\alpha=53^\circ$ 为演示参数，不将其当作由所引文献确定的实验截面。第二步仍采用各向同性近似，不包含自旋关联。

几何背景参见 Han et al., Communications Physics 6, 220 (2023), Methods。

这里需要强调，Geant4 内部并不具备本实验所需的专用反应机制。Geant4 可以处理粒子在材料中的输运，但不会自动产生“指定 $^{14}\mathrm C^\ast$ 激发态、指定分支比、指定角分布、再顺序衰变”的信号事件。因此，本节采用如下近似方法：

首先，在 PrimaryGeneratorAction 中由用户自定义信号反应发生器，生成该事件的真值四动量；然后，将生成的反冲质子、$\alpha$ 和 $^{10}\mathrm{Be}$ 作为 primary 交给 Geant4；最后，由 Geant4 处理这些粒子在靶、死层和探测器中的真实输运与能量沉积。

可复用 4.4 中“给定初态与各核质量，生成一次级联衰变”的部分。不要把 4.4 的整个事件生成循环 连同读出展宽一起搬入 PrimaryGeneratorAction：本节的能量沉积已经由 Geant4 产生，读出展宽应只在后续通道整理时施加一次。若改用了经验散射角抽样，就替换第一步各向同性方向，不要再叠加一套互不相容的方向抽样。

---

3 束流、靶和探测器几何

PPAC 数据提供束流位置与方向的相关分布。若输入树中保存的是已重建的实验轨迹，抽样这些记录已经包含测量分辨，不再额外叠加同一 PPAC 分辨。若从理想束流真值生成 PPAC 读数，才需要加入独立规定的测量误差；例如可用 $\sigma_x=\sigma_y=0.65$ mm 作教学假设，但本节不把它当作已验证的实测参数。

本例的 x,y 定义在靶中心平面，单位 mm；tx,ty 是无量纲斜率。读取实验文件时先核对其参考平面和单位。保持同一条轨迹的四个量一起抽样，才能保留束斑与发散角的相关性。

实验使用多种靶条件，包括 $(\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$ 包含三层 1000 μm 的 BB7 DSSD、三层 SSD 和 $2\times2$ CsI(Tl)。下面仅示范三层 DSSD 的构造：有效面积取 $64\times64$ mm²、每面 32 条，条带间距 2 mm。层间距、死层及后级 SSD/CsI 的具体尺寸必须从实验几何补齐，不能据这段简化代码计算整套装置的效率。

文献中侧向望远镜中心角约为 $30.8^\circ$ 和 $68.5^\circ$，距靶约 170 mm。$T_{1x}$ 包含 50 μm W1 与 300 μm BB7 两层 DSSD，$T_{2x}$ 的 DSSD 为 65 μm BB7；二者均有后续探测器。这里的 $T_1,T_2$ 只表示两类侧向位置，不应误认为完整实验只有两个侧向模块。

---

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"
#include <stdexcept>

class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction {
public:
    PrimaryGeneratorAction();
    virtual void GeneratePrimaries(G4Event*) override;

private:
    TFile* fBeamFile = nullptr;
    TTree* fBeamTree = nullptr;

    double bx, by, btx, bty;
    Long64_t nBeamEntries = 0; // 与 Tree 的 entry 类型一致

    void SampleBeam(double& x, double& y,
                    double& tx, double& ty,
                    double& ekin);
};

构造函数中读取束流树：

fBeamFile = TFile::Open("data/beam_ppac.root");
if (!fBeamFile || fBeamFile->IsZombie())
    throw std::runtime_error("Cannot open data/beam_ppac.root");
fBeamTree = fBeamFile->Get<TTree>("beam");
if (!fBeamTree || fBeamTree->GetEntries()==0)
    throw std::runtime_error("Missing or empty beam Tree");

fBeamTree->SetBranchAddress("x",  &bx);
fBeamTree->SetBranchAddress("y",  &by);
fBeamTree->SetBranchAddress("tx", &btx);
fBeamTree->SetBranchAddress("ty", &bty);

nBeamEntries = fBeamTree->GetEntries();

逐事件读取同一条束流记录，位置乘以单位 mm，能量沿用本练习的 Gaussian 输入：

void PrimaryGeneratorAction::SampleBeam(double& x,
                                        double& y,
                                        double& tx,
                                        double& ty,
                                        double& ekin)
{
    long id = G4RandFlat::shootInt(nBeamEntries);
    fBeamTree->GetEntry(id);

    x  = bx*mm; // 已测得的位置，不重复叠加 PPAC 分辨
    y  = by*mm;

    tx = btx;
    ty = bty;

    double mean  = 317.1*MeV;
    double fwhm  = 1.83*MeV; // 与 4.4 的教学输入一致
    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$。厚靶反应概率可能随能量改变，不能总作均匀假设。给定一条在 $z=0$ 平面记录的轨迹，反应点应满足 $x_v=x_0+t_xz_v$、$y_v=y_0+t_yz_v$，并检查仍位于靶内。

设入射能量为 $E_{\mathrm{in}}$，入射前在靶中走过的路径长度为 $\ell_{\mathrm{in}}$，则反应点处的能量写为 $$ E_{\mathrm{reaction}} = E_{\mathrm{in}} - \left(\frac{dE}{dx}\right)\ell_{\mathrm{in}} . $$

下面用 G4EmCalculator 查询已初始化物理表中的平均阻止本领，作薄靶的一步近似；若沿路径能量变化明显，应分段更新 $dE/dx$ 或实际输运入射束流。这段估计不产生能量歧离。若入射束流已被 Geant4 跟踪到反应点，就直接使用该处动能，不再重复扣除这段能损。

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();

// 以下质量、动能和动量统一用 Geant4 能量单位。
double s = m14C*m14C + mp*mp + 2.0*mp*EaTot;
double sqrtS = std::sqrt(s);

设 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;

if (sqrtS <= mC14star + mOutP) return; // 此候选事件不允许
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 SampleThetaCM()
{
    const double alpha = 53.0*deg;
    while (true) {
        double theta = G4RandFlat::shoot(0.0, CLHEP::pi);
        double w = std::sin(theta)*std::exp(-theta/alpha);
        if (G4RandFlat::shoot() < w) return theta; // w <= 1
    }
}

随后对 $\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));

// 这里 theta 是 14C* 在质心系中相对束流的极角。
pStarVec.rotateUz(beamDir);
G4LorentzVector starCM(pStarVec,
    std::sqrt(pcm*pcm + mC14star*mC14star));
G4LorentzVector recoilCM(-pStarVec,
    std::sqrt(pcm*pcm + mp*mp));
G4ThreeVector betaCM = pa*beamDir/(EaTot + mp);
starCM.boost(betaCM);
recoilCM.boost(betaCM); // 两者现在均在实验室系

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 一致。若复用 4.4 的 GeV 数值，先乘 GeV；不要将 ROOT 中裸数值的 GeV 当成 Geant4 的 MeV。生成器与输运粒子还应采用一致的核质量。若生成器已让 $^{10}\mathrm{Be}^*$ 退激，则交给 Geant4 的是基态 Be，并按模拟范围决定是否输运同时生成的 γ；不能再让同一 Be 重复退激。

出射粒子从反应点向外的能损由 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());

第一种配置用于研究电离能损和多重散射等电磁效应，不包含完整的强子反应。若要研究核弹性散射、非弹性反应造成的损失或本底，需要第二类包含强子过程的配置，并核对质子、α、重离子在相应能区的模型。两种配置的物理范围不同，不是可任意互换的写法。过程清单和能区应查阅 Geant4 Physics List Guide。

质子的 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。

在 ConstructSDandField() 中创建并注册 SiSD，将它关联到硅的 logical volume；在 SiSD::Initialize() 中为每个事件创建并登记 hitsCollection。这样有沉积的 step 才会调用下面的 ProcessHits()，且不同事件不会共用上一事件的 hit。接口见 Geant4：Hits。

实现时，先读取该 step 的总能量沉积；若 edep<=0 则返回。下面将沉积位置近似取为 step 中点，并转换到探测器局域坐标：

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   = 0; // 这里只示范 T0 的三层 DSSD
    hit->layerID = copyNo; // 与上面的 placement 编号 0,1,2 对应
    hit->posGlobal = 0.5*(pre->GetPosition() + post->GetPosition());
    hit->posLocal = touchable->GetHistory()->GetTopTransform()
                             .TransformPoint(hit->posGlobal);
    hit->stripX = int(std::floor((hit->posLocal.x()/mm + 32.0)/2.0));
    hit->stripY = int(std::floor((hit->posLocal.y()/mm + 32.0)/2.0));
    if (hit->stripX < 0 || hit->stripX >= 32 ||
        hit->stripY < 0 || hit->stripY >= 32) {
        delete hit;
        return false;
    }
    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;
}

这段条带映射把一个 step 的全部沉积归到中点所在条带。它适用于 step 横向长度远小于条带宽度的近似；跨越多个条带时应细分沉积或使用条带几何，否则会误计多重性。电荷共享属于后续读出模型，不由这个整数编号自动产生。

---

9 EventAction.cc：从 hit 到实验事件

EventAction 先把局域位置映射到条带，再合并同一通道的多个 step，最后施加一次读出分辨与阈值。为便于理解，下面先说明通道合并的目标，再给出位置到条带的关系。

1. 通道合并

由于同一粒子在同一条带、同一层中可能产生多个 Geant4 step，因此不能把每个 step 都视为一个独立读出。应按探测器、层号和条带号进行聚合。可以定义

struct ChannelID {
    int detID, layerID;
    int side;   // 0: front, 1: back
    int strip;  // 本面条带号
    bool operator<(const ChannelID& other) const {
        return std::tie(detID, layerID, side, strip)
             < std::tie(other.detID, other.layerID, other.side, other.strip);
    }
}; // 需要 #include <tuple> 和 <map>

随后用 std::map<ChannelID,double> 或 std::map<ChannelID,ChannelSignal> 对能量沉积进行累加：

std::map<ChannelID, double> edepMap;
for (auto hit : hits) {
    // stripX/Y 已由局域坐标得到，先确认在 [0,31] 内。
    ChannelID front{hit->detID, hit->layerID, 0, hit->stripX};
    ChannelID back {hit->detID, hit->layerID, 1, hit->stripY};
    edepMap[front] += hit->edep;
    edepMap[back]  += hit->edep;
}
// 两面是同一次沉积的两套读出；求粒子能量时不能把 front+back 相加。

DSSD 是两面各一组条带读出，不是 stripX,stripY 定义的逐像素读出。上例在两面分别累加同一能量沉积。std::tie 只是按四个编号排序，供 std::map 区分通道；它不做任何物理选择。若还需记录通道时间与位置，可保存：

struct ChannelSignal {
    double edep = 0.0;
    double time = 1e99;
    double x = 0.0;
    double y = 0.0;
    double z = 0.0;
    int nSteps = 0; // 累加的 step 数，不是触发条带数或粒子数
};

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。

下面的 30 keV 和 2% 是读出模型的演示参数，不是实测标定；kDSSD/kSSD/kCsI 是应用程序定义的探测器类型。CsI 的真实光输出可能依赖粒子种类，不能普遍用沉积能量的固定比例描述：

double SmearEnergy(double E, int detType)
{
    double sigma = 0.0;

    if(detType == kDSSD || detType == kSSD){
        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;

每个通道先累加沉积，再展宽一次并应用读出阈值。不要对每个 step 分别展宽。条带多重性是阈值后有信号的条带数，不是 step 数或 track 数。Geant4 已产生的能损歧离不能作为同一效应再额外叠加一次。

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})

编译过程为

# 以下命令均在工程根目录运行
cmake -S . -B build
cmake --build build -j4

准备运行宏文件 run.mac：

/run/initialize
/run/beamOn 100000

随后执行

./build/g4sim run.mac

上述命令适用于已补齐用户类、数据路径与完整几何的工程。在只有本页片段的情况下不能直接编译运行；这也不代表已经生成了可验证的实验模拟结果。

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 树的组织方式

输出同时保留真值、通道读数和选择标志；不必把全部 step 写入分析树。下面先展示标量分支的写法，完整多粒子分析还需要后述的逐通道信息，不能仅用望远镜总能量直接计算 Q 值。能量与动量输出分别统一用 MeV、MeV/$c$，位置用 mm。

第一类是真值信息，用于后续将模拟结果与真实输入反应学进行对照；

第二类是望远镜级实验读出信息，用于直接构造 $\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$ 中可能同时有 α 和 Be。下面的每层标量总能量可作监视谱，但不能单独保存两粒子的独立测量。用于实际重建时，还要保存每面、每条被触发条带的能量和编号，再做 front/back 配对、跨层关联与 PID。不要从真值 trackID 直接指定重建的粒子配对。

建议包含：

• t0_E_D1
• t0_E_D2
• t0_E_D3
• t0_E_SSD
• t0_E_CsI
• t0_StripX_D1
• t0_StripY_D1
• t0_mult_D1
• t0_mult_D2
• t0_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_hitT0
• flag_hitT1
• flag_hitT2
• flag_triple

这里的三体符合指 α、$^{10}\mathrm{Be}$ 和反冲质子均被重建：前两者在 $T_0$，质子通常在某一个侧向模块。它不是要求 $T_0,T_1,T_2$ 三台望远镜都触发。flag_hitT0 也不足以证明其中有两种已鉴别的粒子；应在通道关联和 PID 后定义 hasAlpha && hasBe10 && hasProton。

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();
}

每个事件开始时将缓存能量和多重性清零、未命中条带编号设为 −1，事件结束后再统一 Fill()。本例的单一 RootIO 写法用于串行运行；多线程运行需要线程独立输出或 Geant4 的分析合并机制，不能让多个 worker 同时写同一个 TTree。

---

14 ROOT 文件的打开方式

程序运行结束后，将得到 ROOT 文件，例如

g4sim.root

打开 ROOT 文件的方法为

root -l g4sim.root

进入 ROOT 之后，首先读取文件和树：

TFile *f = TFile::Open("g4sim.root");
if (!f || f->IsZombie()) throw std::runtime_error("Cannot open g4sim.root");
TTree *tr = f->Get<TTree>("evt");
if (!tr) throw std::runtime_error("Missing evt Tree");
tr->Print();

这样即可查看树中包含的分支。

如果希望直接在 ROOT 命令行中检查事件数，可写为

tr->GetEntries();

---

15 模拟数据的后续分析

从通道读数走到 4.5–4.6 的粒子四动量，还需要中间一步：完成 DSSD 两面配对、相邻条带合并和跨层关联，用 $\Delta E-E$ 做 PID，再从停止层和各层沉积恢复粒子动能。探测器入口的动能还需按路径修正死层及靶中能损，才可与反应点处的运动学比较。

经这些步骤整理出 α、Be、反冲质子的动能和方向后，才能复用 4.5 的 Q 值公式以及 4.6 的分支选择与不变质量重建。对模拟和实验采用相同定义的观测量与选择；4.6 的数值门限是 toy MC 示例，不能不经检查就照搬到带真实几何的模拟。

依次比较反应真值、通过几何的事件、形成读出的事件、通过 PID/Q 门的事件。这样才能区分物理输入、接受度、分辨和选择各自对最终谱的影响，而不是用修改输入峰形来补偿重建问题。

关于用户类和 hit 的完整接口，参见 Geant4 Application Developers Guide。