宇宙線を空から降らせる
参考
phits/lecture/advanced/cosmicray
マニュアル 5.3.19.9 章
半径200cmの球殻から大気中の宇宙線を2000発降らせている。スペクトルの設定は以下。
緯度(glat) = 33 (度)
経度(glong) = 130 (度)
高度(alti) = 0 (km)
これは佐賀市の場合
例としてミューオンと中性子だけ降らせている。
cell=21が地上。22が地下。22は物質番号 -1 を割り当てて下からの成分を完全にカットした。
鉛ガラス検出器（cell=1）を地上に置いた。
ミューオン (2d_xz_muon.eps)
中性子 (2d_xz_neutron.eps)
energy.eps
地上(cell=21)における粒子のエネルギー分布。
angle.eps
地上(cell=21)から地下(22)に入る時の粒子の角度分布。中性子のほうがミューオンよりも角度が大きい。
deposit.eps
鉛ガラス(cell=1)のエネルギー付与分布。
統計が少ない。解決方法は例えば
maxcasを増やして計算回数を上げる宇宙線が発生する球殻の半径
c1を小さくする鉛ガラスを複数置く
（参考：統計を増やした deposit.eps）
into_detector_dmp.out
検出器(cell=1)に入った粒子の情報をテキストファイル（
into_detector_dmp.out）に書き出している。
-1.300000000000000D+01 5.999999999683743D+00 3.167490567208137D+00 1.400311182461877D+01 -3.162570913876687D-01 -4.350904000259421D-01 -8.430170792767353D-01 5.043608622953946D+02 1.000000000000000D+00 6.261725450254336D+00 7.440000000000000D+02 1.000000000000000D+00
-1.300000000000000D+01 4.270813381084424D+00 5.999999999166689D+00 2.737788156725464D+01 4.870208692571157D-01 -8.333107648138632D-01 -2.615412819295985D-01 1.152322691655301D+04 1.000000000000000D+00 6.276482146981324D+00 8.670000000000000D+02 1.000000000000000D+00
-1.300000000000000D+01 -5.999999999571961D+00 -3.746179009237953D+00 2.902283900132777D+01 4.280392466046081D-01 -5.582329588916357D-02 -9.020344577687874D-01 1.024044951552910D+04 1.000000000000000D+00 5.713130883858525D+00 1.014000000000000D+03 1.000000000000000D+00
-1.300000000000000D+01 1.363282460846900D+00 5.999999999579546D+00 2.125876286168722D+01 -3.618266020839387D-01 -4.204536988867433D-01 -8.320457902765009D-01 7.622341503246878D+03 1.000000000000000D+00 5.976120951789371D+00 1.041000000000000D+03 1.000000000000000D+00
-1.300000000000000D+01 2.324424479812270D+00 5.999999999332699D+00 6.830101307872464D+00 -2.182880564808399D-01 -6.673006902635852D-01 -7.120815354800041D-01 3.581966376533930D+04 1.000000000000000D+00 6.358585399034658D+00 1.170000000000000D+03 1.000000000000000D+00
2.112000000000000D+03 2.987449763408414D+00 -5.999999999100431D+00 1.345475335831164D+01 -4.317182097982893D-01 8.995691624174087D-01 -6.629260408372918D-02 3.916245759037464D+00 1.000000000000000D+00 7.036389105658317D+01 1.272000000000000D+03 1.000000000000000D+00
2.112000000000000D+03 2.605150703877696D+00 5.999999999260335D+00 1.536708328395350D+01 4.251299510657168D-01 -7.396648733591683D-01 -5.216899460651190D-01 8.222506588419546D-01 1.000000000000000D+00 1.502141858288128D+02 1.449000000000000D+03 1.000000000000000D+00
-1.300000000000000D+01 5.999999999524234D+00 -3.724896567140975D+00 1.773088121318878D+00 -4.757660453435825D-01 4.631870654161487D-01 -7.477328483685229D-01 4.150098315192381D+03 1.000000000000000D+00 6.474771844197241D+00 1.475000000000000D+03 1.000000000000000D+00
1.300000000000000D+01 5.999999999386787D+00 1.074583986118774D+00 1.932136694372502D+01 -6.132135948588188D-01 -4.911351716966157D-01 -6.186722316403577D-01 5.864092166797841D+02 1.000000000000000D+00 6.192399315136884D+00 1.859000000000000D+03 1.000000000000000D+00
この場合、各行に 12 列の情報が書かれている。どの情報を書き出すかは
main.inpの以下で指定されている。
[ t-cross ]
file = into_detector.out
part = all
mesh = reg
reg = 1 $ number of crossing regions
non r-from r-to area
1 21 1 1.0
e-type = 2 # e-mesh is log given by emin, emax and ne
emin = 0 # minimum value of e-mesh points
emax = 1e+10 # maximum value of e-mesh points
ne = 1 # number of e-mesh points
unit = 1 # 1:[1/cm^2/source]
axis = eng
output = current # surface crossing flux
epsout = 0
dump = -12 # (D=0) number of dumped data, <0: ascii, >0: binary
1 2 3 4 5 6 7 8 9 10 18 19
dumpデータに関しては例えば以下を参照。
phits/lecture/advanced/sourceB/phits-lec-sourceB-jp.pptx
main.inp
$======================================
[ Parameters ]
icntl = 0 # (D=0) 3:ECH 5:NOR 6:SRC 7,8:GSH 11:DSH 12:DUMP
maxcas = 2000 # (D=10) number of particles per one batch
maxbch = 1 # (D=10) number of batches
negs = 2 # (D=-1) =-1:original, =0:No, =1,2:EGS
e-mode = 2 # (D=0) 0: Normal, 1: Event generator mode Ver.1, 2: Ver.2
irqmd = 1 # (D=0) 0: JQMD legacy version, 1: JQMD-2.0
mdbatima = 3000 # (D=500) max database size of ATIMA
maxbnk = 100000 # (D=10000) maximum bank memory length
iMeVperU = 1 # (D=0) 0:[MeV] or 1:[MeV/u] is unit of tally output
esmax = 100000000. # (D=1.0e6) max energy/u for charged particle range
$======================================
[ source ]
set:c1[200.0] $ Source radius [cm]
totfact = pi*c1**2 # (D=1.0) global factor
infl:{source.inp}
proj = muon+
infl:{source.inp}
proj = muon-
infl:{source.inp}
proj = neutron
$======================================
[ surface ]
1 rpp -6 6 -6 6 0 31.5
2 so c1+0.1
3 pz -1
$======================================
[ material ]
$ 鉛ガラス（密度 = 5.2g/cm3）
$ [Ref] PERFORMANCE OF THE VENUS LEAD-GLASS CALORIMETER AT TRISTAN
$ Nuclear Instruments and Methods in Physics Research A271 (1988) 432-441
mat[1] K 0.00927418 Na 0.00939653 O 0.605591 Pb 0.154165 Sb 0.0010569 Si 0.220517
$======================================
[ cell ]
1 1 -5.2 -1 $ 鉛ガラス検出器
21 0 -2 3 #1 $ 地上
22 -1 -2 -3 $ 地下
99 -1 2
$ ======================================
[ t-track ]
infl:{2d_xz.inp}
file = 2d_xz_muon.out
part = (muon+ muon-)
t-type = 2
tmin = 0
tmax = 10
nt = 20
$ ======================================
[ t-track ]
infl:{2d_xz.inp}
file = 2d_xz_neutron.out
part = neutron
t-type = 3
tmin = 1
tmax = 100000
nt = 20
$======================================
[ t-track ]
file = energy.out
part = muon+ muon- neutron
mesh = reg # mesh type is region-wise
reg = 21
e-type = 3 # e-mesh is log given by emin, emax and ne
emin = 1 # minimum value of e-mesh points
emax = 1.0e+6 # maximum value of e-mesh points
ne = 60 # number of e-mesh points
unit = 1 # unit is [1/cm^2/source]
axis = eng # axis of output
gshow = 1 # 0: no 1:bnd, 2:bnd+mat, 3:bnd+reg 4:bnd+lat
epsout = 1 # (D=0) generate eps file by ANGEL
y-txt = Flux [1/cm^2/sec]
$======================================
[ t-cross ]
file = angle.out # file name of output for the above axis
part = muon+ muon- neutron
mesh = reg
reg = 1 $ number of crossing regions
non r-from r-to area
1 21 22 1.0
e-type = 2 # e-mesh is log given by emin, emax and ne
emin = 0.0 # minimum value of e-mesh points
emax = 1.0e10 # maximum value of e-mesh points
ne = 1 # number of e-mesh points
a-type = -2 # e-mesh is log given by emin, emax and ne
amin = 0.0 # minimum value of e-mesh points
amax = 180.0 # maximum value of e-mesh points
na = 20 # number of e-mesh points
unit = 4 # unit is [1/cm^2/sr/source]
axis = the # theta = angle
output = a-curr # surface crossing flux
epsout = 2 # (D=0) generate eps file by ANGEL
y-txt = Flux [1/cm^2/sr)/sec]
$======================================
[ t-deposit ]
file = deposit.out # file name of output for the above axis
part = all muon+ muon-
mesh = reg # mesh type is region-wise
reg = 1
e-type = 3
emin = 1
emax = 1000
ne = 60
unit = 3 # 3: Number [1/source], 4: Number [1/nsec/source]
output = deposit # total deposit energy
axis = eng
gshow = 1
epsout = 1
$======================================
[ t-cross ]
file = into_detector.out
part = all
mesh = reg
reg = 1 $ number of crossing regions
non r-from r-to area
1 21 1 1.0
e-type = 2 # e-mesh is log given by emin, emax and ne
emin = 0 # minimum value of e-mesh points
emax = 1e+10 # maximum value of e-mesh points
ne = 1 # number of e-mesh points
unit = 1 # 1:[1/cm^2/source]
axis = eng
output = current # surface crossing flux
epsout = 0
dump = -12 # (D=0) number of dumped data, <0: ascii, >0: binary
1 2 3 4 5 6 7 8 9 10 18 19
$ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
$ kf x y z u v w e wt time c1 c2 c3 sx sy sz name nocas nobch no
source.inp
<source> = 1.0 # should be 1.0 for e-type=25
s-type = 9 # mono-energetic axial source
r1 = c1 # radius of sphere [cm]
r2 = c1 # radius of source circle [cm]
dir = iso # inner direction with uniform dis. by analog
e-type = 25 # Cosmic-ray source
icenv = 1 # 0>: Terrestiral GCR, D=0: Free-space GCR, <0:Free-space SEP
solarmod = 0 # Solar activity W-index
glat = 33
glong = 130
alti = 0 $ km
2d_xz.inp
$ ======================================
$ [ t-track ]
$ file = 2d_xz_muon.out
$ part = (muon- muon+)
mesh = xyz
x-type = 2
nx = 200
xmin = -c1*1.1
xmax = c1*1.1
y-type = 2
ny = 1
ymin = -c1*1.1
ymax = c1*1.1
z-type = 2
nz = 200
zmin = -c1*1.1
zmax = c1*1.1
e-type = 2
ne = 1
emin = 0
emax = 1e+10
unit = 1 $ 1: [1/cm^2/source]
axis = xz
epsout = 1
gshow = 3 # 0: no 1:bnd, 2:bnd+mat, 3:bnd+reg 4:bnd+lat
z-txt = Flux [/cm2/sec]
$ Rotating figure by 90 degree
trcl = 0 0 0 0 0 1 0 1 0 1 0 0
x-txt = x [cm]
y-txt = z [cm]