
简介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(fout_{wl}.log) as f: text f.read() # 日志里匹配 Detected photons: N m re.search(rDetected photons:\s*(\d), text) detected int(m.group(1)) efficiencies.append(detected / 50000.0) plt.plot(wavelengths, efficiencies, markero) plt.xlabel(Wavelength (nm)) plt.ylabel(Collection efficiency) plt.grid(True) plt.savefig(eff_vs_wl.png)逻辑说明正则抓日志里的命中数除以总光子数得到效率。参数说明如果日志格式改了正则跟着改50000要跟宏里的beamOn一致。跑完这条曲线再跟实验数据叠一起就能看出是表面模型偏了还是折射率表不准。验证方法拿一个已知收集效率的简单几何比如裸光纤对接闪烁体跑一遍如果仿真值跟理论值差在 5% 以内说明链路是通的。从那以后我每次改光学参数都强制先跑这个基线几何对一遍确认没引入新偏差再上复杂模型。希望帮到你。本文还有配套的精品资源点击获取