九州大学学術情報リポジトリ

全文

(1)

九州大学学術情報リポジトリ

Kyushu University Institutional Repository

コウシトブッシツノソウゴサヨウ二カンスルモンテ カルロ・シミュレーションノプログラミング

上原, 周三

九州大学医療技術短期大学部一般教育助教授

https://doi.org/10.15017/136

出版情報:九州大学医療技術短期大学部紀要. 11, pp.1-16, 1984-03-01. 九州大学医療技術短期大学部 バージョン:

権利関係:

(2)

光子と物質の相互作用に関するモンテカルロ シミュレーションのプログラミング

上 原 周 三*

A Programming of Monte Carlo Simulation

for the lnteraction of Photons with Matter

SHUZO UEHARA

1 緒言

 モンテカルロ法は,乱数を取扱う技法の総称 である。確率的事象に対して,とりわけ有力な 方法であり,粒子の輸送現象の解析にしばしば 用いられる。この方法を用い,実験条件を模擬 して計算機によって実験を行なうことを,モン テカルロ・シミュレーションという。

 光子についてのモンテカルロ計算は,これま でにも多数なされている。4 5 6 8 11 13)この 方法の特長は

①実験的に求めることが不可能な現象や,複 雑な実験システムの場合でも,容易に結果を予 測することができる。

②光子と物質の相互作用の結果生じるあらゆ る現象に対して,多方面に応用できる。

③特定分野の専門的知識がなくても,ごく基 本的な事柄だけから.解を得ることができる。

たとえば線量分布計算がその例である。

などを挙げることができる。

 一方欠点は.統計的ゆらぎを避けることがで きないので,多数回の試行を繰り返すことにな り,そのため演算時間が長くなる点である。し かし,大型計算機が容易に利用できるようにな った現在,これは大きな問題ではなくなった。

 放射線治療や診断の物理的側面は,放射線と りわけ光子と物質の相互作用そのものであり,

*九州大学医療技術短期大学部一般教育助教授

この方法の良い対象である」5)これらに関連す る現象についてs)モンテカルロ・シミュレーシ ョンのFORTRANプログラムを作成した。 大 筋のフローチャートは,文献の多くに示されて いるが,細部は明らかにされていないので,ア ルゴリズムや手法を独自に開発しなければなら なかった。プログラミングの段階で工夫を要し た点,処理方法,サブルーチンなどの技術的問 題を皿で述べる。すでに5種類の異なる課題に ついて,完結したプログラムを作成しているが,

これらの詳細と計算結果は別稿に譲る。本稿で はIVにおいて,若干の留意点を述べるにとどめ

る。

且 TSS端末

 一般に,長い演算時間を要するモンテカルロ 法は,高速で記憶容量の大きい大型計算機によ

って計算される。

 九大大型計算機センター・FACOM M200 およびM382のTSS(Time Sharing Sys−

tem)利用によって,プログラム開発・演算を 行なった。ホストコンピュータ(大型)に300

ボー(ビット/秒)の交換回線を接続し,音響 カプラ,RS232 Cインターフェイスを介して,

TSS端末にはパソコンMZ80Bを用いた。シ

ステム図をFig.1に示す。端末の送受信制御の ためのアセンブラプ・グラム17)は,フ・。ピー ディスクを用いたファイル転送こそ未完成だけ れども,スクリーンエディット機能,プリンタ

(3)

    光子と物質の相互作用に関するモンテカルロ・シミュレーションのプログラミング

制御機能を備えており,MZ80Bはインテリジ   の活用によって,大型計算機の有効利用と作業 エント端未として働く。これとTSSコマンド7)  能率の大幅な向上が実現できた。

九大計算機

ンター   電話局       音響カプラ型計算機

電話

MZ80・B

         RS232C

Fig.1 TSS端末システム

PRINTER

皿 共通に用いられるルーチン及び手法

皿.1 乱数

 モンテカルロ法では,質のよい乱数が要求さ

れる。〔0,1〕の一様乱数には,FACOM科

学計算ライブラリSSI, fiよりRANU2を利用し た。このサブルーチンの出力は一般にN個の要 素からなる1次元配列として得られるが,N=

1と与えておけば.CALLするごとに配列と同 じ順序で1個つつ発生させることができる。こ こで作成したプログラムにおいて,RANU 2が 呼び出される回数は,1光子ヒストリーにつき

20−30回に及ぶ。

 RANN2は正規分布をなす乱数を発生させる サブルーチンであり.検出器のレスポンスを求 める計算において,検出器分解能が有限の場合 のサンプリングに使用される。

 なお,これらのサブルーチンにおける初期パ

ラメタIX、 Nは, RANN2がRANU2を基に

導出されるので,共通に用いられる。

皿.2 座標系と初期条件

 通常,3次元カーテシアン座標系(x,y,

z)を設定し,これを物質固定系と称する。原 点は実験条件に応じて,たとえば線源の位置,

ファントムの中JD ,あるいは検出器の表面など に置く。いったん.座標系を決定すれば,実験 条件に対応させて,寸法.位置などを入力する。

cgs単位が一般的であるが,エネルギーはkeV 単位が便利である。

 光子の進行方向の,各座標軸の正方向につい

n z

NN NN

NN

      ロ

        e l

      :

      o  m

   l   ・

       争\ 1  ,

x隔簡一一_\レ/

Fig.2 座標系と光子の進行方向

y

ての方向余弦を(1,m,n)とする。 Fig.2に 示すように,光子の進行方向とz軸のなす極角

