2022-8-1, 2024-7-24改 IWASE
1

元となるファイル
以下の三つが必要• simple.f• simple.inp• simple.data• （＋ egs5フォルダ ⼀式）

2

フォルダ構成
おさらい• egs5/
VERSION
auxcommons/ cgview-307/ docs/
egs5run*extra_ucodes/ pegs/
run5again* slac730.pdf uccg-210721/
auxcode/
calc/
data/
egs/
include/pegscommons/ samplecodes/ tutorcodes/• egs5/calc/
egs5run*simple/• egs5/calc/simple
simple.f

simple.inp

simple.data3

simple.f
CG体系を⽤いた最⼩構成プログラム• １つの主プログラムと２つの副プログラム（サブルーチン）• program simple• subroutine ausgab• subroutine howfar

4simple.fのダウンロード

5simple.fの内容確認• nmed = 1• medarr(1) = ʻAL• chard(1) = 7.62

! number of medium 物質数
ʻ! 物質名を任意で⼊⼒（24⽂字)

物質番号1はALという名称の物質 (simple.inpで定義される)! 物質1のcharacteristic dimension• call pegs5! simple.inpを元に物質毎のテーブルを作成• iphter(i) = 0 …

! EGSオプション• iqin = 0 …

! 線源⼊⼒6simple.fの内容確認その２• ncases = 1000

! 計算回数（線源発⽣回数）• maxpict = 50

! CGVIEWプロット線源数• call shower! EGSメインルーチンのcall7

simple.inpの作成
物質情報と、物質ごとの粒⼦輸送のパラメータの設定
ELEM

物質名（任意の名称︓左詰24⽂字）
密度効果フラグ(31⽂字⽬から記⼊）&INP IRAYL=1 /END
AL

ALAL
ENER

元素記号（左詰)&INP
AE=0.521,AP=0.010,UE=2.511,UP=2.0/END
PWLF&INP /END
DECK

光⼦の下限エネルギー&INP /END 電⼦・陽電⼦の下限エネルギーに0.511 を⾜した値(MeV)

光⼦の上限エネルギー

電⼦、光⼦とも 0.01 MeV - 2 MeVまでの範囲のデータテーブル作成を意味する

電⼦・陽電⼦の上限エネルギーに0.511 を⾜した値(MeV)8

CG体系の作成（simple.dataの作成）

円柱2
円柱1
計算打ち切り領域
線源・z=0に底⾯をもつ円柱1（ターゲット）・ターゲット周囲に円柱2(周囲の計算領域)
・線源は円柱２の範囲内に置く（simple.f)
9・円柱2の外側は計算打ち切り領域

simple.dataの作成
RCC1RCC

2円筒1円筒2000005

0 0 -100 0 20

底⾯z=0

底⾯z=-10END
Z15cm⾼1010cm半径15220cm⾼ 15cm半径

1+1Z2+2Z3-2

線源

-1END
100z1, z2, z3の物質番号0は真空⼀番最後の領域は⾃動的に計算打ち切り領域に設定される
10z1z3 計算打切領域
z2

