FLUKA实操——质子治疗瞬发辐射剂量计算
国内计算瞬发辐射剂量通常用最大能量和最大流强来计算瞬发辐射的剂量。这种计算方式虽然最为保守，但对于降低总成本不利，尤其是质子治疗中心未来可能出现建设热潮。当前已设计建造的一些质子治疗中心，有的主屏蔽体混凝土厚度达5~6米，加上砌块有些部位甚至达到7~8米，超大体积混凝土对浇筑工艺提出新的挑战和困难，对预埋管线也有着诸多不利因素。为使质子治疗的社会效益最大化，本文建议按照临床实际情况进行设置束流损失，进而计算瞬发辐射剂量。
利用FLUKA等蒙卡软件进行辐射防护设计和屏蔽计算，相比于经验公式，可以相对准确的给出剂量大小和全局三维的剂量分布，这依赖于制作出合理简化，束流模型和屏蔽体模型。这有利于人们优化屏蔽体，降低辐射防护的固定成本。本文将深度讲解如何用FLUKA和Flair获得瞬发辐射的剂量分布。选用的实例是IBA多室系统中的一个旋转治疗室，如下图所示。
先简单介绍一下旋转治疗室及其工作流程，将有助于模型简化。1985年，Loma Linda大学第一次使用旋转机架进行多方向照射，自从该系统发展以来，旋转机架已成为质子治疗系统的标准功能。旋转治疗室可以对患者进行360°多射野照射，患者体验好，是很受欢迎的一种配置。即便是现在的小型化质子治疗系统，Gantry也是必备的，虽然有些只能旋转一定的角度。下图是Gantry的一个示意图。
一次治疗中，患者接收多个方向射野的照射治疗，一般为2~4个。每个照射野治疗头会根据计划系统的处方剂量向加速器请求能量和流强的质子束流。照射过程中，所有的质子会阻停在患者体内，沉积能量杀灭癌细胞。同时也产生大量的瞬发辐射。
由于旋转治疗室的束流传输系统质子损失非常小，这部分束流的磁铁和真空管道以及治疗头可以忽略，因此只包含屏蔽体的几何模型，如下图所示，图中显示的是X-Y平面。可见还是比较复杂的。
前期的FLUKA系列文章讲述过，对于复杂的几何模型，可以采用SimpleGEO来进行建模。大型质子治疗中心医院的建筑非常复杂，虽然一些与辐射无关的结构可以被忽略或被简化，但使用FLUKA直接进行几何建模仍比较困难，而使用NX、Solid works、Inventor等CAD软件制作的3D模型与FLUKA格式不同，无法互用。因此Flair和SimpleGEO是必要的有力工具。下图是Flair制作的几何模型，分别是X-Z和Y-Z平面。
由于照射野通常为2~4个，本例取4个，分别为垂直方向两个射野，如Y-Z平面所示，水平两个射野，如X-Z平面所示。
在很多环评报告中，大家通常用最大能量和最大流强来计算瞬发辐射的剂量。这种计算方式虽然最为保守，但也简单粗暴，对于降低总成本是不利的，尤其是质子治疗中心未来可能出现建设热潮。为使质子治疗的社会效益最大化，本文建议按照临床实际情况进行设置束流损失。
临床治疗时，一般不会用最大单一能量和束流来治疗，质子束的能量和流强根据tps制定的处方剂量分布得到，这个转换过程根据系统的供应商而定。同时也要考虑到病种的分布，不同地区不同医院各不相同。病种的不同带来质子能量深度的不同，几乎不同患者也是不同的。所以考虑实际情况比较复杂，需要供应商和医院统计较长时间，才能得到质子治疗所需要的能量和流强的一个真实统计。
为简单起见，本文取文献[1]表1中的数据来进行计算，红框的部分。可见160 MeV左右能量的质子束使用最多，权重为0.313，230 MeV使用得比较少。流强由第二行得出，一年的总强度是211.76nA.h。
表1 束流损失能量统计
下面要根据上表计算转换系数，以130 MeV为例，一年的总强度是：
130 MeV的质子数为：
FLUKA中辐射剂量用Dose-EQ，默认的单位是pSv每核子，转化为uSv/年就需要乘以一定的系数，这就是转换系数Normal Factor，在Flair出剂量分布图时需要这个系数才能把计算值转换成物理量，这也是初学者比较难搞清楚的地方。
注：当你的统计单位不是nA.h，这里的计算公式不适用，仅做参考。
计算出所有能量对应的转换系数，如下表最后一列。
表2 转换系数
屏蔽计算并不需要很精确的束流参数，可以简单地设置一下束流参数，如下：
由于本例中有两个射束方向并不是沿轴线，需要设置射束的角度，在BEAMPOS中设置Cosx=cosd(60)，并检验一下设置的是否是想要的方向，如下图
注1：有余力的读者也可以直接用Useroutine一次性考虑权重，无需权重叠加出图，可以节省很多计算时间和数据处理的工作量，后期有机会再讲； 注2：输入文件中其他需要注意的地方，如果遇到几何体需要旋转的，可以用transform，但不能嵌套，只能有一层；
注3：在区域REGIN的设置中，不推荐用括号()来处理几何体，尤其是两层括号嵌套后FLUKA极有可能不识别，但是嵌套的括号()对设计几何结构带来极大方便，在不得已的情况下最多用一层；
注4：通常屏蔽体材料是普通混凝土CONstan，成分配比可以参考下图设置。另外，重混凝土CONNiron的材料配比仅供参考。
GTR1每个能量朝90°照射野对应的瞬发剂量分布如下图，可见随着能量增大，瞬发辐射剂量逐渐增强。
特别注意，Flair作图的过程中就需要第3节中计算的转换系数，填入Norm中，如下图所示。
然后，对于每个方向用Gunplot在把能量点进行叠加，这里的叠加是按照表1给的权重进行，具体的代码如下：
最后一行代码中的.dat文件是每个能量在用Flair作图时生成的，包含了几何和剂量信息。权重叠加的结果如下图，介于130 MeV和230 MeV单能分布结果之间。
按照以上步骤，重复计算其他三个方向。
最后，将四个方向等权重叠加，结果如左图(a)，MCNPX的验证结果如右图(b)，可见两个不同蒙卡程序获得的结果符合的很好。
为减少运算时间，可以采用计算机多核并行的方法，这里针对多服务器多核的小型Cluster做简单的介绍。由于FLUKA本身并不具有并行功能，需要用户自行开发并行运算程序。基本原理就是将大粒子数任务拆分成多个小粒子数任务，分配给多个计算机核心同步完成，这样计算时间将随计算机核心数增加而反比例减少。示意图如下：
例如用10个刀片服务器，每个服务器有4个核心。总共40个核心，理论上总计算时间就能减少到单机运算时间的四十分之一。
对于并行运算，也相应要开发运算流程，这个流程将保障因为粒子数拆分不会影响最后的结果。当然如果有需要也要对能量，源项等进行拆分，分步运行，最后以权重叠加。
注意，并行系统并不会自动做数据处理，需要将得到的计算原始数据传回主机，用Flair或者其他程序进行进一步分析才能出结果。以上结果采用的是小型Cluster并行运算。
在很多环评报告中，通常用最大能量和最大流强来计算瞬发辐射的剂量。这种计算方式虽然最为保守，但也简单粗暴，对于降低总成本是不利，尤其是质子治疗中心未来可能出现建设热潮。为使质子治疗的社会效益最大化，本文建议按照临床实际情况设置束流损失。IBA和Varian这类成熟厂商有比较全面的数据，但国内新兴的企业和研究所数据缺乏，可适当借鉴国外的数据先行计算转换系数。