θ,xy面への射影の方位角qとおけば

   1 == sine cosq

   m==sine sin q (1)

   n= cos e

である。光子が等方的に放出される場合の,方 向余弦のサンプリング法について,以下に示す。

 方向余弦のz成分は

   n=2e 一1, (2)

また,方位角は,

   q==z(2rp−1) (2r

となる。(1),(2),(2)ノ式より(1,m,n)は乱

数ξ,ηによって表わされ,サンプりングされ

アこことになる。

 ここで,一様な強度分布を持つ,半径aの円

(4)

上 原 周 奇 形線源が(x,y)面にある場合を考える。出

発時の進行方向は,上記のとおりである。一方,

出発地点は,円内の各点を一様にサンプリング する1)ことによって得られる。ξ,ηを新たな 乱数とすれば

   x=aff cos T(2n一一1)

       (3)

   y=:aff sin x(2ny一一1)

となり,a=0ならば点線源を表わす。ヒスト リーを始めるときは,この例のように出発地点 と進行方向の初期値を与えておかなければなら

ない。

皿.3 自由行程 1.境界がない場合

 光子が物質中を進むとき,反応を起こす点ま での飛行距離を自由行程といい,乱数ξを用い

    R=m Jne/a

      {4)

      tot で与えられる。3)

σ totは物質の全反応断面積をcm司単位であら わした線吸収係数で.したがってRはcm単位で 得られる。全断面積は光電効果σphoto.イン コヒーレント散舌Lσincoh・ コヒーレント散舌L σcohおよび電子対生成σpairの和;

 atot == a photo + aincoh + acoh + apair        (5)

である。

 1keV−8MeVの範開のエネルギーにおけ

る断面積のデータはStorm, Israe118)による テーブルを用いた。このデータは,別の用途に も使用することが多い基礎データなので,随時 読み出せるよう,元素記号を冠したデータセッ

トにファイルした。

 テーブル中のエネルギーの中間の値,Eに対 する断面積σは,両対数内挿により求めた。す なわち・テーブルのエネルギーをEXi_1とE X.,その中間のエネルギーをEとすれば  l  r = ( log E−log EXi−1 )/ ( log EXi     一 log EXi−1)

a =:1 o ( log ai ml +r   (log a,一 log ai 一一1) )

       (6)

 ところでStorm, IsraelのデータはBarn/

Atom単位で与えられている。(4)式の計算のた めに,これをcm  1単位に変換する必要があっ

た。

 一般に物質が分子の場合について,単位の変 換法を述べる。1モルの質量をM(g),密度

をρ(g/cm3),化学式のj番目の元素の原子 数をfjとすれば, その元素についての換算定

数Kは

  Kj=O.6023・p・fj/M (7)

であり,B/A単位で表わされている断面積に これを掛ければ,cm 一1単位に変換することがで きる。いくつかの物質についての換算定数を

T able 1に示す,,

T able l 種々の物質の換算定数;K matter ρ(9/6m㊧ M(9)

K

Ge

5.4 72.59 0.04481

NaI

3.67 149.89 NaO.01475

P 0.01475

water iH20)

1.0 18.02

H O.06687 O 0.03343

 lucite

iC5H802) 1.18 100.12

C O.03549 g O.05679 O 0.01420

 相互作用の型および元素の種類をそれぞれi,

jで指定すれば・分子の全断面積・囎{(cm・ 1)

a 111.Oti(cm一 ) = ii. 17. Kj a 1・ (B/A) (8)

によって得られる。(8}式においてj == 1の場合,

左辺は単原子物質の線吸収係数になる。

2.境界がある場合

 Fig.3のように.物質1における自由行程が 長く,境界を越えて異なる物質皿に入る場合に

は,(4)式はそのまま適用できない。そこで,光 子の進路と境界面の交点を求め,そこを基点に.

あらためてllにおける自由行程を(4)式によって サンプリングする。ただし,エネルギーと方向 は,1における値を保持している。

となる。

(5)

光子と物質の相互作用に関するモンテカルロ・シミュレーシ・ンのプログラミング 1

(Xo・Yo・Zo)

n

SN(X,,YA・ZA)

f(x,y,z) = O

Fig.3 異なる物質問での自由行程

 明瞭な境界がある場合,境界面の一般形はx,

y,zの2次関数

   f(x, y, z)==O (9)

で表わされる。一方.(X。,y。,Z。)を出発 した光子が,行程R進んだときの位置は    x= xo 十R l

   y == yo 十R  rn    z 一一 zo 十R n である。

(10)

 (9),(10)式よりRを消去して解けば,入射点

(XA,yA,ZA)が求まる。いくつかの具体

例.を次に示す。

a)傾きのない平面(Fig.4a)

   ZAM ZI

   ・A一・・+÷(・バ・・)(11)

   y。一y・+昏(・べ・・)

b)傾いた厚さdの平板(Fig.4b)

 板の上面との交点は

    一 nd 十 2 s in ct (1zo  nxo )

       (12)

   ZA 一 rm/sin ct 一n cosa

XA,yAは(11)式と同様である。

c)3次元曲面(Fig。4c)

 Zについての一般解は

   、土一B土穿≡亙(13)

の2交点になる。ただし,B2−AC〈0の場合 には交わらず,1での進行を続ける。単純な境 界面については,解析的な解があり,A,B,

Cは次の通りである。

c−1)円柱(x2+z2 ・a2,一。。〈y〈・。)

  A= 12十 n2

  B=: 12 zo 一ln xo (14)

  C =(1zo )2−2 ln xozo +n2(x20 一 a2)