egs5runで実⾏( ifortの場合)-----------------------------------Compiling (with -save -zero -L/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/lib and O)-----------------------------------egs5job.f(6499): remark #8291: Recommended relationship between field width 'W' and the number of fractional digits 'D' in this edit descriptor is 'W>=D+7'.
102 FORMAT(I6,' NEGATIVE USTEP=',E12.6,' IR,IRNEW,IROLD=',------------------------------------^egs5job.f(22363): remark #8291: Recommended relationship between field width 'W' and the number of fractional digits 'D' in this edit descriptor is 'W>=D+7'.
1510format(2G10.5)-------------------^
egs5job.f(22355): remark #8291: Recommended relationship between field width 'W' and the number of fractional digits 'D' in this edit descriptor is 'W>=D+7'.
1480
format(2G10.5)-------------------^
egs5job.f(24575): remark #8291: Recommended relationship between field width 'W' and the number of fractional digits 'D' in this edit descriptor is 'W>=D+7'.
write(31,'((1H ,7(1PE9.3,A)))') (ELECNI(I),',',I=1,MXRAW)---------------------------^
egs5job.f(24580): remark #8291: Recommended relationship between field width 'W' and the number of fractional digits 'D' in this edit descriptor is 'W>=D+7'.
write(31,'((1H ,5(1PE11.5,A)))') (CAPIN(I),',',I=1,MXRAW)---------------------------^
egs5job.f(24583): remark #8291: Recommended relationship between field width 'W' and the number of fractional digits 'D' in this edit descriptor is 'W>=D+7'.
write(31,'((1H ,7(1PE9.3,A)))') (SCPROF(j,i),',',j=1,31)-----------------------------^----------------------------------------------Does this user code read from the terminal?
(Enter 1 for yes, anything else for no)----------------------------------------------1******************************************************** User code simple/simple.f has been compiled and is starting ********************************************************
Ready for user input:
realusersys========================
egs5run script has ended========================

0m1.288s
0m0.737s
0m0.076s

11

CGVIEWで⾶跡の確認

egs5/calc/egs5job.picを開きます12円柱1にz=-5から光⼦が50本真直に⼊射し5cmのアルミ内で半数程度が散乱している(散乱せず透過した光⼦は軸上を進むのみ）

13

エクセサイズsimple.fにてmaxpict を50から⼤きめの数字にしてみる→ 実⾏してCGVIEWで⾶跡を観察

14

15

maxpictは多すぎるとCGVIEWの動作が重くなるので注意エクセサイズ2
simple.fにて線源種類と位置を変更してみる→ 実⾏してCGVIEWで⾶跡を観察
16

17・電⼦の⼤半はアルミの表⾯で⽌まっている・⼀部は後⽅へ散乱している
・制動放射線（ガンマ線）を若⼲放出している

ベータ線をアルミで遮へい
アルミを薄く、線源位置を戻してみる

1ミリのアルミに1.253MeVのベータ線を⼊射

18

19・電⼦の⼤半はアルミを透過している
・⼀部は後⽅へ散乱している
・制動放射線（ガンマ線）を若⼲放出している

ターゲット背後から放出される放射線を書き出す
CG版tutor121

線源

z3 計算打切領域

z1z2

20

CG版tutor1
tutor1/tutor1.fのausgabの⼀部をsimple/simple.fのausgabへ貼り付けて実⾏

追記

追記

21

iarg=3かつir=3の条件については後で考察

運動エネルギー、粒⼦種、z軸からの⾓度（０のときz軸上)

22

枠線内のデータでヒストグラムを作ってみるgnuplotを⽤いたヒストグラム作成
set style fill solid
filter(x,y)=floor(x/y)*yplot ʻout' u (filter($1,0.02)):(1)smooth frequency with boxes________
↑
ビン幅

23

IARGについて
その他詳細

移動エネルギーカット1
エネルギーカット2

体系外など

エネルギー付与
光電効果の後処理としてコンプトン散乱の後処理としてK-shell EIIの後処理として

24

iarg=3 かつ ir=3とは

2
1

線源

z3 計算打切領域

z1z225

ターゲット背後から放出される放射線を書き出すには

2
1

線源

z3 計算打切領域

z1z2

irold = 1 かつ ir = 2 かつ z = アルミ厚のとき
26surface crossing設定例
領域1から領域2へ移動する粒子

z=0.1d0と倍精度での指定がよい
z=0.1と単精度の指定は条件を満たさない恐れ

27

28

電⼦と光⼦を分けて出⼒する

29

結果
fort.92 (electron)
0.540271512664514
0.694575869707918
0.701306645181710
0.829987733342612
0.203134256948940
0.807720792774279
0.579118950986061
0.886655119575416
0.423448134695322
0.829891512015685
0.615723358845304
0.781277070179265
0.855219289021240
0.846917379926870
0.813741911554502
0.185829413707758
0.691270549925618:

fort.91 (photon)
0.107140484629189
0.102267272684809
0.395397252626522
0.639925780212798
0.247087137482216
0.267966328782001
6.659264538373698E-002
3.107949150083373E-002
7.063016926108848E-002
0.137363839686406
0.120023201528415
9.653519355772798E-002
0.237160049495556
4.857712306131567E-002
4.395142113070239E-002
0.134551177997255
2.291256156484392E-002
:30

31

シンチレータの応答、標的の熱計算CG体系でtutor2
energy deposit量とその分布の計算

z3計算打切領域

z2z1線源↑
領域1内に吸収されたエネルギーを求める32

EGS5の吸収エネルギー
領域

電子1
入射光子

e1

e1’’
e1’

電⼦1が領域中に失うエネルギー= (e1-e1ʼ) + (e1ʼ-e1ʼʼ) + …

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
コンプトン散乱 ~~~~~~~~~~~~~~~~~~~~光電吸収e~2 e2’電⼦2電⼦2が領域中に失うエネルギー= (e2-e2ʼ) + (e2ʼ-e2ʼʼ) + …
33

１ステップとdE
計算の段階ごとのedep

e1

e1’’
e1’

3

edep=e1’-e1’’~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 2 edep=e -e ’~~~~~~~~~~~~~~~~
計算1~~~~~~~~
e~2 e2’
edep=01

1

電子2
１⼊射放射線あたりの吸収エネルギー = edep を⾜し合わせたもの
34

35

fort.93の中身
1.25300002098084
1.25300002098083
0.506240130558479
0.712728508316322
0.558424151272919
0.847693228025980
1.25300002098083
1.25300002098083
0.357737686830670
0.551693375799125
1.25300002098084
0.423012287638223
1.25300002098083
1.04986576403189
1.25300002098083
0.376270297827966
0.445279228206556
0.673881069994773
0.366344901405419
0.336851243980807
0.829551886285513
0.423108508965151
0.340862623567150
0.637276662135530

36エネルギー収⽀を確認する(tutor2)•

線源エネルギー (1) = 吸収エネルギー (2) + 体系外へ逃げたエネルギー (3)

(1)
1 MeV

2
1photon
線源

(2)
吸収エネ
ルギー
•

エネルギー収⽀の取れていない計算はどこかにミスがある可能性

37

3(3)
体系外へ逃げたエネル
ギーメインルーチンの変更点

38

ausgabの変更点

39

結果

(1) 線源エネルギー(1 MeV x 10,000回)(2) 領域1の吸収エネルギー（約34%, 厚みによる）(3) 体系外へ逃げたエネルギー(2) + (3) = (1) となった
⇒ エネルギー収⽀が確認できた
40

測定器の応答計算(tutor3)• 3インチのNaIの応答を計算する（簡易版）• 結晶に直接光⼦を⼊射

41

simple.inp (PEGS⼊⼒:物質とデータテーブル作成）
COMP&INP NE=2,RHO=3.67, PZ=1,1 IRAYL=1 /END
NAI
NAI
NA I
ENER
&INP AE=0.521,AP=0.01,UE=2.511,UP=2.0 /ENDPWLF
&INP /ENDDECK&INP /END

EGS研究会ホームページ
講習会資料 8-7 materials 参照
https://rcwww.kek.jp/egsconf/2015-nm/materials.html

42

simple.data (計算体系）
RCC
RCC
END
Z1Z2Z3END
1 012+1+2-20 0 0 0 0 7.62 3.810 0 -10 0 0 20 15-1043メインルーチン4849
50
51
52
53
54
55
56

nmed=1if(nmed.gt.MXMED) thenwrite(*,*)'error:nmed exceeds MXMED'
stop
end ifcall block_set
medarr(1)='NAI! Initialize some general variables'102 ! source
103
iqin = 0! Incident particle charge - photons
104
ekein = 0.662d0! Incident particle kinetic energy
105
xin = 0.0! Source position106
yin = 0.0107
zin = -5.0108
uin = 0.0! Moving along z axis
109
vin = 0.0110
win = 1.0111
irin = 0! Starting region (0: Automatic search in CG)
112
wtin = 1.0! Weight = 1 since no variance reduction
used
113

44

185
186
187
188
189
190
191
192
193
194
195
cases
196
197
198
199

uf(1)=0.0
vf(1)=0.0
wf(1)=0.0 ! Needed if lpolar(i)=1
de = 0d0
call shower (iqin,etot,xin,yin,zin,uin,vin,win,irin,wtin)
if(de.ne.0) write(93,'(1pe13.5)') dencount = ncount + 1! Count total number of actual

if(iwatch.gt.0) call swatch(-1,iwatch)end do! End of CALL SHOWER loop

ausgab
48
49
50
51
52
53
54
55
56

nmed=1if(nmed.gt.MXMED) thenwrite(*,*)'error:nmed exceeds MXMED'stop
end ifcall block_set
medarr(1)='NAI! Initialize some general variables'102 ! source
103
iqin = 0! Incident particle charge - photons
104
ekein = 0.662d0! Incident particle kinetic energy
105
xin = 0.0
! Source position106
yin = 0.0107
zin = -5.0
108
uin = 0.0
! Moving along z axis
109
vin = 0.0
110
win = 1.0
111
irin = 0! Starting region (0: Automatic search in CG)
112
wtin = 1.0! Weight = 1 since no variance reductionused113

45185
186
187
188
189
190
191
192
193
194
195
cases
196
197
198
199

uf(1)=0.0
vf(1)=0.0
wf(1)=0.0 ! Needed if lpolar(i)=1
de = 0d0
call shower (iqin,etot,xin,yin,zin,uin,vin,win,irin,wtin)
if(de.ne.0) write(93,'(1pe13.5)') dencount = ncount + 1! Count total number of actual

if(iwatch.gt.0) call swatch(-1,iwatch)end do! End of CALL SHOWER loop

46

３インチNaIの応答, Eγ=622keV
100000

Events

10000

1000

100

10

100.1

0.2

0.3
0.4E (MeV)
47

0.5

0.6

0.7

• 2024-7-24改 ページ数追加、p38と39のcmn1/からmaxpictを削除

48