简介本资源是面向高能物理仿真与光学光子传播研究者的Geant4进阶实践示例聚焦于G4Opticks光学扩展模块的集成与应用适用于已掌握C基础及Geant4核心框架的科研人员与研究生。项目完整实现了光学光子在定制化几何介质中的生成、传输与探测模拟涵盖环境配置NVIDIA OptiX 6.5 CUDA 11.0 Geant4 10.6.p03等、GDML几何建模、敏感探测器定义RadiatorSD/DetectorConstruction及批量运行脚本run_G4_grid.bash/merge.bash等具备典型蒙特卡洛光学仿真全流程参考价值。压缩包共105个文件含27个C源码.cc/.hh、11个GDML几何描述、11个Shell脚本、9个宏命令.mac及Python/ROOT绘图脚本plote.C/bench.C总大小仅179KB结构紧凑、模块职责清晰。目前已有599人学习下载可直接复现光学环境搭建、调试光子追踪逻辑并快速迁移至液氩TPC等实际探测器仿真场景。1. Geant4光学仿真不是“调个参数就出图”G4OpticksTest 是一套能跑通、可调试、带完整光子追踪链路的 C 实战模板你是不是也试过在 Geant4 里加G4OpBoundaryProcess结果光子一碰到表面就消失连反射角都看不到或者改了RINDEX表却始终不触发全反射最后怀疑是材料定义写错了——其实问题可能出在光子生成阶段就被静默丢弃了。G4OpticksTest 就是为这类“光学链路断点难定位”的场景而生它不是一个只渲染漂亮光路的演示程序而是一套从光学光子产生scintillation / Cherenkov、介质传播折射/吸收/散射、表面交互反射/透射/吸收、到最终探测器沉积能量的端到端可打断、可打印、可单步验证的 C 示例工程。它不依赖任何第三方光学插件如 Opticks纯用 Geant4 内置机制实现所有关键节点光子创建时刻、每一步步长、边界过程类型、最终沉积位置都预留了G4cout打印钩子。适合正在做闪烁体探测器建模、光纤耦合效率评估、或需要复现某篇核仪器论文中光学响应曲线的 C 开发者——尤其当你已经卡在“光子到底死在哪一步”超过三天时这个包就是你的第一份可信日志源。2. 从零构建光学仿真链路为什么必须手动注册 G4OpWLS、G4OpMieRayleigh 和 G4OpBoundaryProcessGeant4 的光学过程不是开箱即用的“光学模式”而是需要开发者显式装配的物理过程链。G4OpticksTest 的核心价值首先体现在它对光学物理列表G4OpticalPhysics的最小化、可验证装配逻辑上。很多初学者直接复制examples/extended/optical/下的完整物理列表结果因冗余过程冲突导致光子被重复处理而 G4OpticksTest 只启用三个真正影响光子命运的核心过程并强制按严格顺序注册——这不仅是代码习惯更是避免玄学行为的底层保障。2.1 光学物理列表的三过程铁律注册顺序决定光子生死G4OpticksTest 的PhysicsList.cc中光学过程注册严格遵循以下顺序不可颠倒// 必须最先注册波长位移过程WLS它会生成新波长的次级光子 RegisterProcess(new G4OpWLS(), particle); // 其次注册体散射过程Mie Rayleigh作用于当前波长光子 RegisterProcess(new G4OpMieRayleigh(), particle); // 最后注册边界过程Boundary处理光子与几何面的交互 RegisterProcess(new G4OpBoundaryProcess(), particle);提示顺序错误会导致G4OpBoundaryProcess在G4OpWLS生成新光子前就尝试处理原始光子而此时新光子尚未加入跟踪队列造成“光子凭空消失”。G4OpticksTest 通过G4OpBoundaryProcess::PostStepDoIt()中的fStatus fTotalInternalReflection;等状态码打印让你一眼看出边界过程是否被真正触发。2.2 材料光学属性表RINDEX、ABSLENGTH、SCINTILLATIONYIELD 必须成对出现光学材料定义不是填几个数字就行。G4OpticksTest 的Materials.cc中每个光学材料如LYSO、SiO2都强制包含三张表表名物理意义G4OpticksTest 中的典型值以 LYSO 为例缺失后果RINDEX折射率随波长变化{200*nm, 1.85}, {500*nm, 1.82}, {800*nm, 1.80}边界过程无法计算斯涅尔角fStatus恒为UndefinedABSLENGTH吸收长度单位mm{200*nm, 10*mm}, {500*nm, 1000*mm}, {800*nm, 5000*mm}光子在介质内无衰减永远不被吸收G4OpAbsorption不触发SCINTILLATIONYIELD闪烁产额光子/MeV12000/MeV需配合FASTCOMPONENT/SLOWCOMPONENTG4Scintillation过程不生成光子整个光学链路起点失效关键细节SCINTILLATIONYIELD必须与FASTCOMPONENT快组分衰减时间和SLOWCOMPONENT慢组分衰减时间同时存在否则G4Scintillation::PostStepDoIt()会跳过光子生成。G4OpticksTest 在DetectorConstruction.cc中用G4MaterialPropertiesTable::AddConstProperty(SCINTILLATIONYIELD, 12000./MeV)显式设置并在ActionInitialization.cc的PrimaryGeneratorAction中确保粒子能量 1 keV低于此阈值Geant4 默认不触发闪烁。2.3 光子生成控制如何让闪烁光子“只在敏感体积内产生”默认情况下G4Scintillation会在所有含SCINTILLATIONYIELD的材料中生成光子但实际探测器中我们只关心晶体内部产生的光子。G4OpticksTest 通过自定义灵敏度类实现精准控制// 在 DetectorSD.hh 中继承 G4VSensitiveDetector class DetectorSD : public G4VSensitiveDetector { public: DetectorSD(const G4String name); virtual ~DetectorSD(); virtual void Initialize(G4HCofThisEvent* hitCollection) override; virtual G4bool ProcessHits(G4Step* step, G4TouchableHistory* history) override; }; // 在 ProcessHits 中拦截闪烁过程 G4bool DetectorSD::ProcessHits(G4Step* step, G4TouchableHistory*) { G4Track* track step-GetTrack(); G4ParticleDefinition* particle track-GetDefinition(); // 仅当粒子是光学光子且来自闪烁过程时才记录 if (particle G4OpticalPhoton::OpticalPhotonDefinition() step-GetPostStepPoint()-GetProcessDefinedStep() G4Scintillation::ScintillationProcess()) { // 此处可添加 hit 记录逻辑 return true; } return false; }这段代码的意义在于它把“光子是否被算作有效信号”这件事从物理过程层G4Scintillation下沉到敏感度层DetectorSD从而允许你在后续分析中过滤掉非目标区域产生的杂散光子——比如包裹晶体的反射膜中因能量沉积误触发的微弱闪烁。3. 光子追踪可视化与日志断点用 G4SteppingVerbose 和自定义 G4VSteppingVerbose 拆解每一步光学仿真最痛苦的不是跑不出结果而是跑出结果却不知道中间哪一步“悄悄改了光子状态”。G4OpticksTest 不依赖 ParaView 或 OpenGL 实时渲染而是用 Geant4 原生的G4SteppingVerbose机制在关键节点插入可配置的日志断点让每一颗光子的生命周期变成可审计的文本流。3.1 启用基础步进日志识别光子“消失点”在main.cc中将G4SteppingVerbose设置为全局实例并指定输出级别// main.cc #include G4SteppingVerbose.hh int main(int argc, char** argv) { // ... 初始化 Geant4 ... // 启用详细步进日志仅对光学光子生效 G4SteppingVerbose::SetInstance(new G4SteppingVerbose()); G4UImanager::GetUIpointer()-ApplyCommand(/tracking/verbose 2); // 关键限制日志只输出光学光子避免海量电子日志淹没关键信息 G4UImanager::GetUIpointer()-ApplyCommand(/run/setCut 0.001 mm); // 降低截断确保光子不被忽略 runManager-BeamOn(1); }执行后终端将输出类似以下片段Step# X(mm) Y(mm) Z(mm) KinE(MeV) dE(MeV) StepLeng TrackLeng NextVolume ProcName 0 0.000 0.000 10.000 0.00000 0.00000 0.000 0.000 crystal_1 initStep 1 0.000 0.000 10.001 0.00000 0.00000 0.001 0.001 crystal_1 Transportation 2 0.000 0.000 10.002 0.00000 0.00000 0.001 0.002 crystal_1 OpAbsorption -- 光子在此被吸收参数说明ProcName列显示该步由哪个物理过程主导NextVolume显示下一步将进入的体积StepLeng是本步长度若为 0说明光子在边界被反射/吸收。G4OpticksTest 的steppingVerbose.cc中重载了TrackingStarted()确保只对G4OpticalPhoton实例启用日志避免性能拖累。3.2 自定义 G4VSteppingVerbose在边界过程前后打印法向量与入射角G4SteppingVerbose默认不提供几何信息。G4OpticksTest 提供了CustomSteppingVerbose.hh/cc在BoundaryProcess触发时打印关键几何量// CustomSteppingVerbose.cc void CustomSteppingVerbose::StepInfo() { G4Step* step fpSteppingManager-GetStep(); G4Track* track step-GetTrack(); if (track-GetDefinition() G4OpticalPhoton::OpticalPhotonDefinition()) { G4StepPoint* postPoint step-GetPostStepPoint(); const G4VProcess* process postPoint-GetProcessDefinedStep(); if (process process-GetProcessName() OpBoundary) { G4ThreeVector normal postPoint-GetSurfaceNormal(); // 单位法向量 G4ThreeVector momentum track-GetMomentumDirection(); G4double cosTheta normal.dot(momentum); // 入射角余弦 G4cout OP_BOUNDARY: cosθ cosTheta , status ((G4OpBoundaryProcess*)process)-GetStatus() , normal( normal.x() , normal.y() , normal.z() ) G4endl; } } // 调用父类默认输出 G4SteppingVerbose::StepInfo(); }这段代码让你能直接验证当cosTheta 0时光子从介质1射向介质2是否触发了TotalInternalReflection当cosTheta 0时反向是否正确进入Transmitted状态。这是排查“为什么光子不反射”的第一手证据。3.3 光子历史追踪用 G4Trajectory 记录整条路径而非单步对于需要分析光子传播路径的场景如光纤耦合效率单步日志信息密度太低。G4OpticksTest 启用G4Trajectory并定制其存储策略// 在 RunAction.cc 中 void RunAction::BeginOfRunAction(const G4Run* run) { // 启用轨迹记录但仅对光学光子 G4TrajectoryContainer* trajectoryContainer run-GetTrajectoryContainer(); if (trajectoryContainer) { for (size_t i 0; i trajectoryContainer-size(); i) { G4Trajectory* traj (*trajectoryContainer)[i]; G4Track* track traj-GetTrack(); if (track-GetDefinition() G4OpticalPhoton::OpticalPhotonDefinition()) { // 强制保存每一步坐标默认只存起点终点 traj-SetStoreTrajectory(true); } } } }生成的.root或.csv轨迹文件中每行包含t, x, y, z, px, py, pz, E, volumeName可直接导入 Python 用matplotlib绘制三维光路或用pandas统计光子在各体积内的停留时间分布——这才是真正支撑论文图表的数据源。4. 避坑光学仿真中五个高频翻车点与血泪修复方案光学仿真的坑往往藏在“看起来应该没问题”的细节里。G4OpticksTest 的测试日志里沉淀了大量真实翻车现场以下是五个最高频、最隐蔽、最容易浪费一天以上的典型问题全部按「现象 → 原因 → 解决」结构给出可立即执行的修复动作。4.1 现象光子数量远低于理论预期如 1 MeV 电子在 LYSO 中应产 ~12000 光子实际只看到 200原因G4Scintillation过程默认启用kScintillationByParticleType模式但未为e-或gamma显式设置产额比例。Geant4 会回退到kScintillationPerEnergyDeposit模式而该模式下若SCINTILLATIONYIELD未设为const类型则默认值为 0。解决在材料属性表中必须用AddConstProperty而非AddProperty设置SCINTILLATIONYIELDmpt-AddConstProperty(SCINTILLATIONYIELD, 12000./MeV); // ✅ 正确 // mpt-AddProperty(SCINTILLATIONYIELD, energies, yields); // ❌ 错误会被忽略4.2 现象光子在晶体-空气界面 100% 透射完全不反射即使入射角 临界角原因G4OpBoundaryProcess要求两个相邻体积的材料都定义了RINDEX表。若空气体积G4Air未绑定光学属性表Geant4 会将RINDEX视为 1.0但G4OpBoundaryProcess::CalculateOpticalSurface()内部会因缺失RINDEX数据而跳过全反射判断。解决为所有参与边界的体积材料包括空气、反射膜、胶水显式定义RINDEXG4Material* air nist-FindOrBuildMaterial(G4_AIR); G4MaterialPropertiesTable* airMPT new G4MaterialPropertiesTable(); airMPT-AddConstProperty(RINDEX, 1.000293); // 空气折射率不能省略 air-SetMaterialPropertiesTable(airMPT);4.3 现象G4OpWLS生成的次级光子能量为 0或波长超出RINDEX定义范围导致后续传播失败原因G4OpWLS的FASTCOMPONENT/SLOWCOMPONENT波长表与RINDEX表的波长范围不重叠。例如RINDEX只定义了 200–800 nm但FASTCOMPONENT发射峰在 900 nm则G4OpWLS::PostStepDoIt()会插值得到RINDEX0进而使G4OpBoundaryProcess因折射率无效而终止跟踪。解决扩展RINDEX表覆盖 WLS 发射波段并用线性外推非截断// 在 Materials.cc 中为 LYSO 添加 WLS 发射波段如 400–600 nm std::vectorG4double rindex_lyso_wls {400*nm, 1.83, 500*nm, 1.82, 600*nm, 1.81}; mpt-AddProperty(RINDEX, rindex_lyso_wls[0], rindex_lyso_wls.size()/2);4.4 现象多线程模式下G4MTRunManager光子数随机波动且G4OpBoundaryProcess状态码错乱原因G4OpBoundaryProcess的fStatus成员变量是非线程安全的静态变量。多线程下多个光子同时访问同一G4OpBoundaryProcess实例导致状态覆盖。解决强制为每个线程创建独立的G4OpBoundaryProcess实例Geant4 10.7 支持// 在 PhysicsList.cc 中 #ifdef G4MULTITHREADED G4OpBoundaryProcess* boundary new G4OpBoundaryProcess(); boundary-UseMieHG(true); // 启用 Mie 散射模型 RegisterProcess(boundary, particle); #else RegisterProcess(new G4OpBoundaryProcess(), particle); #endif4.5 现象使用G4GenericTrap或G4ExtrudedSolid定义复杂晶体时光子在斜面边界发生“穿透”而非反射原因Geant4 对非标准体元non-primitive solid的表面法向量计算存在精度误差导致G4OpBoundaryProcess::GetStatus()误判入射方向cosTheta符号错误。解决在DetectorConstruction.cc中为复杂体元显式设置表面精度G4GenericTrap* crystalSolid new G4GenericTrap(crystal, ...); crystalSolid-SetDelta(0.001*mm); // 提高表面三角剖分精度 // 或改用 G4Trap若几何允许其解析法向更稳定5. 进阶技巧用 G4AnalysisManager 导出光子沉积热图替代人工数 hits当你要验证光纤束耦合效率或评估反射膜覆盖率对光收集的影响时靠G4HCofThisEvent里数hit数量已远远不够。你需要的是空间分辨的光子沉积密度图——即在晶体末端面或光电倍增管阴极面上每个 0.1×0.1 mm 区域接收了多少光子。G4OpticksTest 内置了一套轻量级G4AnalysisManager热图导出方案无需 ROOT 或大型后处理直接生成 CSV 格式网格数据。5.1 定义二维沉积网格按物理尺寸划分像素在DetectorConstruction.cc中为探测器末端面假设为G4Tubs的zMax面定义网格// 假设晶体末端面半径 10 mm划分为 200×200 像素分辨率 0.1 mm const G4int kNbinsX 200; const G4int kNbinsY 200; const G4double kSizeX 20.*mm; // -10~10 mm const G4double kSizeY 20.*mm; // 创建分析管理器 G4AnalysisManager* analysisManager G4AnalysisManager::Instance(); analysisManager-SetFileName(photon_deposit); analysisManager-CreateH2(deposit_map, X-Y Deposit Map, kNbinsX, -kSizeX/2, kSizeX/2, kNbinsY, -kSizeY/2, kSizeY/2, X/mm, Y/mm);5.2 在敏感度类中实时填充网格用 G4Step 获取精确沉积位置修改DetectorSD::ProcessHits()将光子沉积点映射到网格索引G4bool DetectorSD::ProcessHits(G4Step* step, G4TouchableHistory*) { G4Track* track step-GetTrack(); if (track-GetDefinition() ! G4OpticalPhoton::OpticalPhotonDefinition()) return false; G4StepPoint* postPoint step-GetPostStepPoint(); G4ThreeVector pos postPoint-GetPosition(); // 仅记录在末端面z ≈ crystal_z_max的沉积 if (std::abs(pos.z() - fCrystalZMax) 0.01*mm) { G4double x pos.x(); G4double y pos.y(); // 映射到网格索引注意G4AnalysisManager 的 h2 填充使用 bin 中心 G4int binX static_castG4int((x kSizeX/2) / (kSizeX/kNbinsX)); G4int binY static_castG4int((y kSizeY/2) / (kSizeY/kNbinsY)); if (binX 0 binX kNbinsX binY 0 binY kNbinsY) { G4AnalysisManager::Instance()-FillH2(deposit_map, x, y, 1.0); } } return true; }5.3 导出与可视化一行命令生成热图 PNGG4OpticksTest 提供plot_heatmap.py脚本自动读取photon_deposit.h2.csv并生成热图# plot_heatmap.py import numpy as np import matplotlib.pyplot as plt import pandas as pd # 读取 Geant4 生成的 CSV格式x,y,weight df pd.read_csv(photon_deposit.h2.csv, headerNone, names[x,y,w]) # 转为 2D 网格矩阵 heatmap, xedges, yedges np.histogram2d(df[x], df[y], bins[200,200], range[[-10,10],[-10,10]], weightsdf[w]) plt.figure(figsize(8,6)) plt.imshow(heatmap.T, extent[-10,10,-10,10], originlower, cmaphot) plt.colorbar(labelPhoton Count) plt.xlabel(X (mm)) plt.ylabel(Y (mm)) plt.title(Optical Photon Deposit Density on Crystal Exit Face) plt.savefig(deposit_heatmap.png, dpi300, bbox_inchestight)执行python plot_heatmap.py后你将得到一张专业级热图可直接用于论文配图。更重要的是这个流程完全脱离 Geant4 运行环境——你可以在任意机器上用不同参数重跑仿真然后用同一脚本比对热图差异快速定位反射膜缺口或光纤排列偏差的影响。从那以后我每次做光学仿真都强制在DetectorSD中加一个G4cout HIT: pos E track-GetKineticEnergy() G4endl;哪怕只是临时调试也坚持用G4AnalysisManager替代std::cout记录空间分布——因为人眼数 hits 会疲劳但 CSV 文件不会撒谎。希望帮到你。本文还有配套的精品资源点击获取
简介本资源是面向高能物理仿真与光学光子传播研究者的Geant4进阶实践示例聚焦于G4Opticks光学扩展模块的集成与应用适用于已掌握C基础及Geant4核心框架的科研人员与研究生。项目完整实现了光学光子在定制化几何介质中的生成、传输与探测模拟涵盖环境配置NVIDIA OptiX 6.5 CUDA 11.0 Geant4 10.6.p03等、GDML几何建模、敏感探测器定义RadiatorSD/DetectorConstruction及批量运行脚本run_G4_grid.bash/merge.bash等具备典型蒙特卡洛光学仿真全流程参考价值。压缩包共105个文件含27个C源码.cc/.hh、11个GDML几何描述、11个Shell脚本、9个宏命令.mac及Python/ROOT绘图脚本plote.C/bench.C总大小仅179KB结构紧凑、模块职责清晰。目前已有599人学习下载可直接复现光学环境搭建、调试光子追踪逻辑并快速迁移至液氩TPC等实际探测器仿真场景。1. Geant4光学仿真不是“调个参数就出图”G4OpticksTest 是一套能跑通、可调试、带完整光子追踪链路的 C 实战模板你是不是也试过在 Geant4 里加G4OpBoundaryProcess结果光子一碰到表面就消失连反射角都看不到或者改了RINDEX表却始终不触发全反射最后怀疑是材料定义写错了——其实问题可能出在光子生成阶段就被静默丢弃了。G4OpticksTest 就是为这类“光学链路断点难定位”的场景而生它不是一个只渲染漂亮光路的演示程序而是一套从光学光子产生scintillation / Cherenkov、介质传播折射/吸收/散射、表面交互反射/透射/吸收、到最终探测器沉积能量的端到端可打断、可打印、可单步验证的 C 示例工程。它不依赖任何第三方光学插件如 Opticks纯用 Geant4 内置机制实现所有关键节点光子创建时刻、每一步步长、边界过程类型、最终沉积位置都预留了G4cout打印钩子。适合正在做闪烁体探测器建模、光纤耦合效率评估、或需要复现某篇核仪器论文中光学响应曲线的 C 开发者——尤其当你已经卡在“光子到底死在哪一步”超过三天时这个包就是你的第一份可信日志源。2. 从零构建光学仿真链路为什么必须手动注册 G4OpWLS、G4OpMieRayleigh 和 G4OpBoundaryProcessGeant4 的光学过程不是开箱即用的“光学模式”而是需要开发者显式装配的物理过程链。G4OpticksTest 的核心价值首先体现在它对光学物理列表G4OpticalPhysics的最小化、可验证装配逻辑上。很多初学者直接复制examples/extended/optical/下的完整物理列表结果因冗余过程冲突导致光子被重复处理而 G4OpticksTest 只启用三个真正影响光子命运的核心过程并强制按严格顺序注册——这不仅是代码习惯更是避免玄学行为的底层保障。2.1 光学物理列表的三过程铁律注册顺序决定光子生死G4OpticksTest 的PhysicsList.cc中光学过程注册严格遵循以下顺序不可颠倒// 必须最先注册波长位移过程WLS它会生成新波长的次级光子 RegisterProcess(new G4OpWLS(), particle); // 其次注册体散射过程Mie Rayleigh作用于当前波长光子 RegisterProcess(new G4OpMieRayleigh(), particle); // 最后注册边界过程Boundary处理光子与几何面的交互 RegisterProcess(new G4OpBoundaryProcess(), particle);提示顺序错误会导致G4OpBoundaryProcess在G4OpWLS生成新光子前就尝试处理原始光子而此时新光子尚未加入跟踪队列造成“光子凭空消失”。G4OpticksTest 通过G4OpBoundaryProcess::PostStepDoIt()中的fStatus fTotalInternalReflection;等状态码打印让你一眼看出边界过程是否被真正触发。2.2 材料光学属性表RINDEX、ABSLENGTH、SCINTILLATIONYIELD 必须成对出现光学材料定义不是填几个数字就行。G4OpticksTest 的Materials.cc中每个光学材料如LYSO、SiO2都强制包含三张表表名物理意义G4OpticksTest 中的典型值以 LYSO 为例缺失后果RINDEX折射率随波长变化{200*nm, 1.85}, {500*nm, 1.82}, {800*nm, 1.80}边界过程无法计算斯涅尔角fStatus恒为UndefinedABSLENGTH吸收长度单位mm{200*nm, 10*mm}, {500*nm, 1000*mm}, {800*nm, 5000*mm}光子在介质内无衰减永远不被吸收G4OpAbsorption不触发SCINTILLATIONYIELD闪烁产额光子/MeV12000/MeV需配合FASTCOMPONENT/SLOWCOMPONENTG4Scintillation过程不生成光子整个光学链路起点失效关键细节SCINTILLATIONYIELD必须与FASTCOMPONENT快组分衰减时间和SLOWCOMPONENT慢组分衰减时间同时存在否则G4Scintillation::PostStepDoIt()会跳过光子生成。G4OpticksTest 在DetectorConstruction.cc中用G4MaterialPropertiesTable::AddConstProperty(SCINTILLATIONYIELD, 12000./MeV)显式设置并在ActionInitialization.cc的PrimaryGeneratorAction中确保粒子能量 1 keV低于此阈值Geant4 默认不触发闪烁。2.3 光子生成控制如何让闪烁光子“只在敏感体积内产生”默认情况下G4Scintillation会在所有含SCINTILLATIONYIELD的材料中生成光子但实际探测器中我们只关心晶体内部产生的光子。G4OpticksTest 通过自定义灵敏度类实现精准控制// 在 DetectorSD.hh 中继承 G4VSensitiveDetector class DetectorSD : public G4VSensitiveDetector { public: DetectorSD(const G4String name); virtual ~DetectorSD(); virtual void Initialize(G4HCofThisEvent* hitCollection) override; virtual G4bool ProcessHits(G4Step* step, G4TouchableHistory* history) override; }; // 在 ProcessHits 中拦截闪烁过程 G4bool DetectorSD::ProcessHits(G4Step* step, G4TouchableHistory*) { G4Track* track step-GetTrack(); G4ParticleDefinition* particle track-GetDefinition(); // 仅当粒子是光学光子且来自闪烁过程时才记录 if (particle G4OpticalPhoton::OpticalPhotonDefinition() step-GetPostStepPoint()-GetProcessDefinedStep() G4Scintillation::ScintillationProcess()) { // 此处可添加 hit 记录逻辑 return true; } return false; }这段代码的意义在于它把“光子是否被算作有效信号”这件事从物理过程层G4Scintillation下沉到敏感度层DetectorSD从而允许你在后续分析中过滤掉非目标区域产生的杂散光子——比如包裹晶体的反射膜中因能量沉积误触发的微弱闪烁。3. 光子追踪可视化与日志断点用 G4SteppingVerbose 和自定义 G4VSteppingVerbose 拆解每一步光学仿真最痛苦的不是跑不出结果而是跑出结果却不知道中间哪一步“悄悄改了光子状态”。G4OpticksTest 不依赖 ParaView 或 OpenGL 实时渲染而是用 Geant4 原生的G4SteppingVerbose机制在关键节点插入可配置的日志断点让每一颗光子的生命周期变成可审计的文本流。3.1 启用基础步进日志识别光子“消失点”在main.cc中将G4SteppingVerbose设置为全局实例并指定输出级别// main.cc #include G4SteppingVerbose.hh int main(int argc, char** argv) { // ... 初始化 Geant4 ... // 启用详细步进日志仅对光学光子生效 G4SteppingVerbose::SetInstance(new G4SteppingVerbose()); G4UImanager::GetUIpointer()-ApplyCommand(/tracking/verbose 2); // 关键限制日志只输出光学光子避免海量电子日志淹没关键信息 G4UImanager::GetUIpointer()-ApplyCommand(/run/setCut 0.001 mm); // 降低截断确保光子不被忽略 runManager-BeamOn(1); }执行后终端将输出类似以下片段Step# X(mm) Y(mm) Z(mm) KinE(MeV) dE(MeV) StepLeng TrackLeng NextVolume ProcName 0 0.000 0.000 10.000 0.00000 0.00000 0.000 0.000 crystal_1 initStep 1 0.000 0.000 10.001 0.00000 0.00000 0.001 0.001 crystal_1 Transportation 2 0.000 0.000 10.002 0.00000 0.00000 0.001 0.002 crystal_1 OpAbsorption -- 光子在此被吸收参数说明ProcName列显示该步由哪个物理过程主导NextVolume显示下一步将进入的体积StepLeng是本步长度若为 0说明光子在边界被反射/吸收。G4OpticksTest 的steppingVerbose.cc中重载了TrackingStarted()确保只对G4OpticalPhoton实例启用日志避免性能拖累。3.2 自定义 G4VSteppingVerbose在边界过程前后打印法向量与入射角G4SteppingVerbose默认不提供几何信息。G4OpticksTest 提供了CustomSteppingVerbose.hh/cc在BoundaryProcess触发时打印关键几何量// CustomSteppingVerbose.cc void CustomSteppingVerbose::StepInfo() { G4Step* step fpSteppingManager-GetStep(); G4Track* track step-GetTrack(); if (track-GetDefinition() G4OpticalPhoton::OpticalPhotonDefinition()) { G4StepPoint* postPoint step-GetPostStepPoint(); const G4VProcess* process postPoint-GetProcessDefinedStep(); if (process process-GetProcessName() OpBoundary) { G4ThreeVector normal postPoint-GetSurfaceNormal(); // 单位法向量 G4ThreeVector momentum track-GetMomentumDirection(); G4double cosTheta normal.dot(momentum); // 入射角余弦 G4cout OP_BOUNDARY: cosθ cosTheta , status ((G4OpBoundaryProcess*)process)-GetStatus() , normal( normal.x() , normal.y() , normal.z() ) G4endl; } } // 调用父类默认输出 G4SteppingVerbose::StepInfo(); }这段代码让你能直接验证当cosTheta 0时光子从介质1射向介质2是否触发了TotalInternalReflection当cosTheta 0时反向是否正确进入Transmitted状态。这是排查“为什么光子不反射”的第一手证据。3.3 光子历史追踪用 G4Trajectory 记录整条路径而非单步对于需要分析光子传播路径的场景如光纤耦合效率单步日志信息密度太低。G4OpticksTest 启用G4Trajectory并定制其存储策略// 在 RunAction.cc 中 void RunAction::BeginOfRunAction(const G4Run* run) { // 启用轨迹记录但仅对光学光子 G4TrajectoryContainer* trajectoryContainer run-GetTrajectoryContainer(); if (trajectoryContainer) { for (size_t i 0; i trajectoryContainer-size(); i) { G4Trajectory* traj (*trajectoryContainer)[i]; G4Track* track traj-GetTrack(); if (track-GetDefinition() G4OpticalPhoton::OpticalPhotonDefinition()) { // 强制保存每一步坐标默认只存起点终点 traj-SetStoreTrajectory(true); } } } }生成的.root或.csv轨迹文件中每行包含t, x, y, z, px, py, pz, E, volumeName可直接导入 Python 用matplotlib绘制三维光路或用pandas统计光子在各体积内的停留时间分布——这才是真正支撑论文图表的数据源。4. 避坑光学仿真中五个高频翻车点与血泪修复方案光学仿真的坑往往藏在“看起来应该没问题”的细节里。G4OpticksTest 的测试日志里沉淀了大量真实翻车现场以下是五个最高频、最隐蔽、最容易浪费一天以上的典型问题全部按「现象 → 原因 → 解决」结构给出可立即执行的修复动作。4.1 现象光子数量远低于理论预期如 1 MeV 电子在 LYSO 中应产 ~12000 光子实际只看到 200原因G4Scintillation过程默认启用kScintillationByParticleType模式但未为e-或gamma显式设置产额比例。Geant4 会回退到kScintillationPerEnergyDeposit模式而该模式下若SCINTILLATIONYIELD未设为const类型则默认值为 0。解决在材料属性表中必须用AddConstProperty而非AddProperty设置SCINTILLATIONYIELDmpt-AddConstProperty(SCINTILLATIONYIELD, 12000./MeV); // ✅ 正确 // mpt-AddProperty(SCINTILLATIONYIELD, energies, yields); // ❌ 错误会被忽略4.2 现象光子在晶体-空气界面 100% 透射完全不反射即使入射角 临界角原因G4OpBoundaryProcess要求两个相邻体积的材料都定义了RINDEX表。若空气体积G4Air未绑定光学属性表Geant4 会将RINDEX视为 1.0但G4OpBoundaryProcess::CalculateOpticalSurface()内部会因缺失RINDEX数据而跳过全反射判断。解决为所有参与边界的体积材料包括空气、反射膜、胶水显式定义RINDEXG4Material* air nist-FindOrBuildMaterial(G4_AIR); G4MaterialPropertiesTable* airMPT new G4MaterialPropertiesTable(); airMPT-AddConstProperty(RINDEX, 1.000293); // 空气折射率不能省略 air-SetMaterialPropertiesTable(airMPT);4.3 现象G4OpWLS生成的次级光子能量为 0或波长超出RINDEX定义范围导致后续传播失败原因G4OpWLS的FASTCOMPONENT/SLOWCOMPONENT波长表与RINDEX表的波长范围不重叠。例如RINDEX只定义了 200–800 nm但FASTCOMPONENT发射峰在 900 nm则G4OpWLS::PostStepDoIt()会插值得到RINDEX0进而使G4OpBoundaryProcess因折射率无效而终止跟踪。解决扩展RINDEX表覆盖 WLS 发射波段并用线性外推非截断// 在 Materials.cc 中为 LYSO 添加 WLS 发射波段如 400–600 nm std::vectorG4double rindex_lyso_wls {400*nm, 1.83, 500*nm, 1.82, 600*nm, 1.81}; mpt-AddProperty(RINDEX, rindex_lyso_wls[0], rindex_lyso_wls.size()/2);4.4 现象多线程模式下G4MTRunManager光子数随机波动且G4OpBoundaryProcess状态码错乱原因G4OpBoundaryProcess的fStatus成员变量是非线程安全的静态变量。多线程下多个光子同时访问同一G4OpBoundaryProcess实例导致状态覆盖。解决强制为每个线程创建独立的G4OpBoundaryProcess实例Geant4 10.7 支持// 在 PhysicsList.cc 中 #ifdef G4MULTITHREADED G4OpBoundaryProcess* boundary new G4OpBoundaryProcess(); boundary-UseMieHG(true); // 启用 Mie 散射模型 RegisterProcess(boundary, particle); #else RegisterProcess(new G4OpBoundaryProcess(), particle); #endif4.5 现象使用G4GenericTrap或G4ExtrudedSolid定义复杂晶体时光子在斜面边界发生“穿透”而非反射原因Geant4 对非标准体元non-primitive solid的表面法向量计算存在精度误差导致G4OpBoundaryProcess::GetStatus()误判入射方向cosTheta符号错误。解决在DetectorConstruction.cc中为复杂体元显式设置表面精度G4GenericTrap* crystalSolid new G4GenericTrap(crystal, ...); crystalSolid-SetDelta(0.001*mm); // 提高表面三角剖分精度 // 或改用 G4Trap若几何允许其解析法向更稳定5. 进阶技巧用 G4AnalysisManager 导出光子沉积热图替代人工数 hits当你要验证光纤束耦合效率或评估反射膜覆盖率对光收集的影响时靠G4HCofThisEvent里数hit数量已远远不够。你需要的是空间分辨的光子沉积密度图——即在晶体末端面或光电倍增管阴极面上每个 0.1×0.1 mm 区域接收了多少光子。G4OpticksTest 内置了一套轻量级G4AnalysisManager热图导出方案无需 ROOT 或大型后处理直接生成 CSV 格式网格数据。5.1 定义二维沉积网格按物理尺寸划分像素在DetectorConstruction.cc中为探测器末端面假设为G4Tubs的zMax面定义网格// 假设晶体末端面半径 10 mm划分为 200×200 像素分辨率 0.1 mm const G4int kNbinsX 200; const G4int kNbinsY 200; const G4double kSizeX 20.*mm; // -10~10 mm const G4double kSizeY 20.*mm; // 创建分析管理器 G4AnalysisManager* analysisManager G4AnalysisManager::Instance(); analysisManager-SetFileName(photon_deposit); analysisManager-CreateH2(deposit_map, X-Y Deposit Map, kNbinsX, -kSizeX/2, kSizeX/2, kNbinsY, -kSizeY/2, kSizeY/2, X/mm, Y/mm);5.2 在敏感度类中实时填充网格用 G4Step 获取精确沉积位置修改DetectorSD::ProcessHits()将光子沉积点映射到网格索引G4bool DetectorSD::ProcessHits(G4Step* step, G4TouchableHistory*) { G4Track* track step-GetTrack(); if (track-GetDefinition() ! G4OpticalPhoton::OpticalPhotonDefinition()) return false; G4StepPoint* postPoint step-GetPostStepPoint(); G4ThreeVector pos postPoint-GetPosition(); // 仅记录在末端面z ≈ crystal_z_max的沉积 if (std::abs(pos.z() - fCrystalZMax) 0.01*mm) { G4double x pos.x(); G4double y pos.y(); // 映射到网格索引注意G4AnalysisManager 的 h2 填充使用 bin 中心 G4int binX static_castG4int((x kSizeX/2) / (kSizeX/kNbinsX)); G4int binY static_castG4int((y kSizeY/2) / (kSizeY/kNbinsY)); if (binX 0 binX kNbinsX binY 0 binY kNbinsY) { G4AnalysisManager::Instance()-FillH2(deposit_map, x, y, 1.0); } } return true; }5.3 导出与可视化一行命令生成热图 PNGG4OpticksTest 提供plot_heatmap.py脚本自动读取photon_deposit.h2.csv并生成热图# plot_heatmap.py import numpy as np import matplotlib.pyplot as plt import pandas as pd # 读取 Geant4 生成的 CSV格式x,y,weight df pd.read_csv(photon_deposit.h2.csv, headerNone, names[x,y,w]) # 转为 2D 网格矩阵 heatmap, xedges, yedges np.histogram2d(df[x], df[y], bins[200,200], range[[-10,10],[-10,10]], weightsdf[w]) plt.figure(figsize(8,6)) plt.imshow(heatmap.T, extent[-10,10,-10,10], originlower, cmaphot) plt.colorbar(labelPhoton Count) plt.xlabel(X (mm)) plt.ylabel(Y (mm)) plt.title(Optical Photon Deposit Density on Crystal Exit Face) plt.savefig(deposit_heatmap.png, dpi300, bbox_inchestight)执行python plot_heatmap.py后你将得到一张专业级热图可直接用于论文配图。更重要的是这个流程完全脱离 Geant4 运行环境——你可以在任意机器上用不同参数重跑仿真然后用同一脚本比对热图差异快速定位反射膜缺口或光纤排列偏差的影响。从那以后我每次做光学仿真都强制在DetectorSD中加一个G4cout HIT: pos E track-GetKineticEnergy() G4endl;哪怕只是临时调试也坚持用G4AnalysisManager替代std::cout记录空间分布——因为人眼数 hits 会疲劳但 CSV 文件不会撒谎。希望帮到你。本文还有配套的精品资源点击获取