・一2)円柱(x2+z2一 a2・邊≦・≦暑)

(145式と同じだが,高さが有限なので交わらな いことがある。また,底面から入射する場合の 交点は

  Y, == ±h/2

       1  ・。=・・+正(y。一y・)

        

  ・A =:・・+『五(y。一y・)

となる。

c−3)球(x2+y2+・2=a2)

  A=1

  B= (12十m2)zo−n(lxo十myo)

(15)

  c=(12 十m2) z20  2n zo(lxo 十my.o)

    十n2 (x20 十 y20 −a2 ) (16)

c−4)放物面(x2/a2+y2/b2=2z)

  A  12b2 十m2 a2

  a) 1 II

b)

Z= Zl

x

1

α

d

H

c)

x

1

Zか

z

II

z

Fig.4 光子と境界面しの交点

z

z

y

(6)

       上 原

  B = 12b2 zo 十 m2 a2 zo 十 n2 a2 b2

    −ln b2 xo−mna2 yo (17)

  C= b2 (1zo mnxo )2 十a2(mzo−nyo)2  物質flから逃げ出す場合も,同じ解になる。

ただし,(13)式のように2交点が存在する場合に,

どちらを採用するかは,進行方向も考慮して決 めなければならない。すなわち,z+≧z_だ から,図から明らかなように

  1→fiの場合

   もしn>0ならばZA=z_

     n〈0ならばZA=z+,

  it→1の場合

   もしn>0ならばZA=z十      n〈0ならばZA=z_

となる。これから,(11)式と同様にして,XA,

yAも求めることができる。ここではたまたま,

zについて解いたが,Xあるいはyの解でもか

まわない。

 人体の臓器は,近似的に上記のいずれかの曲 面で模擬することができる。組織(1)中に臓 器(ll)が存在し,それにX線を照射するよう な場合のシミュレーションにおいて,ここで得

られた結果は有用であった。

皿.4 実験室系から物質固定系への座標回転  散乱の理論は,入射方向をz軸にとって散乱 角を与える。これを実験室系(ターゲヅト静止 系)と呼ぶ。一般に入射光子の方向は,反応の たびに異なり,そのつど実験室系は変わる。し

xl

z

        z       ,ノノ

\ミや・95

     t

芒  、

    tN     y

x

       y Fig.5 実験室系と物質固定系の関係

周 三

たがって,光子ヒストリーを追うには,実験室 系での散乱角を,最初に設定した物質固定系で 観測した角度に変換しなければならない。

 物質固定系(x,y,z)で(θ,q),ま

たは(1。,m。ln。)方向に入射した光子がP 点において反応し,実験室系(X!,yノ, zl)で

(ψ,ρ)に散乱される様子をFig.5に示す。

問題は(ψ,ρ)を,(1。,m。,n。)を用

いて(x,y,z)系から見た方向(1,m,

n)に変換することにある。図から明らかなよ うに,この変換はまず,yノ軸のまわりにθ回転 し,次に新しいz軸のまわりにg回転すること によって得られるので,結局これらの回転の重

ね合わせになる,,

 こうして,新しい方向余弦は

 1−1。c。。ψ+⊥(1。n。si。ψ。。,ρ       む

   一m・sir.ψsinρ)

m=m。C。、ψ+⊥(m。n。Si。ψ。。,ρ        

   +1・sinψsinρ)

 n= nocos¢mso sin¢cosP

となる。ここでSo = VT=II一{である。

もしSo :0ならば

   1 == sin ¢ cos p    m = sin {b sin p    n == nocos¢

となる。

 これらの変換式は,

に用いられる。

(18)

(18)/

      サブルーチンKINEMA        KINEMAは各相互作用サブル

ーチンの中でCAI.しされる。後述M.6の方法に よって(ψ,ρ)をサンプリングして,KINE MAの入力とする。 KINEMAは新しい(1,m,

n)を計算してメインルーチンに連れ帰り,今 度はこれらが次の反応の入射方向(1。,m。,

no)になる。

皿.5 反応の型の決定

 (x。,y。,z。)から出発し,行程R進行し て(10)式で表わされる位置で反応を起こす。反応

地点(X,y,Z)が,もし物質の外部ならば 次の光子の試行に移り,内部にあれば反応の型 を選択する。どの反応が起こるかは,全断面積 に対するそれぞれの断面積の比,および乱数に

よって決まる。

(7)

光子と物質の相互作用に関するモンテカルロ・シミュレーションのプログラミング

a) b)

T= a rphoto/crtot

N

    Y

} II1 T X PHOTO

T・%驕・t戦忠。

T=T+ ortncoh/artot

     Y

}ET N

 H

PHOTO

N

    Y

}ST

IN COH

T・T・暢i,/σ島。

T = T+ orcoh/crtot

     Y} S一 T

N

N

}ET Y

 H

PAIR T・T÷啄。t。緒0

COH

PAIR N

1 S一 T Y /rM

U  PHOTO

T=T唾9h緒O

;ST N

Y r 6

 COH

  o  PAIR

Fig.6 反応の型の決定方法。 a)単原     子物質,b)分子(H20の場合)

 単原子物質の場合の選択方法をFig.6aに示 す。図において,反応の順序は順不同である。

選ばれた反応のサブルーチンへ飛んで,反応後 のエネルギーE!および(1,m,n)を連れて メインプログラムへ戻る。Eノが,追跡の終了を 意味するカヅトオフエネルギーEcより小さけ れば,光子はヒストリーを終える。もしE!が Ecより大きければ, E!に対する自由行程のサン

