1. 项目概述为什么要用14MeV中子轰击金刚石先说结论这套模拟解决的是“高能中子入射到金刚石碳靶后靶内到底发生了什么”这个基础探测问题。做核探测、辐射损伤评估、暗物质实验反向标定的人几乎都会撞上这一需求——你手里有一颗天然/人造金刚石前面来了14MeV中子源你需要知道反冲碳核的能谱分布、次级粒子种类、能量沉积位置甚至辐射损伤点数。拿它的原理解释清楚长久以来几个工程问题就能通通落地。14MeV这个能量不是随手选的。聚变装置里D-T反应产生的中子恰好在14.1MeV这基本上是目前实验室里能稳定拿到的“高能中子”标准能量点。模拟这个能量下中子与碳-12的相互作用准备好长期价值巨大探测器的中子标定、聚变堆第一壁材料评估、暗物质实验的核反冲刻度都要用它。从物理角度看14MeV对碳而言是高能中子的范畴。它对应的反应通道不只有弹性散射还会打开非弹性散射、次级粒子产生以及(n,α)、(n,p)类嬗变反应。你要是只用“中子打进去就弹一下”的思维来建模做屏蔽设计或用闪烁体探测器标定时误差会到离谱。用Geant4模拟首要任务就是把每条反应通道的产额、产物、能量流转逐个拆清楚。这篇内容面向的是两类读者一类是做中子物理实验的学生和工程师刚接触蒙特卡洛模拟与Geant4需要一份能跑的参考配置另一类是在做金刚石探测器、暗物质直接探测核反冲标定的同路人。你能从这篇里拿到物理过程的拆解、物理列表的选择逻辑、实际跑的配置代码以及几个我踩过才能发现的大坑。2. 物理过程拆解14MeV中子在金刚石里会发生些什么2.1 弹性散射是主角反冲碳核能谱有硬上限在14MeV能量下中子在碳靶中的主要能量损失通道是中子的弹性散射n,n。这就是经典二体碰撞问题中子质量为1碳原子核质量为12能量守恒加动量守恒求解后反冲碳核能量可以用一个简洁的截断公式估算E_R [4·A·cos²θ / (A1)²] · E_nA是靶核质量数12θ是中子散射角E_n 14MeV。代入A12你立刻得到最大反冲能量E_R,max [4×12/(13)²] × 14 ≈ 3.98 MeV注意这个数——它是金刚石探测器中“中子本地能量沉积”的硬上限。之所以强调这一点是因为暗物质实验中核反冲标定恰恰需要一条已知的反冲粒子能谱而14MeV中子在碳上的弹性反冲在3.98MeV存在锐利截止边。实测时只要看见这个截止边说明你的电子学增益和中子源能量可信很好地标定。反冲碳核在金刚石中的射程很短几乎都在微米量级内耗散所以这一谱形直接对应探测器的“核反冲能量响应”。模拟里弹性散射的全部意义在于追踪反冲碳核的完整输运。反冲碳核带几个MeV能量会持续电离、位移晶格原子最终能量在几十微米长度内散尽。宏观表现就是金刚石探测器输出一个快速脉冲电子学上看到幅值谱就对应上述反冲能谱。2.2 非弹性散射与伽马产生4.44MeV特征峰到场14MeV远超碳第一激发态阈值入射中子会容易把碳-12打成激发态这就是C-12(n,n)C-12*非弹性散射通道。碳-12的第一激发态是4.44MeV从激发态退激时会发射特征伽马射线。实际测金刚石辐照中子场时探测器周围若放一个NaI或HPGe伽马谱仪最醒目的特征是什么4259keV到4440keV之间的峰——碳的4.44MeV伽马。这条峰的强度反演后可以直接反推入射中子通量。这是中子诊断里的常规操作模拟的价值在于把非弹截面和伽马产额比精确算出来。注意这里有一个纯几何的关键天然金刚石是碳-12和碳-13的混合碳-13占比1.1%。对于14MeV中子碳-12非弹性为主通道但精确模拟时仍然要开启同位素丰度比否则几百分之一的伽马谱结构会偏。实际上Geant4的HP模型里默认不会替你掺入碳-13的丰度构建材料时要显式写同位素混合比例这一点放在物理列表部分细说。2.3 嬗变反应阈值在6MeV附近Be-9与α产物值得关注高能中子能撬开碳核这是许多人容易忽略的一件事。反应式n ¹²C → ⁹Be α质量亏损算完后这个反应的Q值大约是-5.7MeV所以入射中子动能必须超过约6.18MeV阈值考虑出射库仑势垒才能开启。14MeV落在通道开启区间它贡献的反应份额虽然比弹性散射低一个量级但它产生一个长程α粒子和一个中等质量的Be碎核。对辐照损伤而言α粒子和Be碎核的射程比反冲碳核长会带来不同的损伤形貌对探测器应用来说嬗变产物也可能沉积在金刚石内部形成长期稳定毒化。这几件事光靠“弹性散射反冲”的简单模型全都会漏掉只有完整的蒙特卡洛物理列表才能把这些次级产物轨道都追踪出来。2.4 次级粒子的级联在无损检测里表现为“二次损伤”这里要引出一个新手常做错的假设只模拟“中子与碳的一次碰撞”。实际上一个14MeV中子在靶里通常会经过多次碰撞每碰撞一次丢掉几分之几的动能如果打在薄的靶上平均碰撞次数不足一次但靶厚大于中子平均自由程时级联效应就重要了。金刚石的碳原子数密度是1.76×10²³ cm⁻³14MeV中子与碳的弹性散射截面约1.3×10⁻²⁴ cm²1.3barn左右平均自由程就是λ 1 / (N × σ) ≈ 1 / (1.76×10²³ × 1.3×10⁻²⁴) ≈ 4.4 cm一个2cm直径的圆柱金刚石靶单个中子平均只发生0.4次弹性碰撞但若用5cm厚的金刚石堆比如探测器阵列叠加平均碰撞次数就到1次以上级联效应不可被忽略。模拟就是要正确处理中子在多次散射间能量的逐渐降低过程因为几条非弹道的截面值随能量急变不能拿平均值代替。3. 建模实操Geant4环境、几何构造与物理列表选型3.1 几何与材料定义同位素丰度不要偷懒金刚石几何很简单我建议先用一个平板靶尺寸5cm × 5cm × 0.5cm密度3.52 g/cm³成分C-12自然丰度占98.9%C-13占1.1%材料定义用G4Element加同位素质量数的写法别直接用G4NistManager::FindOrBuildElement(C)草草了事。原因在上面讲非弹伽马部分已经说过——1.1%的C-13在精确伽马谱分析里并非无关紧要。原子数密度的验证计算每立方厘米含碳原子数 ρ·N_A / A ≈ 3.52 × 6.022×10²³ / 12 ≈ 1.767×10²³。这个数会在后面的截面、自由程、单位通量计算里经常用到建议模型里打印输出做自检。3.2 粒子源单能中子束的位置、方向与统计设置粒子源用G4ParticleGun最省事也可以用G4GeneralParticleSource。粒子种类填G4Neutron能量填14MeV位置放在靶面上方1cm处方向沿-Z轴轰击靶面。这里有一个统计物理的坑粒子源的均匀半径会影响注量均匀性。我在参考配置里使用半径1cm的平面圆形源让束流截面完全落在5cm×5cm靶面内这样就避免边缘漏束的几何不确定性。运行粒子数至少跑到10⁵至10⁶个中子。10⁵时弹性散射统计在薄靶条件下已经比较光滑但嬗变产物这类低概率通道需要10⁶以上才够看。跑10⁷时单线程时间开销很大下面物理列表和线程优化部分专门展开。3.3 物理列表选型QGSP_BERT_HP还是Shielding如何取舍物理列表是整个模拟里影响力最大的单一选择。对14MeV中子与碳靶我会直接剔除纯“BIC”或“BERT”这类不含高精度中子数据的构造因为它们在中子能量小于20MeV时使用的模型例如LEP在中子弹性散射截面和各向异性角度分布上精度不如带HP(High Precision)的列表。实际建议先在两种之间做选择物理列表适用场景14MeV中子-碳模拟表现优点缺点QGSP_BERT_HP高能物理、强子量能器弹性散射与嬗变产物覆盖齐全高能与低能衔接好常用于探测器响应能区数据文件大初始化慢FTFP_BERT_HPLHC量能器模拟出身与QGSP结果差异极小厚度大时粒子输运更稳健相对QGSP无明显优势Shielding屏蔽设计、剂量计算截面库较全推荐用于工程屏蔽评估低能中子的截面处理精细强子模型偏工程化与QGSP在探测器响应上略有差异个人习惯是标定谱形、分析反应产物能谱时用QGSP_BERT_HP评估屏蔽效能和次级伽马剂量时用Shielding。两个结果在14MeV碳靶上差异通常在百分之几以内不涉及敏感决策时可以互换。但这里要强调一个重点跑14MeV中子模拟HP版的数据文件G4NDL是必须装的。没有HP低能中子弹性散射截面会默认用低精度外推你将看到反冲谱在1-4MeV区间出现不正常的平坦化这是最典型的“物理列表降级”症状。3.4 打开关键过程弹性散射、非弹性、核反应、原子位移QGSP_BERT_HP 已经默认启动了中子与核的弹性、非弹过程。需要额外注意的是“反冲核”在靶内的完整输运。有些教程里人为把“中子与碳的相互作用”看成中子自己输运而反冲碳核交给用户自带的射程公式算。这不行。要在G4Step里勾选第二粒子的种类如果反冲的是碳-12离子则设置G4IonBinaryCascadePhysics或参考物理列表中的离子非弹处理保证碳离子本身的能量沉积连续可追踪。实际做DPA每原子位移评估时你还要额外打开“原子位移”的计数机制通常借用G4AtomicDisplacer或自定义Score。简单做法使用G4MultiFunctionalDetector在靶体上挂一个“能量沉积分值器PS”同时挂一个Tracker来统计所有次级粒子种类和生成位置。这样既能拿反冲碳核的能谱又能拿完整次级粒子列表一步到位。4. 记录与输出Scorer、输出文件与统计误差控制4.1 分体记录反冲能谱、伽马能谱、嬗变产额把靶体本身划分成若干层例如沿深度方向切10层每层挂一个PS记录的物理量包括动能沉积对应探测器的“反冲能量响应”反冲碳核数目与能量用G4PSNofStep或G4PSPassageCell限制只统计从碳核第一步离开碰撞格点的粒子伽马产额统计从靶体发出的所有伽马记录能量用于4.44MeV峰分析嬗变产物计数基于Step的末态粒子种类判断分别累计Be-9、α、p等这套多维记录有一大骗人的地方如果只记录总能量沉积而不看粒子种类你很难区分“中子被弹性散射产生的反冲碳核”和“中子被非弹散射产生并退激的伽马再转换为电子”这两个截然不同的过程。在暗物质探测器标定里这区别至关重要所以务必在Score里跟一条“粒子身份”属性。4.2 并行与统计误差10⁶个中子的运行时间估算单线程跑10⁶个14MeV中子轰击薄靶通常要几十分钟到几小时量级取决于机器与物理列表。Geant4多线程MT模式下可以划分8或16个线程并行每个线程独立跑一段事件最终自动合并。绝大多数情况下16线程能把10⁶事件压缩到几分钟以内。这里有两个细节MT模式里必须在RunAction初始化阶段为每个线程固定随机种子否则不会得到可重复的结果薄靶里碰撞稀疏单事件耗时极短但10⁶事件下Flux级别的平滑度才OK统计误差控制经验值10⁵事件下1-4MeV反冲区间的相对误差可以降到5%左右10⁶事件下则普遍进入1%上下。逐点报相对误差有助于判断哪一段能谱需要更多统计量——通常高能截止边和低能起始平台是最耗统计量的建议单独加大事件量而不是全局粗暴通吃。4.3 输出格式与后处理CSV/ROOT怎么选记录数据直接写CSV是最不折腾的格式。但对能谱分析来讲ROOT文件加TH1F直方图会更省事。个人的建议是CSV适合你想用Python马上处理ROOT适合你想在C环境里反复对撞谱形。CSV里每行放事件号、粒子名称、动能、位置、方向、所属过程名、来源父粒子。这样后期做反冲事件筛选时一条Python循环就能搞定import pandas as pd df pd.read_csv(output.csv) carbon_recoil df[(df[particle]C12) (df[process]hadElastic)] # 查看反冲能谱 print(carbon_recoil[kinetic_energy].describe())用ROOT时直接把每层PS数据填入TH1F命名建议带depth和quantity后缀例如h_recoil_energy_layer0。设置好Binning能谱区间0到5MeVbin宽1keV比较合适。如果想看截止边的精细结构bin宽取0.1keV也值得。5. 实际运行结果分析反冲谱、伽马谱与DPA初步估算5.1 反冲碳核能谱的截止边特征我按上述配置跑了一次10⁶中子的模拟平板金刚石靶0.5cm厚QGSP_BERT_HP。对弹性散射反冲碳核做能量分布统计后最醒目的结构就是3.98MeV截止边。实际模拟里你看不到陡峭的垂直截断而是一条带尾的斜坡下降这是因为入射束并非完全单能、碳-13的存在和次级反冲的多次叠加贡献都会把边界抹平一点。但“能谱在4MeV附近突然下跌”这个特征总归是不骗人的。这个截止边的位置对入射中子能量的响应非常敏感入射中子能量偏移0.5MeV截止边就会移动约0.14MeV所以实验标定里它本身就是一把标尺。5.2 次级伽马谱与4.44MeV峰强度伽马谱统计结果里最强的一条峰就是4.44MeV的全能峰来源是碳-12第一激发态退激。峰下方还会有约4.1MeV到4.5MeV的连续底这是反冲碳核等带电粒子初级轫致辐射和少量逃逸伽马的残余。如果你要把这条峰用于中子通量诊断就必须同时计算“每入射中子产生的4.44MeV伽马个数”——这在模拟里把PS计数除以入射粒子总数即可。参考值在不同几何下不同但0.5cm金刚石靶大致能看到每中子产4.44MeV伽马的几率在10⁻³量级。这个数量级决定了现场测量时中子源强度至少要达到10⁷ n/s以上才能抬出一个像样的峰。这是模拟对实验设计最直接的贡献——动手前先确认计数率够不够。5.3 从反冲能谱到DPA一个贴近工程的计算框架做辐照损伤评估时光看能量沉积还不够需要换算成位移原子数DPA。最简化的框架是Kinchin-Pease模型DPA N_d / N Σ(0.8 · E_dep / (2 · E_d)) / N其中E_d是位移阈能碳的金刚石一般取12-20eV模拟里能调。E_dep是每个初级反冲原子沉积在损伤级联里的能量。把反冲碳核的能谱乘以Kinchin-Pease系数再除以靶内碳原子总数得到体积平均的DPA。实战里我跑出来的量级参考1×10¹⁴ n/cm²的14MeV中子注量在0.5cm金刚石靶内产生的DPA大约在10⁻⁶到10⁻⁵这个量级。这个数看起来小但长期运行聚变装置第一壁或中子探测器窗口累积后会明显劣化金刚石的电学性能。这个估算流程是“模拟-位移损伤-工程寿命”链条里相对可靠的第一步后续要更严谨的话再用分子动力学做级联细节放行修正。6. 常见坑与排查技巧实战中反复踩过的雷区6.1 HP数据文件缺失导致的“假谱形”症状弹性反冲谱在低能端异常平坦高能截止边消失或过度平滑。排查先确认物理列表是QGSP_BERT_HP而不是QGSP_BERT再确认G4NDL环境变量指向了正确的数据目录。我第一次跑模拟时数据目录配错结果10MeV以上的弹性截面被旧模型高估反冲谱整体抬高这个坑花了两天才排查出来。6.2 统计量不足把统计涨落当物理结构症状能谱出现一组看似“等间距”的小峰仔细一查都落在同一个bin里纯属涨落。排查对每个bin算相对误差。相对误差大于20%的bin在图上看不出明显结构很多教程里“小峰”都是这个来源。解决办法很简单把该区间的bin合并或对重点能区加算事件数。不要一上来就换物理模型。6.3 几何优化不够多线程也救不了耗时症状10⁶事件的运行耗时超过你午休时间。排查对靶体外的世界体积做G4Region裁剪把探测区只限制在金刚石周围关闭空气里的中子物理细节确认每次事件里从G4Navigator返回的几何步数是否过大。在靶体外包一个半径10cm的“死区”把粒子杀灭在范围内对速度提升立竿见影。再出不了才考虑换物理列表到Shielding它通常在低能截面上做了一点优化初始化速度也更快。6.4 随机数种子不一致导致无法复现症状同一配置跑两遍结果却有一定漂移。排查在RunAction里给每个线程用master_seed thread_index的方式重设G4Random。另外改变线程数也会改变结果因为事件分配到各线程的批次不同。如果想严格复现必须固定线程数并保存G4RandomStatus。实测量能谱需要多日大数据采集比对时这些细节会成为硬需求。7. 结尾要说的话这套模拟还能怎么接续扩展一个人实际跑过之后才能感受到Geant4给这个14MeV-碳靶问题提供的答案不只是“一个能谱图”而是一整套可插拔流程换靶材、换能量、换物理量输出框架都不动。我这套配置直接改成SiC靶只要在材料定义里加上硅元素反冲能谱就自动覆盖硅反冲成分改成12MeV入射截止边自动移下来很好用。接下来准备拿来做的扩展是“伴生伽马时间谱”给靶体加一个时间切片PS统计中子入射后0-100ns内伽马到达探测器的分布这能直接对接未来快中子诊断系统的波形仿真。同样的框架再挂一个G4RadioactiveDecay物理列表还能进一步追踪嬗变产物Be-9的衰变链评估长寿命放射性积累。这些扩展本质上都是同一个模拟骨架的小改动14MeV轰击金刚石只是开了个头后续能跑的方向还很多。