14MeV中子轰击金刚石的Geant4模拟全流程解析

“14MeV轰击金刚石”这个题目第一反应不该是“打开Geant4就开跑”而是先想清楚你想从模拟里拿到的是探测器的脉冲高度谱还是材料辐照损伤的初级离位原子PKA分布这两个目标对应的输出量完全不同物理列表的侧重点也不同。我自己做聚变中子诊断探测器模拟那会儿在这个问题上绕了很大一圈。这篇就按我完整复现的路线写从环境配置、几何构造、物理列表选型、粒子源设置到最终谱图解读全部拉通中间会穿插截面估算和几个实打实的坑适合做金刚石中子探测器、CVD金刚石辐照实验以及快中子能谱分析的同学直接参考。1. 14MeV中子在碳材料里的物理图景先搞清楚你在算什么14MeV这个能量不是随便挑的。氘氚聚变反应D T → ⁴He n放出的中子动能约14.1MeV所以磁约束聚变装置、惯性约束点火装置以及便携式中子发生器出射中子基本都是这个单能峰。单能快中子有个好处反应道清晰便于模拟结果和ENDF/TENDL核数据逐项对比。金刚石这边单晶CVD金刚石因为载流子迁移率高、禁带宽度5.5eV、暗电流极低是做快中子飞行时间谱仪和闪烁体替代方案的常见材料。但14MeV中子在碳上不是只发生一种作用这点很多人建模时没意识到。14MeV中子在碳-12上的主要反应道可以分成三类弹性散射n ¹²C → n′ ¹²C基态。出射的反冲碳核带电在金刚石晶格中电离这是探测器最主要的信号来源。碳的A12弹性反冲最大动能T_max 4A/(A1)² × E_n代入可得约3.98MeV这个数值后面看谱时会经常用到。非弹性散射n ¹²C → ¹²C* n′¹²C第一激发态退激发出4.44MeV的伽马。非弹的阈能约4.8MeV14MeV下截面已经比较可观。核反应¹²C(n,α)⁹Be的Q值约-5.7MeV反应阈能约6.5MeV14MeV时能顺利发生¹²C(n,n′)3α也在该能量区间打开¹²C(n,p)¹²B等通道弱一些但在精确模拟时不能完全无视。这些反应道放在模拟里意味着什么呢如果只把金刚石当成一个“能量沉积体”来算你会丢掉所有反应产物信息如果只记录总沉积能你无法区分弹性散射和(n,α)反应对谱形各自的贡献。所以我的做法是先明确输出量。做探测器模拟输出Edep谱、反冲碳核能谱、次级粒子能谱做辐照损伤评估输出PKA能谱和位移损伤截面积分。这两套目标在Geant4里就是同一个程序的不同输出分支但一开始就要在SensitiveDetector的字段设计上区分开。额外提醒一点金刚石本质上是碳低Z对伽马探测效率天然低。4.44MeV伽马很可能直接逃逸出毫米级晶体在沉积谱上只会留下康普顿平台的尾巴。这个物理事实决定了实验谱和模拟谱对比时伽马分量应该很小如果你模拟出了大片伽马峰那基本是几何或物理列表设置出了问题。2. 环境与数据准备G4NDL是14MeV模拟的灵魂Geant4本身是工具包不是开箱即用的软件。版本上我建议直接用11.x系列我本地的实测环境是11.210.7也能跑但老版本在ROOT输出和多线程合并机制上不够省心。编译前想清楚要不要OpenGL可视化——批量跑数据的场景下可视化模块纯属拖累cmake阶段关掉能省很多依赖mkdir build cd build cmake -DGeant4_DIR/path/to/geant4/lib/Geant4-11.2 \ -DGEANT4_BUILD_MULTITHREADEDON \ -DGEANT4_USE_OPENGL_X11OFF \ -DGEANT4_INSTALL_DATAON \ .. make -j4数据文件是14MeV中子模拟的命门。Geant4的数据包里G4NDLNeutron Data Library必须完整下载它包含了从热中子到20MeV的中子诱导反应截面数据正是HP高精度中子输运模型的底料。14MeV落在它的能量覆盖范围内所以弹性散射、非弹性散射、(n,α)、(n,p)这些反应道都由G4NDL里的ENDF/B、JEFF等评价核数据驱动。如果只装默认数据包而漏了G4NDL程序启动时会打出一堆警告然后中子直接“穿墙”所有反应都是零这个坑我见得太多了。安装数据后建议显式导出环境变量避免cmake自动配置的路径在换机器后失效export G4NEUTRONXSDATA/opt/geant4-data/G4NDL.4.7 export G4ENSDFSTATEDATA/opt/geant4-data/G4ENSDFSTATE.2.3 export G4PARTICLEXSDATA/opt/geant4-data/G4PARTICLEXS.1.1 export G4PHOTONEVAPORATIONDATA/opt/geant4-data/G4PhotonEvaporation.5.7一个小经验跑任何中子模拟前先跑一个100个中子的最小测试看第一个相互作用发生在哪个过程。用G4Step的GetPostStepPoint()-GetProcessDefinedStep()打印过程名确认有弹性散射和非弹反应出现再放开了跑大数据。这一步能替你排查掉80%的“模拟结果全为零但不知道哪里错”的情况。3. 几何构造与敏感探测器尺寸选择不是随手填的金刚石体块我建议直接用G4Box典型的单晶CVD探测器尺寸是5mm × 5mm × 1mm。为什么是1mm从宏观截面看金刚石原子密度约1.76×10²³ atoms/cm³14MeV中子对碳的总截面大致1.2~1.3barn宏观截面约0.21/cm平均自由程约4.7cm。换句话说1mm厚晶体一次穿越被中子击中的概率大概2%。这个数对探测器效率概念很关键——金刚石探测器对快中子本质上是低效探测器厚度增加对效率贡献是线性的。材料定义用G4NistManager简洁直接auto nist G4NistManager::Instance(); auto carbon nist-FindOrBuildElement(C); auto diamond new G4Material(diamond, 3.515*g/cm3, 1); diamond-AddElement(carbon, 1);世界体我用真空G4_Galactic避免空气对出射低能带电粒子产生不必要的电离能量。敏感探测器的实现核心就一个类class DiamondSD : public G4VSensitiveDetector { public: DiamondSD(G4String name) : G4VSensitiveDetector(name) {} G4bool ProcessHits(G4Step* step, G4TouchableHistory*) override { auto edep step-GetTotalEnergyDeposit(); auto track step-GetTrack(); auto particle track-GetDefinition()-GetParticleName(); auto process step-GetPostStepPoint()-GetProcessDefinedStep() ? step-GetPostStepPoint()-GetProcessDefinedStep()-GetProcessName() : G4String(None); // 填充ntuple: edep, particle, process, posX, posY, posZ return true; } };这里我给每个step记录的字段包括能量沉积、粒子名、过程名、沉积位置。为什么要记粒子名和过程名因为后续分析时要按“反冲碳”“alpha”“伽马”分类切片。如果你只记一条“总能量沉积”的直方图后面想拆反应道就得重跑一遍极浪费时间。从第一次运行就尽量把谱系信息存全这是模拟程序设计的习惯问题。如果还关心位置分辨或损伤分布可以把金刚石切成多层薄片每层一个逻辑体分别绑定SD。但注意切片数越多step处理开销越大初学者先整体一块跑通再考虑细分。4. 物理列表选型FTFP_BERT_HP和QGSP_BIC_HP到底该选谁物理列表是模拟结果可信度的根基。14MeV中子必须在高精度中子输运HP模式下跑Geant4里最常用的两个参考物理列表是FTFP_BERT_HP和QGSP_BIC_HP。我实测下来这个场景FTFP_BERT_HP更合适。FTFP_BERT_HP的构成是20MeV以上用FTFPFritiof弦模型处理高能强子相互作用20MeV以下切换到HP高精度中子输运低能区用BERTBertini级联作为补充。QGSP_BIC_HP则集成了Binary Cascade更多面向离子束治疗模拟对重离子的次级粒子细节有优势但在中子与轻核碳的反应道覆盖上FTFP_BERT_HP是文献里更“标准”的选择。取物理列表的过程很简单#include G4PhysListFactory.hh auto factory new G4PhysListFactory(); G4VModularPhysicsList* physics factory-GetReferencePhysList(FTFP_BERT_HP);还需要显式加上电磁过程吗FTFP_BERT_HP默认包含标准电磁过程包反冲碳核、alpha、质子、伽马产生后的电离和输运都会被处理。不需要额外添加。有一个概念必须说清楚HP模型的精度上限是20MeV。14MeV正好落在G4NDL数据覆盖区间所以弹性散射的角分布、非弹性散射的激发函数、(n,α)反应道的截面都由ENDF评价数据直接驱动。如果你手滑选了FTFP_BERT不带HP后缀14MeV中子会走参数化模型截面和反应产物基本没法看。这个后缀值几个小时的排查时间值得刻在脑门上。还有一件事容易被忽略HP模式下低能中子的输运速度完全取决于截面数据查表的效率14MeV单能束还好如果是宽谱中子源跑起来会明显偏慢。这时关闭所有与中子无关的可视化、调低输出频率是基本操作后面会再展开说。5. 粒子源设置与计算规模事件数怎么定才够统计意义粒子源在这个题目里就是一个G4ParticleGun没有太多花样auto particleGun new G4ParticleGun(1); auto neutron G4ParticleTable::GetParticleTable()-FindParticle(neutron); particleGun-SetParticleDefinition(neutron); particleGun-SetParticleEnergy(14.0*MeV); particleGun-SetParticlePosition(G4ThreeVector(0, 0, -0.6*mm)); particleGun-SetParticleMomentumDirection(G4ThreeVector(0, 0, 1));想模拟束斑分布就在Position里引入随机量例如在半径2mm的圆内均匀取样生成x、y再传给Gunauto r 2.0*mm * std::sqrt(G4UniformRand()); auto theta 2.0*M_PI * G4UniformRand(); G4double x r * std::cos(theta); G4double y r * std::sin(theta);事件数不能拍脑袋。回到第3节的截面估算1mm厚金刚石对14MeV中子的作用概率约2%。跑1×10⁶个中子真正发生核作用的只有约2×10⁴次其中弹性散射占大部分(n,α)可能只占百分之几。如果你要的是反冲碳谱的平滑曲线1×10⁶事件可以做要单独看弱反应道1×10⁷不嫌多。多线程不要一上来就开满。我习惯先用单线程跑2000事件确认程序稳定性和截面数据正常再用MT模式跑全量。Geant4多线程下的随机数种子是每线程独立的事件分配也由框架完成不需要自己写人工划分。提速的一个小技巧在于生产阈值production cut)。默认的低能电磁阈值会让低能伽马和电子产生大量无意义step对中子诱导反应的结果影响很小。如果只关心碳反冲核和核反应产物可以把cut值设到0.1mm甚至0.5mm来砍掉软伽马和低能电子速度能快30%~50%。但如果你的目标是精确模拟探测器脉冲高度谱这一步要谨慎因为低能沉积被切掉会直接压低谱的低能段。6. 数据输出与结果解读反冲谱、α谱和4.44MeV伽马我用G4AnalysisManager直接输出ROOT文件。推荐的ntuple字段设计如下字段类型用途edepdouble每个step的能量沉积particlestring当前step的粒子名称processstring产生该step的过程名posX/posY/posZdouble沉积位置trackIDint关联粒子轨迹跑完后重点看三个东西。第一反冲碳核能谱。弹性散射的碳反冲能量从0到约3.98MeV连续分布。由于质心系散射角分布不是各向同性谱形不会是一条水平线而是在低能段偏高、高能端缓慢下降后在3.98MeV附近出现截止。这个截止值是对14MeV能量的直接验证——如果截止位置明显偏离先检查入射能量是否真的设成了14MeV。碳反冲在金刚石探测器里的信号是最主要的对应实验上就是快中子引起的核反冲脉冲。第二(n,α)反应产物。¹²C(n,α)⁹Be的Q值约-5.7MeV反应释放的动能使得alpha粒子能量在几个MeV区间。把ntuple里particlealpha的step单独挑出来画能谱能看到一个宽峰结构。这个alpha信号在探测器里沉积效率接近100%因为alpha射程远小于1mm晶体厚度。但注意和反冲碳信号区分alpha的粒子径迹电离密度和碳不同实验上常利用脉冲形状甄别而模拟里直接用粒子种类字段切片即可。第三4.44MeV伽马。非弹性散射退激伽马虽然能量高但在毫米级金刚石里沉积概率很低。因为金刚石是低Z材料光电吸收截面小伽马主要以康普顿散射方式损失能量留在晶体里的往往只有几十到几百keV的电子能量。所以你在Edep谱里会看到一个小而缓的康普顿平台而不是尖锐的光电峰。这个平台就是非弹成分存在的指纹想用来标定的话得把探测器做厚或加高Z包壳纯金刚石很难直接看到4.44MeV峰。有一点必须特别强调实验脉冲高度谱不能直接等于模拟的Edep谱。金刚石探测器对能量沉积的响应有载流子产生统计涨落Fano因子、电荷收集不完全、电子学噪声等影响模拟峰通常比实验峰窄。严谨的对照方法是在模拟Edep谱上叠加高斯展宽展宽参数由实验噪声水平决定。很多同学第一步就栽在这里把模拟谱和实验谱直接叠图然后怀疑物理列表错了实际只是少了展宽这一步。7. 实测中的三个坑数据路径、多线程合并与低能cut这段把实际操作中真正浪费过我时间的坑拎出来说比任何教程里的“注意事项”都有价值。第一个坑G4NDL数据路径不对程序不报错但结果全错。现象是运行日志里能看到PhysicsList加载了HP模型但中子打到金刚石上全部穿透Edep全是0。原因大概率是G4NEUTRONXSDATA环境变量没生效或指向了不完整的数据目录。检查方法很简单在程序的初始化阶段打印一遍中子与碳的首个相互作用过程如果全是Transportation而没有Elastic直接去查环境变量。这个坑的隐蔽性在于它不像“找不到文件”那样被系统直接拒绝Geant4会在数据缺失时静默降级到参数化模型。第二个坑MT模式跑完ROOT输出文件里一堆空直方图。新版G4AnalysisManager在多线程下会自动处理各线程的直方图合并但有一个前提你必须在EndOfRunAction里正确调用WriteFile()并且不要在Worker线程里直接操作主线程的ROOT文件。我遇到的情况是RunAction里忘了写WriteFile导致只在部分线程触发时输出看起来就是数据“丢了一半”。解决办法是先在初始化里显式设置输出格式和文件名再在RunAction的EndOfRunAction里统一调用G4AnalysisManager::Instance()-WriteFile();跑小规模测试时顺便验证一下ROOT总计数是否等于粒子源发射的总数能快速发现合并问题。第三个坑production cut设置得过激进把真实信号也切没了。我做剂量评估时为了提速把cut设到1mm结果反冲碳的低能部分能量沉积被大幅低估。原因是低能碳核和alpha的射程本来就只有几微米而cut参数主要限制的是伽马和电子但间接影响了电磁簇射的后继沉积。教训是cut值可以优化但不能脱离你的物理目标乱调。建议分两个阶段推进——先不调cut跑完整参考结果后面做参数扫描时再研究cut对速度与精度的影响比。最后说一句个人体会这个模拟的核心与其说是把Geant4跑通不如说是在跑通之后能否把每个物理量对应到实验可测信号上去。我强烈建议在正式大规模模拟前先用1000个事件跑一遍把每个粒子产生过程名打印出来对照ENDF截面逐项确认反应道存在。这个习惯帮我排掉了至少五个看起来很玄学的问题比如非弹伽马完全消失、反冲碳高能端截止偏移、alpha计数异常偏少等等。养成这个“最小验证”的习惯再复杂的模拟也能少走一大半弯路。