プリング,相互作用の選択,という具合に繰返 されることになり,これを多重過程と呼ぶ。

 物質が分子ならば,皿.3の方法によって分子 の全断面積を求め,Fig.6aと同様の判定を すべての構成元素ごとに,直列に付け加えてい けばよい。H20を例にとりFig.6bに示す。

皿.6 相互作用のサブルーチン

 各相互作用サブルーチンは,求めたい物理量 によって多少の違いはあるが,基本的にはどの

ような対象にも共通である。

1.光電効果

 ここでは,K殼の光電効果だけを考慮した。

フローチャートをFig.7に示す。 まず,入射 エネルギーEが・エッジエネルギーekより小 さい場合には,すべての光子エネルギーが電子 に与えられるとした。また,原子番号の小さい 元素に関しては・ekが小さくしたがってKX

線エネルギーも小さいので,光子エネルギー一一 E はすべて吸収されることになり,結局同じ形式 にまとめることができる。反応後の光子エネル ギーをE!,光電子が得る運動エネルギーをEe

とすれば    Ei= O

      (19)

   Ee=E

となる。

(8)

       上 原  次にE>ekのとき・起こりうる過程は・オ ージェ電子放出か,あるいはK電子とKX線の 同時放出である。いずれが起こるかは,原子の 螢光収率(OK.12)と乱数の大小によって決まる。

オージェ過程ならば(19)式と同じになる。

 一方・KX線の場合には・KαとK vOに分け Kαの強度比18)をratioとおく。

ratio == IKct /(IKa + IK,t9 ) (20)

これを乱数と比較して,ξ≦ratioならば,Kα 放出になる。KX線放出の場合のE!, Eeは

E =EKαまたはEK,S (21)

         またはE−E

 Ee =: E 一一 E

       Kfi        K ct

となる。なおKX線の角分布は等方的なので,

放出方向は②,(2)ノ式と同じである。これらE,,

(θ,g)を連れてメインプログラムに戻り,

自由行程のサンプリングからあとのルーチンを 繰り返すことになる。

E = O

Ee= E

周 三

2,インコヒーレント散乱

 以下に示す散乱公式の表記法は,文献3に拠 った。1個の静止した自由電子に対するコンプ

トン散乱の,散乱角θにおける微分断面積は,

Klein−Nishinaの式

K(α,θ)一・・z(妥)・(薯+釜+

       cos2e 一1 ) (22)

であらわされる。ここで

 a = E/ 511

 ai= ct/( 1+ ct (1 一一 cos O)) == Ei/ 511  ro = 2. 81794 1)〈 10 一  i 3 cm (23)

Y

EくeK

N

Y

; = coK

N

Y

vl・

2S ratio N

である。(22)式は,低エネルギーの極限では Thomsonの断面積((27)式)に移行する。

 実際には,軌道電子は束縛されており,原子 構造を考慮に入れる必要がある。その場合,散 乱はインコヒーレント散乱とコヒーレント散乱 に分けられる。多数の自由電子による散乱は,

非干渉性なのでインコヒーレント散乱と呼ばれ,

これの微分断面積は,原子番号Zの原子につい

 ai.,.h ( Z,(x,e)=S(Z,v)K( ct,e)

       (24)

になる。ここでVは原子への運動量移行を表わ すパラメタであり,EをkeV単位にとれば

  ・一1髪、・i・号A−i (25)

E = EKa

Ee=E EKa

E = EKB

Ee=E−EKB

s

co s e=2 S−1

Fig.7 光電効果サブルーチン

と書ける。Scattering Factor S(Z,v)は,

光子から運動量を受け取った軌道電子が,原子 を励起あるいは電離する確率をあらわす。S(Z,

v)は,前方のK(α,θ)を小さくする効果

がある。

 θを求めれば(23)式よりEノが決まり,した がって反跳電子の運動エネルギーEeは

   Ee :E 一  E  (26)

より求められる。θはコンプトン散乱に対する K。h。の方法20)によってサンプリングし,次に 佐藤・加藤16)と同様の方法により,S(Z,v)

の効果を取り入れた。Fig.8にフローチャート を示す。Hubbe 11ら9)のテーブルより,0≦v

≦50の範囲でS(Z,v)/Zデータを作成し

(9)

     光子と物質の相互作用に関するモンテカルロ・シミュレーションのプログラミング た。中間のvに対しては,直線近似の内挿によ

って求めた。高エネルギー光子についてはv>

50となるが,そのときは最大値S(Z,50)

/Z=1を代用した。20keVにおける酸素原

子についての微分断面積の,サンプリングによ る値と(24)式より求めた理論値は良く一m一致し

アこ。  ( Fig● 9 )

1,, 1,,g3

Y

   2a.1

}sS

N

y=1.2a;2

1・s4(嬉)

2a.9

N

Y

  1.2a

y=

  1・2a;2

s・5去(ナ(1一翌)2)

Y

a = a/y

cose=i− 求De

N

v. 一!L,i !:一sig:2gcose

12 .4 2

