简介：G4OpticksTest 是一份面向高能物理仿真与光学光子传输开发者的 C++ 示例工程，重点演示如何在 Geant4 框架内调用光学元件完成光学光子的创建与传播，并正确配置光学环境。它适合已具备 Geant4 基础、希望将 GPU 加速光线追踪引入探测器模拟的中高级研究人员参考。资源包共 105 个文件，约 179KB，以 27 个 hh 头文件与 27 个 cc 源文件为核心，辅以 11 个 gdml 几何描述、11 个 sh 与 4 个 bash 脚本、9 个 mac 宏命令及少量 c、py、xml 等配置与辅助文件，覆盖探测器构建、敏感探测器定义与批处理运行等模块。工程依赖 NVIDIA OptiX 6.5、CUDA 11.0、Geant4 10.6.p03、Boost、Xerces-C、OpenMesh、GLFW 等外部库，并附带环境变量脚本供读者按本机路径调整。目前已有 599 人学习，可帮助读者理解光学光子传播链路、GPU 加速集成方式与依赖排错思路。
1. 光学光子仿真为什么值得单独拆一个工程
做探测器光学仿真的同行大概都有过这种体验：机械结构、灵敏面积、反射率都调好了，结果光收集效率对不上实验数据，回头一查，问题出在光子追迹那一段——要么是表面模型选错，要么是波长相关的折射率没喂对。Geant4 本身有一套完整的光学光子物理过程，但真要把「产生光子、在介质界面折射反射、被吸收或探测到」这条链路跑通，光靠官方例子往往不够直观。G4OpticksTest 这个示例的价值就在这儿：它把 Geant4 内部的光学元件建模和光子传播流程单独拎出来，做成一个可以编译、可以改参数、可以看每一步结果的最小工程。适合两类人——刚接触光学仿真的新手拿来跑通第一条链路，以及已经在用 Geant4 但光学部分总是「玄学对不上」的老手拿来当调试基线。下面按「它到底怎么组织 → 怎么编译跑起来 → 参数怎么调 → 坑在哪」的顺序拆。
2. 光学光子传播在 Geant4 里是怎么落地的
2.1 光学物理过程与光学表面模型的绑定关系
Geant4 的光学光子不是靠一个「光学物理」包全包圆的，它拆成两层：一层是物理过程，比如G4OpAbsorption、G4OpRayleigh、G4OpBoundaryProcess，负责光子跟物质、界面的相互作用；另一层是光学表面，用G4OpticalSurface描述界面的粗糙度、反射类型、反射率模型。这两层必须显式绑定，否则光子到了界面要么直接穿过去，要么按默认镜面反射走，结果跟真实探测器差得远。
G4OpticksTest 里把这两层的关系写得很清楚：每个光学界面先建G4OpticalSurface，设SetModel、SetType、SetFinish，再通过G4LogicalSkinSurface或G4LogicalBorderSurface挂到对应的 logical volume 上。物理过程则在PhysicsList里注册，光子追迹时由G4OpBoundaryProcess去查当前界面的表面属性。这个顺序不能反，反了表面属性读不到，光子行为就退化成默认值。
常见做法是：先确定界面属于「介质-介质」还是「介质-涂层」，前者用 border surface，后者用 skin surface。G4OpticksTest 的示例结构里两种都给了，方便对照。
2.2 光子产生与追迹的代码骨架
光子源在示例里用G4ParticleGun或G4GeneralParticleSource产生，粒子类型设成opticalphoton，能量按波长换算。下面这段是典型的初级光子源配置，抄的时候注意单位：
// 创建光子源，发射 opticalphoton G4ParticleGun* photonGun = new G4ParticleGun(1); photonGun->SetParticleDefinition( G4OpticalPhoton::OpticalPhotonDefinition()); // 设定波长 450 nm，换算成能量 G4double wavelength = 450. * nm; G4double energy = 1239.84193 * eV * nm / wavelength; photonGun->SetParticleEnergy(energy); // 发射方向沿 +z，位置在 scintillator 中心 photonGun->SetParticlePosition(G4ThreeVector(0., 0., 0.)); photonGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., 1.));逻辑说明：OpticalPhotonDefinition()拿到的是光学光子粒子定义，它跟普通 gamma 不同，不参与电磁级联，只走光学过程。能量换算用1239.84 eV·nm / λ，这个常数记牢，改波长时直接套。参数说明：SetParticleEnergy必须跟材料的光学属性表里的能量点对得上，否则查折射率时会插值到边界外，Geant4 会警告甚至返回默认值。
追迹部分由G4OpBoundaryProcess在 stepping 里自动处理，但要在SteppingAction里记录光子到达探测面的状态，示例里用G4OpBoundaryProcess::GetStatus()判断是折射、反射还是吸收。这一步是后面统计收集效率的关键。
2.3 材料光学属性表的填法
光学仿真翻车最多的地方就是材料属性表。Geant4 要求用G4MaterialPropertiesTable给每个材料挂折射率、吸收长度、反射率等，而且能量点必须单调递增。示例里给闪烁体和光导分别建了表：
// 闪烁体折射率，按能量点给 G4MaterialPropertiesTable* mpt = new G4MaterialPropertiesTable(); const G4int nEntries = 2; G4double photonEnergy[nEntries] = {2.0 * eV, 3.5 * eV}; G4double rIndex[nEntries] = {1.58, 1.60}; mpt->AddProperty("RINDEX", photonEnergy, rIndex, nEntries); material->SetMaterialPropertiesTable(mpt);逻辑说明：RINDEX是 Geant4 约定的关键字，不能拼错。能量点从低到高排，示例里只给两个点做线性插值，实际项目建议按实测数据给 5 到 10 个点。参数说明：如果材料同时要参与光学和普通物理，属性表可以共用，但RINDEX只影响光学过程，不会干扰电磁过程。
提示：属性表里没给
ABSLENGTH的材料，Geant4 默认光子不吸收，长距离追迹时会出现光子「永远不死」的情况，统计结果偏乐观。
3. 把示例编译跑起来：从 CMake 到可视化
3.1 编译环境与依赖检查
G4OpticksTest 依赖 Geant4 本体，编译前先确认geant4-config --version能正常输出，并且环境变量G4INSTALL、G4LEDATA这些已经 source 过。常见做法是用 CMake 组织，示例里给的CMakeLists.txt核心就三行：
find_package(Geant4 REQUIRED) include(${Geant4_USE_FILE}) add_executable(G4OpticksTest main.cc src/*.cc) target_link_libraries(G4OpticksTest ${Geant4_LIBRARIES})逻辑说明：find_package会去Geant4_DIR找配置，如果报找不到，多半是没 sourcegeant4.sh。参数说明：Geant4_USE_FILE会自动带上 include 路径和编译选项，不用手写-I。编译命令：
mkdir build && cd build cmake .. -DGeant4_DIR=$(geant4-config --prefix)/lib/Geant4-<版本> make -j$(nproc)把<版本>换成实际版本号。跑起来后如果提示找不到G4OpticksTest可执行文件，检查add_executable里的源文件路径是否跟实际目录一致。
3.2 运行宏文件与可视化配置
示例通常带一个run.mac，里面控制粒子数、发射位置、是否开可视化。最小可跑宏：
/run/initialize /gun/particle opticalphoton /gun/energy 2.75 eV /run/beamOn 1000逻辑说明：/run/initialize必须在设置粒子之前，否则物理表没建好。/run/beamOn 1000表示发 1000 个光子，统计收集效率时这个数至少上万才有统计意义。参数说明：如果开了可视化，/vis/open OGL之后加/vis/scene/add/trajectories能看到光子轨迹，但轨迹多了会卡，调试阶段建议先关可视化跑批。
3.3 输出结果与收集效率统计
示例里一般在SteppingAction或EventAction里记录到达探测面的光子数。常见做法是在探测体里设一个G4VSensitiveDetector，命中时计数。跑完 1000 个光子后，终端会打印类似Detected photons: 342的行，收集效率就是 342/1000。如果这个数明显偏离预期，先查表面模型和折射率表，再查光子源位置是否在灵敏体积内。
4. 参数怎么调：波长、表面模型与统计
4.1 波长扫描与折射率插值
光学仿真里波长一变，折射率和吸收长度都跟着变，所以单点跑没意义，要做扫描。示例里可以用宏循环或者改ParticleGun能量，常见做法是写个 shell 循环：
for wl in 400 450 500 550 600; do energy=$(echo "1239.84193 / $wl" | bc -l) sed -i "s|/gun/energy .*|/gun/energy ${energy} eV|" run.mac ./G4OpticksTest run.mac > out_${wl}.log done逻辑说明：bc -l做浮点除法，sed替换宏里的能量行。参数说明：波长单位是 nm，能量单位 eV，换算常数 1239.84 别记错。跑完对比不同波长的收集效率，能看出材料色散对结果的影响。
4.2 表面模型选型：polished、ground、unified 的差别
G4OpticalSurface的SetModel有三个常用值：glisur、unified、LUT。glisur是老模型，只支持 polished 和 ground 两种 finish，适合快速验证；unified支持更细的粗糙度参数，比如SigmaAlpha，适合模拟真实抛光面；LUT查表法，适合有实测 BRDF 数据的场景。示例里默认给的是unified，因为它在粗糙度和反射率之间平衡得比较好。
选型理由：如果只是验证几何光路，glisur够用且快；如果要跟实验的角分布对比，必须上unified，并且把SigmaAlpha调到跟实测粗糙度一致。常见误用是拿glisur的 ground 去模拟喷砂面，结果反射方向偏得离谱。
4.3 统计误差与光子数下限
光子追迹是蒙特卡洛，统计误差按1/sqrt(N)降。跑 1000 个光子，相对误差大概 3%；要降到 1% 以下，至少 10000 个。示例里默认beamOn 1000只是演示，真做效率对比时我一般直接上 50000。参数说明：光子数上去后运行时间线性涨，可以用多线程G4MTRunManager加速，但要注意随机数种子别重复。
5. 避坑与排查：光学仿真里最容易翻车的五件事
5.1 现象：光子全部直接穿过界面，没有反射
原因：G4LogicalBorderSurface或G4LogicalSkinSurface没挂上，或者挂的时候 volume 名字写错。解决：在DetectorConstruction里打印每个 surface 挂载的 volume 名，确认跟G4LogicalVolume的实例名一致。常见做法是挂完后用G4LogicalBorderSurface::GetSurface()反查一遍。
5.2 现象：运行时报「RINDEX not found」
原因：材料属性表里关键字拼错，或者能量点没按升序排。解决：检查AddProperty的第一个参数是否全大写，能量数组是否单调递增。如果材料是从 GDML 读的，确认 GDML 里matrix的energy顺序。
5.3 现象：收集效率随光子数剧烈波动
原因：光子数太少，统计涨落大；或者光子源位置在边界上，部分光子初始就在灵敏体积外。解决：把beamOn提到 10000 以上，光子源位置往体积中心挪，并检查SetParticlePosition的单位是不是 mm。
5.4 现象：开了可视化后程序卡死
原因：轨迹点太多，OpenGL 渲染跟不上。解决：调试阶段用/vis/scene/add/trajectories时加/vis/scene/endOfEventAction accumulate 10限制累积事件数，或者干脆关可视化跑批，用日志分析。
5.5 现象：多线程下结果跟单线程对不上
原因：随机数种子没按线程隔离，或者共享的计数器没加锁。解决：用G4MTRunManager时给每个线程设独立种子，计数器用G4Accumulable或原子操作。常见做法是先在单线程下调通，再开多线程对比前 100 个事件的结果。
6. 进阶：用光学光子做波长分辨的收集效率曲线
把前面几章串起来，最终要拿到的是一条「收集效率 vs 波长」曲线。做法是：固定几何和表面模型，扫 400 到 700 nm，每个波长跑 50000 个光子，统计探测面命中数。下面这段是后处理脚本的骨架，读日志、算效率、画曲线：
import re import matplotlib.pyplot as plt wavelengths = [400, 450, 500, 550, 600, 650, 700] efficiencies = [] for wl in wavelengths: with open(f"out_{wl}.log") as f: text = f.read() # 日志里匹配 "Detected photons: N" m = re.search(r"Detected photons:\s*(\d+)", text) detected = int(m.group(1)) efficiencies.append(detected / 50000.0) plt.plot(wavelengths, efficiencies, marker='o') plt.xlabel("Wavelength (nm)") plt.ylabel("Collection efficiency") plt.grid(True) plt.savefig("eff_vs_wl.png")逻辑说明：正则抓日志里的命中数，除以总光子数得到效率。参数说明：如果日志格式改了，正则跟着改；50000要跟宏里的beamOn一致。跑完这条曲线，再跟实验数据叠一起，就能看出是表面模型偏了还是折射率表不准。
验证方法：拿一个已知收集效率的简单几何（比如裸光纤对接闪烁体）跑一遍，如果仿真值跟理论值差在 5% 以内，说明链路是通的。从那以后我每次改光学参数，都强制先跑这个基线几何对一遍，确认没引入新偏差再上复杂模型。希望帮到你。
本文还有配套的精品资源，点击获取
