5.3.18. ユーザー定義ソース¶
s-type=100のソースタイプです。usrsors.f にプログラムを書き込むことにより、ユーザー定義のソースを利用することができます。
以下に示すソースのパラメータをソースセクションに指定した場合は、ユーザーソースで定義された値よりも優先します。
以下のものが指定されたとき、ユーザー定義データがある場合もユーザー定義データより優先します。
値 |
説明 |
下限 \(x\) 座標[cm]。 |
値 |
説明 |
上限 \(x\) 座標[cm]。 |
値 |
説明 |
下限 \(y\) 座標[cm]。 |
値 |
説明 |
上限 \(y\) 座標[cm]。 |
値 |
説明 |
下限 \(z\) 座標[cm]。 |
値 |
説明 |
上限 \(z\) 座標[cm]。 |
値 |
説明 |
スピンの方向ベクトルの \(x\) 成分。 |
値 |
説明 |
スピンの方向ベクトルの \(y\) 成分。 |
値 |
説明 |
スピンの方向ベクトルの \(z\) 成分。 |
値 |
説明 |
(D=1.0) |
入射粒子の \(z\) 軸方向からの方向余弦。 all を指定した時は、等方分布。 data を指定した時は、 a-type サブセクションが必要。 |
値 |
説明 |
(D= all) |
入射粒子の方位角[degree]。 = all; 0から360の範囲でランダムに決定。 |
値 |
説明 |
(D=0.0) |
入射粒子方向の立体角範囲[degree]。 = \(-1\) ; \(\cos^2\) のbiasがかかった分布。 |
値 |
説明 |
（単色の場合）入射粒子のエネルギー[MeV/n]。 e-type = & （分布をもつ場合）入射粒子のエネルギー[MeV/n]。 |
値 |
説明 |
（初期値無し） |
（分布をもつ場合）入射粒子のエネルギー[MeV/n]。これとe0のどちらかの指定が必要。 |
値 |
説明 |
(D=1.0) |
ソース粒子のウエイト。 |
値 |
説明 |
(D=1.0) |
ソース粒子のウエイトの規格化定数。 t-type = & (D=0) 時間分布。 |
値 |
説明 |
(D=all) |
領域を限定する. |
値 |
説明 |
(D=1000) |
領域限定の際の最大再試行回数。 |
値 |
説明 |
(D=なし) |
座標変換番号もしくは座標変換定義。 |
インプットファイルで set により指定したユーザー定義定数 c1-c99 は、プログラム内で cval(1-99) により参照できます。 ただし、複数回指定した場合は最後に定義した値が有効となります。
デフォルトで usrsors.f に入っているプログラムは、以下のようなものです。
最初のコメント部に、ソースに必要な変数と、粒子指定に必要なkf codeのリストがあります。
次に、プログラム内で使用できる一様乱数と、ガウス乱数のfunctionの説明があります。
プログラムの最初の部分は、初期化の際のファイルのオープンの例があります。
最後の部分に必要な変数の値の一つの例が書かれています。
この例題を参考にユーザー定義ソースを書いてください。
ソースセクションで指定した変数は、優先しますから、その場合は、ここで定義する必要はありません。
1: ************************************************************************
2: subroutine usrsors(x,y,z,u,v,w,e,wt,time,name,kf,nc1,nc2,nc3,
3: & sx,sy,sz)
4: * sample subroutine for user defined source. *
5: * variables : *
6: * x, y, z : position of the source. *
7: * u, v, w : unit vector of the particle direction. *
8: * e : kinetic energy of particle (MeV). *
9: * wt : weight of particle. *
10: * time : initial time of particle. (ns) *
11: * name : usually = 1, for Coulmb spread. *
12: * kf : kf code of the particle. *
13: * nc1 : initial value of counter 1 *
14: * nc2 : initial value of counter 2 *
15: * nc3 : initial value of counter 3 *
16: * sx,sy,sz : spin components *
17: *----------------------------------------------------------------------*
18: * kf code table *
19: * kf-code: ityp : description *
20: * 2212 : 1 : proton *
21: * 2112 : 2 : neutron *
22: * 211 : 3 : pion (+) *
23: * 111 : 4 : pion (0) *
24: * -211 : 5 : pion (-) *
25: * -13 : 6 : muon (+) *
26: * 13 : 7 : muon (-) *
27: * 321 : 8 : kaon (+) *
28: * 311 : 9 : kaon (0) *
29: * -321 : 10 : kaon (-) *
30: * kf-code of the other transport particles *
31: * 12 : nu_e *
32: * 14 : nu_mu *
33: * 221 : eta *
34: * 331 : eta' *
35: * -311 : k0bar *
36: * -2112 : nbar *
37: * -2212 : pbar *
38: * 3122 : Lanmda0 *
39: * 3222 : Sigma+ *
40: * 3212 : Sigma0 *
41: * 3112 : Sigma- *
42: * 3322 : Xi0 *
43: * 3312 : Xi- *
44: * 3334 : Omega- *
45: *----------------------------------------------------------------------*
46: * available function for random number *
47: * unirn(dummy) : uniform random number from 0 to 1 *
48: * gaurn(dummy) : gaussian random number *
49: * for exp( - x**2 / 2 / sig**2 ) : sig = 1.0 *
50: ************************************************************************
51: implicit real*8 (a-h,o-z)
52: *-----------------------------------------------------------------------
53: parameter ( pi = 3.141592653589793d0 )
54: data ifirst / 0 /
55: save ifirst
56: character filenm*50
57: *-----------------------------------------------------------------------
58: * example of initialization
59: *-----------------------------------------------------------------------
60: if( ifirst .eq. 0 ) then
61: c filenm = 'input.dat'
62: c inquire( file = filenm, exist = exex )
63: c if( exex .eqv. .false. ) then
64: c write(*,*) 'file does not exist => ', filenm
65: c call parastop( 887 )
66: c end if
67: c open(71, file = file(i), status = 'old' )
68:
69: c close(71)
70: ifirst = 1
71: end if
72: *-----------------------------------------------------------------------
73: * example for 3 GeV proton with z-direction
74: *-----------------------------------------------------------------------
75: x = 0.0
76: y = 0.0
77: z = 0.0
78: u = 0.0
79: v = 0.0
80: w = 1.0
81: e = 3000.0
82: wt = 1.0
83: time = 0.0
84: name = 1
85: kf = 2212
86: nc1 = 0
87: nc2 = 0
88: nc3 = 0
89: sx = 0.d0
90: sy = 0.d0
91: sz = 0.d0
92: *-----------------------------------------------------------------------
93: return
94: end