霧≦S(Z・V

z N

Y

E

cos e

    E a(1 一一 cose)

E.=

 e 1・ a(1 一 cose)

Fig.8 インコヒーレント散乱サブルーチン

3. コヒーレント散乱

 コヒーレント散乱は,電子が束縛されている とき,低エネルギーで主要な反応である。この 過程では,光子はエネルギーを失わず,束縛電

子との衝突によって散乱角のみ変化する。すべ ての電子が似たように振る舞うので,別々の電 子で散乱されても,散乱光子は干渉性を持つこ

とになる。

(10)

上 原 周 三

8

Z ト 6

cr

w 4

8 i2

o

g SAMPLING 一 THEORY

O ATOM

20 keV

    O O.!5

  cos e

サンプりングによる断面積と理 論値の比較

1

30

20   n   o

  工   m

  刀   m   z 10 H

一1 一〇.5

Fig・ 9

o

O.1 8

O.1 5

  N O.1

N隔

   0.05

o

O.Ol O.1

Fig・ IO

.一一一・一一一一一一一一・一

   1 10

   V2

Atomic Form Factorの積分分・

布関数

10 0

(11)

    光子と物質の相互作用に関するモンテカルロ  微分断面積は,Thomson散乱の断面積

   T(e)==nr2,(i十cos2 e) (27)

を用いて

   o,.h(z,at,e)=F2(z,v)T(e)

       (28)

で表わされる。

 F(Z,v)はAtomic Form Factor と呼 ばれ,軌道電子が,反跳運動量vを受けとる確 率を表わす。F2(Z,v)/Z2は一般eこ,後 方Thomson断面積を小さくする効果を持つ。

 散乱角のサンプリングはCarter・Cashwe11 の棄却法3)によった。あらかじめHubbe11ら9)

によるF(Z,v)データを用いて,積分分布

関数

A(z,v2) = JVi2(z,v)z−2 dv2

         0

      (29)

を数値計算し,v2に対するA(Z,v2)のデータ テーブルを作成する。Fig.10に酸素原子につ いてのA(Z,v2)の例を示す。 A(Z,v2)

は原子番号が大きくなるに伴ない,増加するが 大体v2〜100ぐらいで飽和するので, v2>100 に対しては一定値A(Z,100)で近似した。

 任意のEに対するv2の上限は,(25)式より

  v2 一一(E/12・4)2 (30)

   max

で与えられる。テーブルより・v ltaxに対応す るA(Z,v2 )を読みとり,これに乱数ξ

      max

を掛けてA*(Z,v2)とする。次に,今度は 逆に,A*(Z,v2)に対応するv*2を読みと

る。一方(25),(30)式より          v *2

  cos e=1 一2一=F一一一 (31)

         VmaX

であり,読みとったv*2を(31)式に代入すれば θをサンプりングできる。さらにこのθに対し,

乱数ηが

  ry fl 一li一(1+cos2e) (32)

ならば散乱角として採用し,そうでなければ最 初の乱数発生からやり直す。

・シミュレーションのプログラミング

 Fig.11 コヒーレント散乱サブルーチン

Vm,・( 2    E

12.4

)2

A(Z、V読q、)

ξ

A畳iZ.V2)=;・A(Z,V諦α)

V畳2

        V畳2

モ盾唐?ニ1顧2濡、

η

マ1圭(繭

N

cose

 θおよびE!==Eをもってメインに帰る。反応 した電子は,運動エネルギーを貰わないのでEe

=0である。フローチャートをFig.11に,20 keVにおける酸素原子についての,サンプリ ングによる微分断面積と,(28)式による理論 値との比較をFig.9に示す。

4.電子対生成

 E;≧2MeVでは,電子対生成が重要な過程 になってくる。この反応のしきい値は,陰陽電 子対の静止質量の2倍,つまり1022keVであ る。したがって,入射エネルギーのうち,1022 keVは対生成に費やされ,残りが電子対の運 動エネルギーの和になる。それが2個の電子に 等分配されると近似すれば,各電子の運動エネ

(12)

ルギーEeは

   Ee :一e (E 一一 1022 )

上 原 周 三

(33)

で表わされる。

 生成された陽電子は,物質中では,わずかな 距離進んでやがて電子と結合する。その時,電

子対の質量=1022keVが消滅し,2本の消滅

r線が生ずる。それらは共に511keVをもって,

互いに逆方向に放出されるが,角分布は等方的

である。

ID=O

 PHOTO  INCOH  COH

PAIR1

rpc・

・Eげ璽ヨ

E ・・51コ

PAIR 2

1r = r

   pc

Ee

;」  z

E = 511

e2 = el. rt 92 = Yl. rt

Fig.13 電子対生成サブルーチン

ID =1 PAIR 1

        N E〈 Ecutoff

 Y

       N  ID =1

Y

lD=2

PAIR 2

REPEAT TERMINATE

HISTORY

Fig.12 消滅γ線の識別法

 便宜上,消滅γ線に番号を付けた。第一光子 は,他の相互作用における二次光子と同様に取 扱うことができる。一方,第二光子については,

第一光子と出発位置が同じで方向は逆向きであ るから,第二光子のヒストリーを追跡するため には,対生成地点と第一一光子の放出方向を記憶 しておか;ねばならない。番号を識別するパラメ タIDを導入し,第一光子の追跡が終了すれば,

第二光子サブルーチンに飛ぶようにする。メイ ンプログラムにおけるこの処理方法をFig.12 に示す。サブルーチンは,Fig.13に示されて

いるとおり,分離してPAIR1,PAIR2と名

付けた。

皿.7 二次電子によるエネルギー付与

 各相互作用の結果,得られた二次電子の運動 エネルギーEeは,エネルギー損失や線量分布 を求めるとき,重要になる。もしEeが小さく て,物質中での山程が1mm以下程度に短かいな らば,エネルギーはイベント発生地点に局所的 に付与されると近似してよい。しかし,入射光

子が数MeVになるとEeも大きくなり,例え

ば水中における飛程は数㎝に及ぶので,電子に 対する物質の阻止能を考慮に入れたエネルギー 付与を計算し,分布を求めなければならない。

電子の物質中での振舞は,厳密にはモンテカル ロ法で追跡すべきであるが,ここでは単純に決 定論的方法によった。Eeは光電効果,インコ ヒーレント散乱および対生成において,それぞ

れ(19)・(21),(26),(33)式に求められている。

また,放出方向はGagge ro 8)、に倣い,それぞれ 以下の値を用いた。

1)光電子

 入射光子の方向となす極角をβとすれば,平 均値〆9は光子エネルギー・一 E(keV)の関数に なる。文献8に図示されているcos 19を,SSL

皿のLESQ 1によって多項式近似すれば

      う

   …β→CiEi   (34)

       i二;o となる。ここで

(13)

光子と物質の相互作用に関するモンテカルロ・シミュレーションのプログラミング  Co=8.6823E一一1 Ci=2J5933E 3

 C2=一3.2779E−6 C3=1.9971 E−9(35)

 C4=一一5.7750E−13 Cs=6.3570E−17 である。

2)インコヒーレント散乱の反跳電子

ru.6.2で求めた散乱光子の散乱角θと,α =・ E

/511を用いて

cos 0=

(1−cosθ)(1斗一α)2・

2+2a(1−cos e)+a2(1一一cos o)

(36)

3)電子対生成の陰陽電子

 入射光子の方向となす極角をβとすれば,

Bethe−Ashkinにより

   ,(9=1022/(E−1022) (37)

 一方,方位角方向については,いずれも等方 的と:考えてよい。ただし,対生成の陰陽電子対 はπだけ差がある。

 電子は直進する,と仮定すれば,光子と同様 にKINEMAを用いて,これらの方向を物質固 定系での進行方向(le・me・ne)に変換でき

x,y,z

EeS 3019

N

Y ENERGY

A CCUM.

x =x・ Ar.te y =y. Ar  Me

z=z.4r・n       e

△E癖治 r

ENERGY ACCUM

    .

N   ムE蓄Ee Y AE=E

   e

ENERGY

ACCUM.

Ee=Ee・ムE N geE Q y

RETURN

Fig. 14 二次電子によるエネルギー付与

サブルーチン

る。次に,出発点(X・,y。,Z。)から進行 方向に沿って微小な行程dr進めると,

   xe=xo +dr le

   Ye==Yo +dr me (38)

   ze==zo +dr ne

に到達する。その間のエネルギー付与dEは

   dE;(dE/d・)E・d・(39)

       となる・ここで(dE/d・)E,は・エネル ギーEeの電子に対する物質の阻止能を表わす。

また,エネルギーは・初めのEeよりdEだけ 下がって

   E各=Ee−dE      (40)

(14)

上 原 周 三

になる。

 行程をさらに進め,そのたびに(38),(39)式 を得て,エネルギー付与積算サブルーチンAC DOSへ飛ぶ。 E6駕300 keVならば・水中で の飛程は1mm以下になり,無視できる。結局 E旭≦300keVに下がるまで,同じ過程をくり 返す。最終的な行程は,最初のEeに対する飛程 に等しくなる。阻止能はPages14)らによるデー タを用いた。フローチャートをFig.14に示す。

W プログラムの具体例

N1.検出器レスポンス

 光子を検出器で測定するとき,得られる出力 を検出器レスポンスといい,入射エネルギーや 計測の幾何学的条件,検出器素材・形状などに 依存する。高純度Ge,Ge(Li),NaIの各光 子検出器のレスポンスをシミュレートした。こ れらの物質における電子の三二は短かいので,

エネルギーは局所的に付与される。エネルギー 付与をEdで表わせば,Ed=Eeとなる。

 検出器中で多重過程の起こる時間は,検出器

・増幅器系の電荷収集時間より,はるかに短か いので,入射光子1個によって生じた数回のエ ネルギー付与は,その光子が消滅するまで積算 されて,1個のパルスが形成される。計算上で は,Ec以下になるまでヒストリーを繰り返し.

Edを積算することに対応する。こうして得ら れた積算値Edを波高Edのパルスとみなして・

マルチチャネル波高分析器と同じように波高分 析を行なう。すなわち,Edをあらかじめ設定 したエネルギーチャネルに振り分け,該当する チャネルにカウント1を加える。波高分析した あと,次の光子の試行を始めるが,その都度Ed

==0にりセットしておくのが肝要である。

 Ge検出器には不感部があり,これは低エネ ルギー光子のスペクトルにかなりの影響を与え る。不感部におけるイベントは,有感部と変わ りなく起こるが,エネルギー付与は加算しては ならない。この不感部の処理は,各サブルーチ ンの中に組み込んだ。また低エネルギー光子に 対しては,入射窓による減衰を無視できない。

ベリりウムの薄い板でできた入射窓中でのヒス

トリーも厳密に追跡した。こうして全検出効率,

全エネルギー吸収ピーク効率,エスケープピー ク等のレスポンスを得ることができた。

W.2.γ線小三選線量分布

 放射線治療に用いられるr線下線源のうち,

セシウム管について線量分布をシミュレートし た。白金製の管に密封されたセシウム線源を,

水槽の中心に置いた。662keVr線の,白金

中での減衰は無視できないので,白金における

ヒストリーも皿.3に述べた手法を用いて,厳密 に追跡した。管内部で消滅してしまう光子もあ るし,散乱によってエネルギーが小さくなって 水中に放出される光子もある。

z

・dθ

!!

、ぐ,

v

1

1

1

1e

Cs TUBE

Fig.15 極座標における体積要素

 次に,ファントムの各位置におけるエネルギ ー付与を積算して,分布を求める。積算の効率 を上げるため,管軸(z軸)に関する対称性に 着目して,体積要素としてFig.15のようなz 軸のまわりの円形殼を考える。動径r.極角θ とすれば体積要素Vは

   v(r,e)==2nr2sine dr ,de (41)

となる。ここでdr,dθはチャネル分割の間

隔である。

(15)

光子と物質の相互作用に関するモンテカルロ・シミ  二次電子の水中での飛程は,最大2mmになる

が,簡単のため局所的付与とした。エネルギー 付与地点(X,y,Z)を極座標で表わすと,

   R ==  v/ x2 十 y2 十 z2

      (42)

   @ = Cos 一1 Hft

になる。R,⑰を該当する(r,θ)チャネル に振り分け,そのチャネルでのエネルギー付与 を積算する。ある体積要素における積算エネル ギーをEd(keV),字源強度を1(mCi),

試行回数をNとすれば,線量率D(rad/hour)

D( r,e) == 2. 1312・1・Ed/(t/gOgO ・N・

       v(r,e)) (43)

によって得られる。分母の100/85は,内部転 換の割合を考慮に入れた因子である。(r,θ)

を矩形分布表示に変換すれば,解析的計算値2)

との比較が容易にできる。

N. 3.電子リニアックによる高エネルギーX線 の線量分布

 線量を計算するためには,入射X線のスペク トルがわかっていなければならない。電子の制 動放射X線スペクトルは近似的に

1 ( h y )=const   { 4 ( i nv Si k) 一3 一iii itl一

     × !n (g/ )} (44)

        hン

を利用した.19)(H。itl,,一Wy。,dの式)

ここで1(hのは,X線エネルギー・hンに対す る相対強度,Elはリニアックで加速された電 子エネルギーを表わす。簡単化のため,フィ ルタはないとすれば(44)式がそのまま使える。

Elは6MeV程度であるから,スペクトルは

高エネルギー光子成分を含み,したがって二次 電子のEeも大きい。この場合には,水中にお ける電子の飛程は最大3cmになり,局所的付与 とすることはできない。

 このため皿.7の方法にしたがってb飛程に沿 った各位置におけるエネルギー付与を計算した。

外照射であるから,カーテシアン系(x,y,

   ユレーションのプログラミング

z)が便利である。水ファントムを3次元メッ シュに分割し,体積要素は立方体とした。前節 と同様,各メッシュにおける積算線量を求める。

深部線量曲線はもちろんのこと,ファントム中 に小さい金属片が放置されている場合の線量分 布の変化も見い出された。

N.4.散乱線

 散乱線のシミュレーションにおいては,物質 から逃げ出す光子に注目し,散乱線(直接線も 含む)が逃げ出す際のエネルギーEノとs散乱角 度θの2次元スペクトルを求める。エネルギー についてはIV.1と同様,Eノをエネルギーチャネ ルに振り分け,次にそのE!について散乱角分布 を求める。一次光子が偏極していなければ,z 軸に関して対称なので,9については一様であ

り,散乱方向はθだけの関数になる。散乱方向 の範囲,0。から180。までを間隔dθで分割し て角度チャネルを設定しておく。θ方向の散乱 強度は,θを該当する角度チャネルに振り分け,

カウントすることによって得られる。この操作 はFig.16に示す斜線部,つまり球面上の立体 角2πsinθdθへの入射光子を,すべてカウ ントすることに対応している。こうして2次元 スペクトルPH(E!,θ)を求めることができ

る。また

QH(E・・θ)一2懸湯)

(45)

より,単位立体角当りのスペクトルQH(E!,

θ)が得られる。

/一;;;

    e

/づ\プ

て多、

 t SCATTERER

Fig.16 散乱角θにおける立体角

一一 z

(16)

       上 原  散乱体が小さく,角度チャネルの間隔に比べ

てエスケープ地点のばらつきが無視できる場合 は,上記の方法でよい。しかし,ある特定の場 所に検出器を置いて散乱線を測定する場合,散 乱体の形状が大きいと,いろいろの方向をもっ

た散乱線を含むことになるので,θだけでは決 まらず,エスケープ地点も考慮に入れた,別の アルゴリズムが必要になる。

N5. X線写真

 X線直接撮影系をシミュレートし,X線写真 像を得るプログラムである。像を形成するフィ

ルム黒化度は照射量(R)の関数であるから,

フィルム面への照射量を求めればよい。フィル ム面を細かいメッシュに分割し,像の分解能を

上げる。

 被写体を透過した光子が,フィルム面のどの 位置に到達するかが重要である。あるメッシュ へ到達した光子のエネルギーを波高分析して,

3次元スペクトルPH(x,y,E!)を得る。

これから,照射量とスペクトルを関係付ける光 子フルエンス/Rデータ10)を用いて,メッシュ 毎の照射量が計算できる。

 試行回数に制限があるので,得られる照射量 の絶対値は,フィルムの黒化に必要な100mR 程度よりはるかに小さい。フィルムの特性に合 わせて,照射量の最大値を100mR程度に規格 化し,相対的照射量に対する黒化度を算出すれ ば,像が形成される。

V 要約

 光子と物質の相互作用について,モンテカル ロ・プログラムを開発した。重要と考えられる 基本的事項や,相互作用のサブルーチンなどの 詳細を述べ,フローチャートを図示した。これ

らを結合して,放射線物理・医学において興味 あるいくつかの例に応用し,プログラム作成と 計算を行なった。予備的な計算結果は,プログ ラムの正確さを証明している。複雑な現象や,

実験的には観測困難なデータも,計算機実験に よって,予測できるようになった。これらの他

にも,核医学,X線CT,ポジトロンCT,放

射線防護など,10MeV以下の領域の光子にま

周 三

つわるあらゆる現象に対して,ここで開発した 方法やサブルーチンは有効であり,対象に応じ てプログラムを部分的に変更するだけで適用可 能になる。

 ただひとつ,残された不備な点は,高エネルギ ー:次電子による制動放射の取扱いである。す なわち,電子が物質中を進むとき,主にイオン 化によってエネルギーを失なうが,数MeV以上 の高エネルギー領域では,制動放射によって失 なう割合が増加する。したがって,二次電子ル ーチンの中に,連続エネルギー光子の処理を追 加しなければならないことになる。この点につ いては,エネルギーと個数をサンプリングする ためのZerbyの方法20)を用いて,現在改良中で

ある。

 光子だけでなく,中性子や荷電粒子もモンテ カルロ法で取扱える対象である。これらの粒子 は,反応のチャネルが多数であり,一層複雑な ヒストリーをたどるので,プログラミングは光 子ほど簡単ではないし,要する演算時間も長く なる。しかし,マイクロドシメトりや,粒子線 治療・診断についての精密な検討に,モンテカ ルロ法は役立っであろうし,そのためのプログ ラム開発は意義がある。

 文  献

(1) Booth, T.: Lecture Note on the Mon te   Carlo Method. Univ・ of New Mexico,

  Los A旦amo.s,且980.

(2} Breitman, K. E.:Dose−Rate Tables

         137

      Cs Sources Sheathed   for Clinical

  in Platinum. British Journal of

  Radiology, 47 : 657 一 664 , 1974.

(3) Carter, L. L. and Cashwell, E・ D.:

  Particle 一Transport Simulation with   the Monte Carlo Method. Energy   Research and Development Administra一一   tion, Oak Ridge, 1977.

(4} Chan,H・ P・ and Doi,K・:The Validi−

  ty of Monte Carlo Simulation in   Studies of Scattered Radiation in   Diagnostic Radiology・ Phys. Med.

(17)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

Biol・

C hen ,

Simulation Studies of Detectors Used in the Measurement of Diagnos−

tic X−Ray Spectra. Med. Phys., 7:

627 一 635, 1980.

Ellett, W・ H. :Specific Absorbed Fractions for Photon Point Sources within a Scattering Medium. Phys・

Med. Biol・, 14 :615一 626, 1969.

富士通:FACOM OS IV/F4,TSS

コマンド文法書.1982.

Gaggero, G. :Monte Carlo Calculations for the Photofractions and Energy Loss Spectra of Ge(Li)Semiconductor Detectors. Nucl・Instr. and Meth., 94

: 481 一492, 1971.

Hubbe11, J. H. et al・ :Atomic Form Factors, lncoherent Scattering Functions, and Photon Scattering Cross Sections. 」. Phys. Chem. Ref.

Data, 4 : 471 一 538,1975・

Johns, H. E. and Cunningham,J. R.:

The Physics of Radiology. Charles C Thomas Publisher, Illinois, 1974.

Koblinger,L. and Zarand, P・:Monte Carlo Calculations on Chest X−Ray Examinations for the Determination of the Absorbed Dose and Image Quality・ Phys・ Med・ Biol., 18:518−

531, 1973.

Lederer, C. M・,Hollander, J一 M. and Perlman, 1. : Table of lsotopes.John Wi ley & Sons, New York, 1968・

中森伸行,金森仁志:モンテカルロ法によ る多重散乱線を考慮したX線スペクトルの 計算(1)水ファントムへの単色X線入射.

日本放射線技術学会雑誌,36:1−9,

1980.

Pages, L. et al.: Energy Loss,Range,

and Bremsstrahlung Yield for 10 keV to 100 MeV Electrons in Various

光子と物質の相.互作用に関するモンテカルロ・

,28 : 109 一 129 , 1983.

 C. S. et al.:Monte Carlo

(15)

(16)

(17)佐藤巳之

(18)

(19)

(20)

シミュレーションのプログラミング Elements and Chemical Compounds.

Atomic Data, 4 :1 一127 , 1972・

Raeside, D. E.:Monte Carlo Princip−

les and Applications. Phys. Med・

Biol・, 21 :181 一197, 1976.

佐藤健児,加藤二久:結合エネルギーの効

.果とコヒーレント散乱について.日本医学 放射線学会物理部会誌,2:31−47,

1982 .

       吉,上原周三:投稿予定.

Storm, E. and lsrael, H. 1.:Photon Cross Sections from 1 keV to 100 MeV for Elements Z=1 to Z=100.

Nuclear Data Tables, A 7:565 一681,

1970.

山崎文男:放射線.共立出版.1973.

Zerby, C. D.:A Monte Carlo Calcu lation of the Response of Gamma−

Ray Scintillation Counters. ln Methods in Computational Physics,

Vol.1,ed, Alder, B., Fernbach, S and

Rotenberg, M.,89−134. Academic Press, New York, 1963.

Updating...

参照

Updating...

関連した話題 :

Scan and read on 1LIB APP