Geant4を用いた電磁カロリメータの性能評価

Geant4 を用いた
ECAL シミュレーション
信州大学理学部物理科学科
高エネルギー物理学研究室
学籍番号 11S2025C
氏名伊佐見将
指導教官
竹下徹
長谷川庸司
1
概要
高エネルギー物理学では、高エネルギーにまで上げられた粒子同士の衝
突によって未知の粒子を見つけ出そうとしている。その実験に用いられる
測定器の一つにカロリメータがある。カロリメータは二種類に分類され、
光子や電子のエネルギーを測定する電磁カロリメータ（Electromagnetic
Calorimeter、略称 ECAL）とハドロンのエネルギーを測定するハドロン
カロリメータ（Hadron Calorimeter、略称 HCAL）がある。実際に作成
して実験したいが、時間も費用も掛かってしまうので、シミュレーション
を行うことで測定に最適なカロリメータを見つける必要がある。今回は
ECAL について実験モデルを考え、粒子を入射させた場合にどのような
挙動をするのかを Geant4 というシミュレーションソフトを用いて研究し
た。
前半では、粒子が入射する角度によって位置分解能（エネルギー重み半径
とその RMS）とエネルギー分解能を鉄 (Fe)、鉛 (Pb)、タングステン (W）
の三種類の物質で確認し、タングステンが最適であることが分かった。
後半では、まず、入射位置を変更した場合や方位角を変更した場合に異
なった結果が得られるかを確認した。入射位置や方位角による影響はあ
まり見られなかった。最後に実用型モデルについて位置分解能、エネル
ギー分解能に加え、モリエール半径を求めることで、最適な物質を探索し
た。結果として、位置分解能やモリエール半径が良いものを作製したい
のであればタングステンを用い、エネルギー分解能が良いものを作製し
たいならば鉄を用いた ECAL を作製すべきであるということがわかった。
2
目次
1
序章
1.1 目的 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 高エネルギー実験用カロリメータ . . . . . . . . . . . . . .
2
電子と物質との相互作用
2.1 電磁シャワー . . . .
2.2 エネルギー損失現象
2.2.1 荷電粒子 . . .
2.2.2 制動放射 . . .
2.2.3 対生成 . . . .
2.3 モリエール半径 . . .
3
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
解析方法・結果
3.1 Geant4 . . . . . . . . . . . . . . . .
3.2 シミュレーションモデル . . . . . .
3.3 角度毎のシミュレーション . . . . .
3.3.1 位置分解能 . . . . . . . . .
3.3.2 エネルギー分解能 . . . . . .
3.4 入射位置を変更した場合 . . . . . .
3.5 方位角を変更した場合 . . . . . . .
3.6 実用型モデルでのシミュレーション
3.6.1 位置分解能 . . . . . . . . .
3.6.2 エネルギー分解能 . . . . . .
3.6.3 モリエール半径 . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
5
5
6
.
.
.
.
.
.
7
7
9
9
10
11
11
.
.
.
.
.
.
.
.
.
.
.
13
13
14
15
15
17
22
24
26
26
27
28
4
まとめ
30
5
今後の展望
31
謝辞
32
参考文献
33
補遺
34
3
プログラム
41
4
序章
1
1.1
目的
本研究の目的は、電磁カロリメータの性能を評価し、吸収層に用いる
物質として最適なものを提案することである。そのために以下の順でシ
ミュレーションを行う。
1. 角度がついて入射した場合（入射角を変更した場合）のシミュレー
ションを行い、位置分解能、エネルギー分解能を算出し、比較する。
2. 入射位置や方位角を変更した場合に得られるデータに違いがみられ
るかどうか検証する。
3. 実際に使用が期待されるモデルを用いてシミュレーションを行い、
位置分解能、エネルギー分解能、モリエール半径のデータを取得し、
総合してどの物質が最適であるかを結論付ける。
5
1.2
高エネルギー実験用カロリメータ
高エネルギー実験では、加速器を用いて粒子を高エネルギーまで加速
させ、測定器の中心で粒子同士を衝突させている。実際の実験装置の例
として CERN にある ATLAS 実験の測定器（図 1）を載せる。
図 1: ATLAS 実験で用いられている測定器
カロリメータは入射粒子のエネルギーを測定する役割を担っており、２
つのカロリメータの違いは次のとおりである。
● 電磁カロリメータ（ECAL）
電子、陽電子、光子などのエネルギーを測定する粒子検出器。電磁
シャワーを発生させて粒子のエネルギーを内部で検出する。
● ハドロンカロリメータ（HCAL）
ハドロンのエネルギーを測定する粒子検出器。高エネルギーのハド
ロンが入射したときにハドロンシャワーを発生させて内部でエネル
ギーを検出する。
カロリメータは吸収層と検出層の二層を組み合わせたレイヤーを何層に
も重ねて作成されている。吸収層には物質量の大きい物質が使用されて
いる。検出層はタイル型やストリップ型など様々な形状のものが提案さ
れている。カロリメータは吸収層に物質量の大きな物質を用いているた
めに測定装置の大半の重量を占めているため、可能な限りコンパクトに
してコストを抑えることが求められている。
6
2
電子と物質との相互作用
本研究では、電子を電磁カロリメータに入射させてシミュレーション
を行った。
2.1
電磁シャワー
高エネルギーの電子が物質中に入射すると、物質の原子核と相互作用
して制動放射を起こし、高エネルギーの γ 線を放出する。その γ 線が原
子核と相互作用を起こすと、電子‐陽電子対を生成する。この電子と陽
電子も制動放射を起こす。. この反応が何度も繰り返され、粒子の広がり
がシャワーのように見えることから電磁シャワーと呼ばれる。その様子
を図 2 に示す。
図 2: 電磁シャワー
このシャワーが起きることによって、ECAL 内部にエネルギーが落とさ
れるのである。ここで述べた制動放射や対生成については、2.2.2 と 2.2.3
で説明する。
次に、シャワーの発達について考える。シャワーの発達を表すための有
7
用な長さの単位に放射長 X0 がある。放射長は次頁の式 (1) のように表さ
れる。
X0 =
7164A
√
[mm]
ρZ(Z + 1) ln (287/ Z)
(1)
Z は物質の原子番号、A は質量数、ρ は密度である。電磁シャワーの過
程を単純化すると、電子や陽電子は 1X0 進むごとに γ 線を放出し、γ 線
は 97 X0 ごとに対生成していく。そのため放射長は、電磁シャワーの縦方
向の広がり（深さ）に大きく関連している。コンパクトな電磁カロリメー
タを目指すには、式 (1) からわかるように Z が大きくなければならない。
また、電磁シャワーの終端について考える。電子や陽電子の制動放射や γ
線の対生成には、閾値となるエネルギー値がある。このエネルギーは臨
界エネルギー Ec と呼ばれ、Ec を下回ると変換が起こらなくなる。Ec は
下の式 (2) で表される。
610
[M eV ]
(2)
Z + 1.2
Z が大きいほど Ec は小さくなることがわかる。このことも電磁カロ
リメータの素材として Z の大きな物質が最適であることの根拠となって
いる。
Ec =
8
2.2
エネルギー損失現象
電子が物質に入射したときには、次の３つのエネルギー損失現象が起
きると考えられる。
2.2.1
荷電粒子
電子は負の電荷をもつ荷電粒子である。荷電粒子が物質中を進むと、主
に物質内の電子と反応する。これは原子核との散乱と電子との散乱の違
いによって起こる。前者は弾性散乱となるため、持っているエネルギー
の一部しか移行しない。しかし、後者の場合は非弾性散乱となるためエ
ネルギー移行が大きくなる。したがって、荷電粒子のエネルギー損失の
大部分は電子との反応によるものであると言える。十分大きな運動エネ
ルギーを持っていれば、軌道に沿って原子を電離するか、励起してエネ
ルギーを与える。励起された原子や分子は基底状態に戻るために γ 線を
放出する。
原子が電離するときには、物質の電離特性を表す変数である阻止能 S(T )
に従う。S(T ) は運動エネルギー T の荷電粒子が軌跡の単位長さ当たりに
失う運動エネルギーとして式 (3 ) として定義されている。
dT
= nion I¯
(3)
dx
nion は単位長さ当たりに荷電粒子が作る電子・イオン対の数、I¯ はその
平均の電離エネルギーである。I¯ は原子番号が大きいときには 10Z と近
似できる。式 (3) の負の符号は、粒子のエネルギー減少を表しており、運
動エネルギーの変化 dT が負であることを示している。
阻止能は電磁相互作用しか含まないので高い信頼度で計算することが可
能である。ハンス・ベーテとフェリックス・ブロッホは相対論的粒子に対
して式 (4) を導いた。
S(T ) = −
S(T ) =
4πQ2 e2 nZ
2mc2 γ 2 β 2
[ln(
) − β2
mβ 2 c2
I
(4)
ここで m は電子の静止質量、β は粒子の速度 v と真空中の光速 c の比
1
はローレンツ因子 (1 − β 2 )− 2 、Q は粒子の電荷 ze、Z は物質の原
子番号、n は単位体積当たりの原子数である。高エネルギーの粒子や β 崩
壊で放出する電子に対しては相対論的補正が大きいので、この式を用い
て阻止能を算出する。また、電離損失は式 (4) の β −2 因子のために粒子
v
、γ
c
9
の速度が小さいときはその運動エネルギーに大きく左右される。しかし、
今回の場合は高エネルギーであるため、β ≈ 1 となり、それほど大きな
影響はないと考えられる。
2.2.2
制動放射
電子は、原子中の電子や核の作る強い電場により大きな加速を受ける。
これによって電磁波の放射が起こる。この現象を制動放射 (bremsstrahlung)
と呼ぶ。したがって、物質を通過する電子の全エネルギー損失は、前項
の荷電粒子と合わせて模式的に次のように表せる。
(
dT
−
dx
)
(
=
tot
dT
−
dx
)
(
+
ion
dT
−
dx
)
(5)
brem
高エネルギー電子の制動放射と電離損失の比は、近似的に次の形で表
される。
(
)
dT
dx
( )brem
dT
dx ion
≈
TZ
1200mc2
(6)
ここで、Z は物質の原子番号、m は入射電子の静止質量、T はその運
動エネルギーである。高エネルギーでは、電離損失は密度効果のため飽
和して一定となり、
Z
(M eV /(g/cm2 ))
(7)
A
と近似されるため、式 (5) において制動放射の項が大部分を占めている。
式 (6) から高エネルギーでのエネルギー損失は電子のエネルギーに比例し
ていることがわかる。そこで、放射長を導入する。放射長は、運動エネ
ルギーが元の 1e になるまでに電子の進む距離として定義される。式 (5) と
式 (6) を用いて次のように示される。
Smin ≈ 3.5
(
dT
dx
)
= −
brem
T
A
、ここで X0 ≈ 170 2 (g/cm2 )
X0
Z
(8)
上の式を T0 から T まで積分すると、次のように T が得られる。
T = T0 e
10
− Xx
0
(9)
式 (9) から、高エネルギー電子は物質中での制動放射により、指数関数
的に運動エネルギーを失っていくことがわかる。
2.2.3
対生成
光子 (γ 線) が十分なエネルギーを持つと、物質に吸収され、正負の電
荷をもった対が作られる。本研究では、光子から電子―陽電子対が生成
される。
γ → e− + e+
(10)
このような変換はいかなる保存則を破らないときのみ可能である。しか
し、単一の光子は質量を持たないため、エネルギー・運動量の保存則が成立
せず、有限質量の粒子の対が作られない。このことを理解するために、光子
が非常に小さな（電子よりもはるかに小さな）質量を持つと仮定する。光子
の静止系で見ると系のエネルギーはほぼ 0 に近い。一方、終状態のエネル
ギーは最低でも 2 粒子の静止質量の和となるので、始状態よりも必ず大きく
なる。したがって、式 10) のような変換が起こるのは、物質中に限られるこ
とがわかる。物質中では、反跳原子核がエネルギー・運動量の保存則を満た
すようにエネルギーを吸収する。電子―陽電子対生成のエネルギーの閾値
は、どちらも同じ質量なので、hν ≈ 2mc2 = 2×0.511MeV ≈ 1.022MeV
である。対生成の断面積は、物質の原子番号を Z とすると Z 2 に比例し、
閾値から急激に立ち上がる。本研究のような 100MeV を超える高エネル
ギーでは、電子―陽電子対生成の断面積は飽和に達し、一定の平均自由
行程で与えられる。その値は、式 (11) のように物質中の電子の放射長に
ほぼ等しくなる。
Xpair = (µpair )−1 ≈
2.3
9
X0
7
(11)
モリエール半径
モリエール半径 ρM は、シャワーの横方向の広がりを表す。電磁カロリ
メータ内部に入射された全エネルギーのうち 90 ％が含まれる半径のこ
とで、その値は放射長と臨界エネルギーを用いて式 (12) のように定義さ
れる。
11

ρM
√

Es 
4π
= X0 ·
Es = me c2
≈ 21M eV 
Ec
α
(12)
この値が小さいほど、良いカロリメータであると言える。その理由と
して以下の図 3 を用いる。
図 3: モリエール半径の大きさによる違い
図 3 のように隣の入射点までの距離が同じ場合、シャワーの半径が大
きくなると重なりができてしまう。そうなるとエネルギー検出したとし
てもどちらのシャワーによるものなのか区別がつかなくなるため、モリ
エール半径は小さい方がよい。
12
3
3.1
解析方法・結果
Geant4
Geant4 とは、粒子が物質に入射した際に生じる複雑な振る舞いや反応
（制動放射や対生成等）をコンピュータ上で再現することができる、モン
テカルロシミュレーションのソフトウェアのことである。ソースファイ
ルを自ら書き換えることで、カロリメータの形状や使用する物質、入射
粒子のパラメータ等も設定することができるので、様々な条件でのシミュ
レーションを行うことが可能である。シミュレーションの例として、吸
収層にタングステンを使用して 10GeV の電子を 10 回入射させたイメー
ジが下の図 4 である。
図 4: Geant4 によって得られるイメージ
左が正面から見たとき、右が横から見たときのイメージである。狭い
範囲で行われるイベントでも、上図のように拡大してみることもできる。
そして、このようなイメージと同時に粒子が入射した位置とその位置で
検出したエネルギーを記録することができる。そのデータを用いて性能
評価を行っていく。
13
3.2
シミュレーションモデル
本研究で用いたモデルは以下のように作製する。
図 5: 電磁カロリメータ
吸収層には鉄（Fe)、鉛 (Pb)、タングステン (W) の 3 つの場合をシミュ
レーションする。それぞれの物質の原子番号 Z は、Fe=26、Pb=82、W=74
である。そして、検出層にはどの場合でもポリスチレン製のシンチレー
タを使用し、縦と横を 100 分割して 1cm×1cm のタイル状にする。
また、入射粒子の入射方向を入射角 θ と方位角 ϕ で表す。入射角と方位
角は下の図 6 のように定義することにする。
図 6: 入射角 θ と方位角 ϕ
14
3.3
角度毎のシミュレーション
入射角 θ は 0◦ ∼75◦ まで 15◦ ずつ変更し、方位角 ϕ は 45◦ で固定する。
入射エネルギーは 1∼5、10、20、30、40、50GeV と変更していき、入射
場所は (x, y) = (0, 0)、イベント数は 3000 回でシミュレーションする。吸
収層と検出層の厚さは、それぞれ 10mm、5mm に設定する。
3.3.1
位置分解能
本研究では、位置分解能としてエネルギー重み半径 Rweight と RMS を
を採用する。エネルギー重み半径 Rweight は次の式（13）で定義される。
∑
Rweight =
iRi Ei
iEi
∑
(13)
ここで Ri は入射軸からのタイルまでの距離、Ei はそのタイルで検出
したエネルギーである。エネルギー重み半径を用いて評価を行う理由は、
実際の距離で見たときに入射軸から遠く離れているタイルを除外するの
ではなく、そのタイルでの検出エネルギーの大きさによってデータの良
し悪しを決めるためである。測定は以下のように行う。
＜位置分解能の測定方法＞
1. 各入射エネルギーにおいて、1 イベント毎に各タイルでのエネルギー
重み半径を求めていき、縦軸にイベント数、横軸にエネルギー重み
半径をとったヒストグラムを作成する。
2. そのヒストグラムから Mean 値と RMS が得られるため、各入射エネ
ルギーにおいて縦軸にエネルギー重み半径、横軸に入射エネルギー
をとったグラフにプロットする。
それぞれの吸収層での結果は次頁のようになった。左がエネルギー重み
半径、右が RMS についてのグラフである。
15
図 7: Fe
図 8: Pb
図 9: W
16
図 7、8、9 の左側のグラフより、エネルギー重み半径は入射エネルギー
に依存しないことがわかる。また、図 7、8、9 の右側のグラフを入射エネ
ルギーの小さい方から見ていくと、はじめは傾きが大きいがだんだん小
さくなっていくことがわかる。そして、どの吸収層の場合でも入射角が
0◦ の場合がエネルギー重み半径、RMS ともに値が小さくなっていること
がグラフから読み取れる。
3.3.2
エネルギー分解能
一定の入射エネルギーを持った粒子をカロリメータに入射させると、1
イベント毎に異なったエネルギーがカロリメータ内部で検出される。そ
のエネルギーの誤差 σE を検出したエネルギーの平均値 EM ean で割ると、
エネルギー分解能 σE /EM ean となる。
また、一般的に検出される光量を N とすると、光量は電子や γ 線の持っ
ているエネルギー E と比例する (N ∝ E)。このことからエネルギー E を
√
測定することができ、その測定誤差は N である。本研究においては、
検出層で検出されるエネルギーが入射エネルギー E0 と比例していると言
√
える。したがって、測定誤差は E0 となるため、エネルギー分解能は次
の式（14）と表すことができる。
σ
EM ean
= √
a
⊕
E0 [GeV ]
b=
v
)2
u(
u
t √a
E0
+ b2
(14)
ここで、電磁カロリメータでは入射粒子のエネルギーを可能な限り正確
に精度良く測定したいため、エネルギー分解能は小さい方が良い。その
ため、a と b の値は小さくなればよい。測定方法は以下のとおりである。
＜エネルギー分解能の測定方法＞
1. 各入射エネルギーにおいて、検出層で検出したエネルギーを 1 イベ
ント毎にとり、ヒストグラムにする。そのヒストグラムをガウス関
数でフィットして、EM ean と σE を求める。
√
2. 各入射エネルギー毎の σE /EM ean と 1/ E0 を縦軸にエネルギー分
√
解能、横軸に 1/ E0 をとってプロットし、式 (14) でフィットして a
と b の値を求める。
結果は次頁からのようになった。（theta0 は入射角 0◦ のことである。）
17
図 10: 吸収層が鉄の場合の角度毎のエネルギー分解能
18
図 11: 吸収層が鉛の場合の角度毎のエネルギー分解能
19
図 12: 吸収層がタングステンの場合の角度毎のエネルギー分解能
20
図 10∼12 のグラフから、シミュレーションで得られたデータと式（14）
がよく合っていることがわかる。したがって、エネルギー分解能を求め
られていると言える。
ここで次に進むにあたり、変更した点がある。
それはイベントの数で、3.3 までは 3000 回としてシミュレーションして
きた。しかし、3.4 からは 1000 回としている。これは私の研究の失敗点
であるが、はじめにイベント数 1000 回のシミュレーションを行った際に、
良い値（実際値を入れる）が得られなかったため、イベント数を 3000 回
にして再度行ったところ、求めたい値（実際の値）が得られたため、角度
毎のシミュレーションにおいては 3000 回で行った。その後、長谷川先生
よりご指摘を頂いて、イベント数が 1000 回と 3000 回の場合の違いがあ
るかどうかエネルギー重み半径の値を用いて検証を行った。吸収層にタ
ングステンにして、１ GeV の電子を入射角 0◦ で入射させ、5 回分のデー
タを取った。結果は表 1 のようになった。
表 1: イベント数 1000 回と 3000 回でのエネルギー重み半径の比較
イベント数
1000 回
3000 回
1 回目 2 回目 3 回目 4 回目
11.79
11.9
11.75 11.82
12.05 11.96 11.76 11.87
(単位はすべて mm である)
5 回目
11.96
11.86
イベント数 1000 回では 11.85±0.08、3000 回では 11.90±0.10 となった。
この検証結果から、イベント数による違いはそれほど見られないことが
わかったため、以後のシミュレーションではイベント数 1000 回で行った。
21
3.4
入射位置を変更した場合
今までの解析では入射位置を電磁カロリメータの中央 ((x, y) = (0, 0))
としていた。そこで、入射位置を (x, y) = (10.5, 10.5) にした場合を考え
る。吸収層にはタングステンを用い、10GeV の電子を 1000 回入射させる。
入射角は θ = 0◦ である。変更前と変更後のエネルギー重み半径、EM ean
で比較すると、以下のようになった。
図 13: 入射場所変更前 (左) と変更後 (右) のエネルギー重み半径
図 14: 入射場所変更前 (左) と変更後 (右) の EM ean
それぞれの値をまとめると、次頁の表 2 のようになった。
22
表 2: 入射場所変更前と変更後の各測定値の比較
エネルギー重み半径 [mm]
EM ean [MeV]
変更前
12.1 284.3 変更後
12.93 281.8 大きな違いがみられないため、入射場所変更による影響はないと言える。
23
3.5
方位角を変更した場合
3.3 では、方位角による影響はないとして入射角を変更した場合をシ
ミュレーションを行った。しかし、実際にその通りかどうかを検証して
いなかったため、確認をする。方位角を ϕ = 0◦ 、15◦ 、30◦ 、45◦ の 4 つの場
合での位置分解能とエネルギー分解能を用いて議論することにする。こ
のとき、入射角は θ = 30◦ に固定する。検証結果は以下の通りになった。
(phi0 は ϕ = 0◦ のことである。)
図 15: 方位角毎のエネルギー重み半径
24
図 16: 方位角毎の Energy deposit の Mean 値
それぞれの値をまとめると、下の表 3 のようになった。
表 3: 方位角毎の各測定値の比較
方位角
エネルギー重み半径 [mm]
EM ean [MeV]
0◦
42.67 263.5 15◦
42.28 266.5 30◦
42.71 264.3 上記の結果から、方位角による影響はないと言える。
25
45◦
42.29 265 3.6
実用型モデルでのシミュレーション
ここからは実際に使用が考えられている、吸収層 3mm、検出層 2mm
での電磁カロリメータについてシミュレーションを行っていき、位置分解
能、エネルギー分解能の２つと新たな比較の指標としてモリエール半径
を導入することとする。入射場所 (x, y) = (0, 0) に 10GeV の電子を θ = 0
で入射して比較を行う。
（方位角は関連してこないため、特に設定してい
ない。)
3.6.1
位置分解能
それぞれの吸収層での結果は、図 17 のようになった。
図 17: それぞれの吸収層でのエネルギー重み半径 (左) と RMS（右）
まず、エネルギー重み半径だが入射エネルギーに依らないので、平均
値を取ると表 4 のようになった。
表 4: 吸収層毎のエネルギー重み半径の平均値
吸収層の物質
エネルギー重み半径の平均値 [mm]
Fe
17.2
Pb
W
14.2 11.3
表 4 より、最もエネルギー重み半径が小さくなるのはタングステンを用
いた場合であるといえる。また、RMS についても入射エネルギーの小さ
い方から見ていくと、20GeV まではほぼタングステンの場合が最も小さ
くなっているが、それ以上ではあまり明瞭な違いがみられなくなってい
る。したがって、位置分解能はタングステンを用いた場合が最も小さい。
26
3.6.2
エネルギー分解能
それぞれの物質での結果は、図 18 のようになった。
図 18: それぞれの吸収層でのエネルギー分解能
それぞれの物質での a と b の値をまとめたものが以下の表 5 である。
表 5: 吸収層毎の a とｂの値
吸収層の物質
a [％]
b [％]
Fe
6.7 ± 0.2
1.1 ± 0.2
Pb
11.5 ± 0.2
0.6 ± 0.2
W
16.3 ± 0.3
0.6 ± 0.4
表 3 からわかるように、b の値は誤差まで含めるとそれほど差がみられ
ない。一方、a の値は鉄の場合が最も小さいことがわかる。したがって、
エネルギー分解能は鉄を用いた場合が最も小さくなる。
27
3.6.3
モリエール半径
それぞれの吸収層の物質での放射長 X0 と臨界エネルギー Ec 、そして
モリエール半径の予想値は表 6 となる。
表 6: それぞれの物質での放射長と臨界エネルギーと得られる予想値
吸収層の物質
放射長 X0 [mm]
臨界エネルギー Ec [MeV]
予想値 [mm]
Fe
18 22
16.9
Pb
5.6 7.4
16.0
W
3.5 8.0
9.3
本研究では、モリエール半径を求めるために際して自身でプログラム
を作成した。求めるまでのアルゴリズムは以下のとおりである。
＜モリエール半径の測定方法＞
1. シミュレーションによって得られたデータから、各タイルまでの距
離 Ri とそのタイルで検出したエネルギー Ei のデータを取得する。
2. 同様のシミュレーションのデータを用いて、1 イベント毎にカロ
リメータ内部に落ちた全エネルギーのデータを取得する。
3. 付録にあるモリエール半径を求めるプログラムを用いて、1 で得た
データを 1 イベント毎に Ri の小さい順にソートを行う。
4. ソートが終わったら、Ri が小さい順に落ちたエネルギーを足しあげ
ていき、2 で得られた全エネルギーの 90 ％に達したときに終了し、
その際の Ri の値を出力する。この Ri がそのイベントにおけるモリ
エール半径となる。
5. 1000 イベント分のモリエール半径を求めたら、その平均値と誤差を
算出する。
上記の方法を行った結果は次頁の図 19 のようになった。
28
図 19: それぞれの吸収層でのモリエール半径
今回のシミュレーションの結果と予想値を並べたものが表 7 である。
表 7: 吸収層毎のエネルギー重み半径とその予想値
吸収層の物質
シミュレーション値 [mm]
予想値 [mm]
Fe
37 28
Pb
29 27
W
20 16
本研究での予想値は、純粋な物質ではなく検出層と交互に並べられて
いるので、吸収層と検出層の比 3 : 2 を考慮して、 53 倍している。表 3 よ
り、予想値とは離れていることがわかるが、モリエール半径はタングス
テンを用いた場合が最も小さくなる。
29
4
まとめ
本研究では、Geant4 を用いて電磁カロリメータを作製し、吸収層の物
質として鉄、鉛、タングステンを用いてシミュレーションを行った。
● 角度毎のシミュレーション
位置分解能は入射エネルギーに依らず、入射エネルギーが大きくな
れば一定になることがわかった。そして、θ = 0◦ のときに最小値を
取ることがわかった。
● 入射位置をした場合と方位角を変更した場合
どちらの場合でも、得られるデータに大きな違いは見られなかった
ため、入射位置や方位角による影響はない。
● 実用型モデルでのシミュレーション
エネルギー重み半径はタングステンのときに最小値をとった。RMS
についても 20GeV までは、ほぼタングステンを用いたときに最小
値をとったが、それ以上では違いがみられなかった。
エネルギー分解能は、鉄を用いた場合が最も良くなった。
モリエール半径は、タングステンを用いると最も小さくなった。し
かし、予想される値とは差が生じていた。
以上の結果から、実用型モデルを使用する際は、位置分解能を良く
する場合、あるいはモリエール半径を小さくしてイベント毎の検出
したエネルギーの測定を正確に行う場合には、タングステンを用い
た電磁カロリメータを作製し、エネルギー分解能を良くする場合に
は、鉄を用いた電磁カロリメータを作製すべきである。
30
5
今後の展望
● タイルの大きさの変更
本研究では、RMS の値が 20GeV 以上では違いがみられなかったり、
モリエール半径の測定値が予想値と離れてしまったりしていた。こ
の原因として考えられることは、タイルの大きさである。今回は検
出層を 1cm × 1cm のタイルに分割してシミュレーションを行った。
使用したプログラム上、図 20 の左のように実際の入射した場所が
違っても、同じタイル内であれば、同じタイル番号を返すように作
成されていた。これを解決するには、図 20 の右のように、タイル
の大きさを小さくしてシミュレーションをすることが最も簡単な改
善方法である。また、測定値を算出するプログラムを再度作り直す
方法もある。
図 20: 電磁カロリメータ
● 電磁カロリメータの縦（深さ）方向のシミュレーション
本研究では、電磁カロリメータの横方向までしか議論ができなかっ
た。縦 (深さ) 方向まで議論することができれば、鉄とタングステン
のどちらを使用すべきかがわかると考えられる。縦方向に関連して
くるパラメータは、放射長 X0 である。表 6 からわかるように、鉄
の放射長は 18mm である。一方、タングステンの放射長は 3.5mm
と短い。この関係がそのまま電磁カロリメータの厚さになると考え
ると、鉄で作製した場合の厚さは、タングステンで作製した場合の
約 5.1 倍となる。実際にシミュレーションを行えれば、最適な物質
がタングステンであることの根拠の一つとなると考えられる。
31
謝辞
本研究に取り組むにあたって、竹下徹先生、長谷川庸司先生、研究員の
小寺克茂さんには丁寧にご指導をして頂きました。毎週の報告会では適
切な助言を頂けたため、これから何をしていけばよいか、どのように改善
すればよいかをしっかりと確認することができました。特に、長谷川先生
には何度も質問に伺ってしまいましたが、そのたびに丁寧に対応をして
頂きました。また、大学院生の寺田怜真さん、土本航也さんには Geant4
やプログラミングのことについて多くのことを御教授頂きました。そし
て、モリエール半径の算出プログラム作成の際には、物性理論研究室の
田中優也さんに助力を頂きました。
お世話になった皆様に、この場を借りて御礼申し上げます。
最後に、大学に通わせてくれた両親に心より感謝致します。
本当に有難うございました。
32
参考文献
[1] 山田作衛、相原博昭、岡田安弘、坂井典佑、西川公一郎編
『素粒子物理学ハンドブック』（朝倉書店、2010）
[2] Ashok Das、Thomas Ferbel 著末包文彦、白井淳平、湯田春雄訳
素粒子・原子核物理学の基礎 ∼実験から統一理論まで∼ （共立出版、2011）
[3] 長島順清著『素粒子物理学の基礎 I』（朝倉書店、1998）
[4] 中尾知博『Geant4 によるカロリメータの性能評価』（2013）
[5] 牧野正樹『Jupiter によるカロリメータのシミュレーション』
（2007）
図 1：http://atlas.kek.jp/sub/photos/ATLAS/PhotoATLAS.html
33
補遺
位置分解能とエネルギー分解能のグラフについて、プロットをするデー
タを求めるために使用したグラフがあった。ここでは、実用型モデルに
関するグラフを載せることにするが、角度毎のシミュレーションでも同
様のグラフが得られている。実用型モデルの条件は、吸収層 3mm、検出
層 2mm、入射粒子は電子で ECAL の中央 (x, y) = (0, 0) に入射角 θ = 0◦
で入射させる。イベント数は 1000 回である。
位置分解能
位置分解能の測定方法 1 で、入射エネルギー毎の Mean 値と RMS を求
めるヒストグラムを作成していた。それぞれの物質でのグラフは図 21 か
ら図 23 のとおりである。
図 21: 吸収層が鉄の場合の入射エネルギー毎のエネルギー重み半径
34
35
図 22: 吸収層が鉛の場合の入射エネルギー毎のエネルギー重み半径
36
図 23: 吸収層がタングステンの場合の入射エネルギー毎のエネルギー重
み半径
37
38
エネルギー分解能
エネルギー分解能のグラフ作成の際に、まず入射エネルギーと検出さ
れるエネルギーに線型性があるかどうかを確認するためのグラフを用い
ていた。それぞれの物質でのグラフは以下のとおりである。
プロットしている値と誤差については、エネルギー分解能の測定方法 1 で
のガウスフィットによって求められた値を使用している。
図 24: それぞれの吸収層での E − E0 グラフ
上の図から線型性があることがわかったので、エネルギー分解能と入
射エネルギーの関係をグラフにする。その結果は次頁のグラフとなった。
39
図 25: それぞれの物質でのエネルギー分解能と入射エネルギーの関係
この図を見ても、どのような関係があるかが明確ではなかったため、本
√
研究では横軸を 1/ E0 にしたグラフを結果として載せた。
40
プログラム
本研究で使用したプログラムを載せる。
G4DetecterConstruction.cc
このプログラムを用いて、カロリメータの形状やそれぞれの層に用い
る物質、検出層の分割数などを設定することができる。
#include "B4DetectorConstruction.hh"
#include "G4Material.hh"
#include "G4NistManager.hh"
#include
#include
#include
#include
#include
#include
"G4Box.hh"
"G4LogicalVolume.hh"
"G4PVPlacement.hh"
"G4PVReplica.hh"
"G4GlobalMagFieldMessenger.hh"
"G4AutoDelete.hh"
#include
#include
#include
#include
"G4GeometryManager.hh"
"G4PhysicalVolumeStore.hh"
"G4LogicalVolumeStore.hh"
"G4SolidStore.hh"
#include "G4VisAttributes.hh"
#include "G4Colour.hh"
#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
G4ThreadLocal
41
G4GlobalMagFieldMessenger* B4DetectorConstruction::fMagFieldMessenger = 0;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
B4DetectorConstruction::B4DetectorConstruction()
: G4VUserDetectorConstruction(),
fAbsorberPV(0),
fGapPV(0),
fCheckOverlaps(true)
{
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
B4DetectorConstruction::~B4DetectorConstruction()
{
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
G4VPhysicalVolume* B4DetectorConstruction::Construct()
{
// Define materials
DefineMaterials();
// Define volumes
return DefineVolumes();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
void B4DetectorConstruction::DefineMaterials()
{
// Lead material defined using NIST Manager
G4NistManager* nistManager = G4NistManager::Instance(); //使
42
用する物質を定義
nistManager->FindOrBuildMaterial("G4_Pb");
nistManager->FindOrBuildMaterial("G4_Fe");
nistManager->FindOrBuildMaterial("G4_W");
nistManager->FindOrBuildMaterial("G4_GLASS_LEAD");
nistManager->FindOrBuildMaterial("G4_POLYSTYRENE");
// Liquid argon material
G4double a; // mass of a mole;
G4double z; // z=mean number of protons;
G4double density;
new G4Material("liquidArgon", z=18., a= 39.95*g/mole,
density= 1.390*g/cm3);
// The argon by NIST Manager is a gas with a different density
// Vacuum
new G4Material("Galactic", z=1., a=1.01*g/mole,density=
universe_mean_density,kStateGas, 2.73*kelvin, 3.e-18*pascal);
// Print materials
G4cout << *(G4Material::GetMaterialTable()) << G4endl;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
G4VPhysicalVolume* B4DetectorConstruction::DefineVolumes()
{
// Geometry parameters //カロリメータの形状を設定する
G4int nofLayers = 100; //層の数
G4double absoThickness = 10.*mm; //吸収層の厚さ
G4double gapThickness = 5.*mm; //検出層の厚さ
G4double calorSizeXY = 100.*cm; //層の大きさ
G4double layerThickness = absoThickness + gapThickness;
G4double calorThickness = nofLayers * layerThickness;
43
G4double worldSizeXY = 10. * calorSizeXY;
G4double worldSizeZ = 1.2 * calorThickness;
// Get materials
G4Material* defaultMaterial = G4Material::GetMaterial("Galactic");
G4Material* absorberMaterial = G4Material::GetMaterial("G4_W");
//吸収層の物質を設定
G4Material* gapMaterial = G4Material::GetMaterial("G4_POLYSTYRENE");
//検出層の物質を設定
if ( ! defaultMaterial || ! absorberMaterial || ! gapMaterial ) {
G4ExceptionDescription msg;
msg << "Cannot retrieve materials already defined.";
G4Exception("B4DetectorConstruction::DefineVolumes()",
"MyCode0001", FatalException, msg);
}
//
// World
//
G4VSolid* worldS
= new G4Box("World",
// its name
worldSizeXY/2, worldSizeXY/2, worldSizeZ/2); // its size
G4LogicalVolume* worldLV
= new G4LogicalVolume(
worldS,
// its solid
defaultMaterial, // its material
"World");
// its name
G4VPhysicalVolume* worldPV
= new G4PVPlacement(
0,
// no rotation
G4ThreeVector(0.,0.,0.), // at (0,0,0)
worldLV,
// its logical volume
44
"World",
0,
false,
0,
fCheckOverlaps);
//
//
//
//
//
its name
its mother volume
no boolean operation
copy number
checking overlaps
//
// Calorimeter //カロリメータについて
//
G4VSolid* calorimeterS
= new G4Box("Calorimeter",
// its name
calorSizeXY/2, calorSizeXY/2, calorThickness/2); // its size
G4LogicalVolume* calorLV
= new G4LogicalVolume(
calorimeterS,
// its solid
defaultMaterial, // its material
"Calorimeter");
// its name
new G4PVPlacement(
0,
G4ThreeVector(),
calorLV,
"Calorimeter",
worldLV,
false,
0,
fCheckOverlaps);
//
//
//
//
//
//
//
//
no rotation
at (0,0,0)
its logical volume
its name
its mother volume
no boolean operation
copy number
checking overlaps
//
// Layer //層について
//
G4VSolid* layerS
= new G4Box("Layer",
// its name
calorSizeXY/2, calorSizeXY/2, layerThickness/2); // its size
45
G4LogicalVolume* layerLV
= new G4LogicalVolume(
layerS,
// its solid
defaultMaterial, // its material
"Layer");
// its name
new G4PVReplica(
"Layer",
layerLV,
calorLV,
kZAxis,
nofLayers,
layerThickness);
//
//
//
//
//
//
its name
its logical volume
its mother
axis of replication
number of replica
witdth of replica
//
// Absorber //吸収層について
//
G4VSolid* absorberS
= new G4Box("Abso",
// its name
calorSizeXY/2, calorSizeXY/2, absoThickness/2); // its size
G4LogicalVolume* absorberLV
= new G4LogicalVolume(
absorberS,
// its solid
absorberMaterial, // its material
"Abso");
// its name
fAbsorberPV
= new G4PVPlacement(
0,
G4ThreeVector(0.,
absorberLV,
"Abso",
layerLV,
// no rotation
0., -gapThickness/2), // its position
// its logical volume
// its name
// its mother volume
46
false,
0,
fCheckOverlaps);
// no boolean operation
// copy number
// checking overlaps
//
// Gap //検出層について
//
G4int ndivx=100; //横の分割数
G4int ndivy=100; //縦の分割数
G4VSolid* gapS
= new G4Box("Gap",
// its name
calorSizeXY/ndivx/2, calorSizeXY/ndivy/2, gapThickness/2); // its size
G4LogicalVolume* gapLV
= new G4LogicalVolume(
gapS,
// its solid
gapMaterial,
// its material
"Gap");
// its name
for(G4int ix = 0;ix < ndivx; ++ix)
{
for(G4int iy = 0; iy < ndivy; ++iy)
{
fGapPV
= new G4PVPlacement(
0,
// no rotation
G4ThreeVector(calorSizeXY/ndivx*(ix-(ndivx-1)/2.),calorSizeXY/
ndivy*(iy-(ndivy-1)/2.) , absoThickness/2), // its position
gapLV,
// its logical volume
"Gap",
// its name
layerLV,
// its mother volume
false,
// no boolean operation
ix*100+iy,
// copy number
fCheckOverlaps); // checking overlaps
}
47
}
//
// print parameters
//
G4cout << "\n------------------------------------------------------------"
<< "\n---> The calorimeter is " << nofLayers << " layers of: [ "
<< absoThickness/mm << "mm of " << absorberMaterial->GetName()
<< " + "
<< gapThickness/mm << "mm of " << gapMaterial->GetName() << " ] "
<< "\n------------------------------------------------------------\n";
//
// Visualization attributes
//
worldLV->SetVisAttributes (G4VisAttributes::Invisible);
G4VisAttributes* simpleBoxVisAtt= new G4VisAttributes(G4Colour(1.0,1.0,1.0));
simpleBoxVisAtt->SetVisibility(true);
calorLV->SetVisAttributes(simpleBoxVisAtt);
//
// Always return the physical World
//
return worldPV;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
void B4DetectorConstruction::ConstructSDandField()
{
// Create global magnetic field messenger.
// Uniform magnetic field is then created automatically if
// the field value is not zero.
G4ThreeVector fieldValue = G4ThreeVector();
fMagFieldMessenger = new G4GlobalMagFieldMessenger(fieldValue);
48
fMagFieldMessenger->SetVerboseLevel(1);
// Register the field messenger for deleting
G4AutoDelete::Register(fMagFieldMessenger);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
必要なデータへの変換プログラム
中尾知博さんの卒業論文に掲載されていたプログラムを参考に、本研
究に合うように書き換えた。以下に載せるプログラムで、エネルギー重
み半径や 1 イベント毎の Enegy deposit を求めることができる。
{
gROOT->Reset();
const Double_t XYSIZE = 1000.; //カロリメータのサイズ
const int N_DIV =100; //分割数
const int N_EVENT =1000; //イベント数
const Char_t CONDITION[50] = "e-xW_0_10GeV_placechange";
Double_t Cell_SIZE = XYSIZE / (float)N_DIV;
Char_t ifname[100];
sprintf(ifname, "/home/isami/myg4work/B4-build/B4a
/Output_%s.dat",CONDITION);
ifstream ifs(ifname);
if(!ifs){
cout << "ERROR!:file is not found" << endl;
exit();
}
string line;
Int_t eventNo=0;
Int_t absID, gapID;
Double_t edep;
Int_t tileID;
49
Int_t xID, yID;
Double_t xcoord, ycoord;
typedef struct{
vector<Double_t> R;
vector<Double_t> EDEP;
Double_t total_edep;
Double_t Rweight;
}Radius;
Radius m_radius[N_EVENT];
m_radius[0].total_edep = 0;
while(!ifs.eof()){
getline(ifs,line);
if(line.find("-----") == string::npos){
istringstream iss(line);
iss >> absID
>> gapID
>> edep
>> tileID;
xID = (Int_t)(tileID / N_DIV); //x 方向のタイル番号
yID = tileID % N_DIV; //y 方向のタイル番号
if(absID == 1){
xcoord = Cell_SIZE * (xID-2.5) + Cell_SIZE * 0.5 -XYSIZE / 2.;
ycoord = Cell_SIZE * (yID-2.5) + Cell_SIZE * 0.5 - XYSIZE / 2.;
m_radius[eventNo].R.push_back(sqrt(pow(xcoord, 2.)+
pow(ycoord,2.))); //入射軸からの距離
m_radius[eventNo].EDEP.push_back(edep);
//検出したエネルギー
m_radius[eventNo].total_edep += edep;
}
}
else{
++eventNo;
m_radius[eventNo].total_edep = 0;
absID = 0;
}
50
}
ifs.close();
----------------------------------------------------------------Double_t RW=0, tmpE=0;
for(int i=0; i<N_EVENT; i++){
RW = 0;
for(int k=0; k<(int)m_radius[i].R.size(); k++){
RW += m_radius[i].R[k] * m_radius[i].EDEP[k];
}
m_radius[i].Rweight = RW / m_radius[i].total_edep;
}
Char_t ofname2[100];
sprintf(ofname2, "/home/isami/myg4work/Rweight_data-graph
/Rweight_%s.dat", CONDITION);
ofstream ofs2(ofname2);
ofs2 << "Rweight" << endl;
for(int i=0; i<N_EVENT; i++){
ofs2 << m_radius[i].Rweight << endl;
}
ofs2.close();
----------------------------------------------------------------}
上述したものは、エネルギー重み半径を求めるものである。点線の内側
を下のように変更すると、エネルギー分解能を求めるために必要な 1 イ
ベント毎の検出した全エネルギー（書き換え 1）、モリエール半径を求め
るために必要な入射軸からの距離 Ri とその場所での検出エネルギー Ei
（書き換え 2）を求めることができる。
（書き換え 1）
Char_t ofname2[100];
sprintf(ofname2, "/home/isami/myg4work/Rweight_data-graph
/Rweight_%s.dat", CONDITION);
ofstream ofs2(ofname2);
51
ofs2 << "TOTAL" << endl;
for(int i=0; i<N_EVENT; i++){
ofs2 << m_radius[i].total_edep << endl;
}
ofs2.close();
（書き換え 2）
Char_t ofname2[100];
sprintf(ofname2, "/home/isami/myg4work/Rweight_data-graph
/Rweight_%s.dat", CONDITION);
ofstream ofs2(ofname2);
ofs2 << "R-EDEP" << endl;
for(int i=0; i<N_EVENT; i++){
for(int j=0; j< m_radius[i].R.size(); ++j)
ofs2 << ‘‘-----’’ << endl;
}
ofs2.close();
モリエール半径
#include
#include
#include
#include
#include
<iostream>
<fstream>
<sstream>
<string>
<vector>
#define l 1000 //イベント数
using namespace std;
int main(){
const char name[50] = "5mmx5mm";
char RE[100];
sprintf(RE,"R-EDEP_%s.dat",name);
ifstream ifs(RE);
52
if(!ifs){
cout << "ERROR!:The file is not found." << endl;
return -1;
}
/* 1 イベントで落ちたエネルギーのデータを配列に格納し、すべての和
を求める */
string buf;
double R[l];
for(int i=0; i<=l; i++){
R[i] = 0.;
}
for(int n = 0; n <= l; n++){
double total=0.;
//ifs から
vector<double> r,e;
double bufe, bufr;
while(!ifs.eof()){
getline(ifs,buf);
if(buf == "-----")
break;
istringstream iss(buf);
iss >> bufr
>> bufe;
r.push_back(bufr);
e.push_back(bufe);
total += bufe;
}
cout << "total " << total << endl;
cout << "ソート開始" << endl;
double swap;
//ソート開始半径を昇順で並び替え
for(int j = 0; j < r.size(); j++){
for(int k = j + 1; k < r.size(); k++){
53
if(r[j] > r[k]){
swap = r[j];
r[j] = r[k];
r[k] = swap;
swap = e[j];
e[j] = e[k];
e[k] = swap;
}
}
}//ソート終了
cout << "ソート終了" << endl;
double sum = 0.;
double rate, a;
for(int i = 0; i < r.size(); i++){
sum += e[i];
rate = sum/total;
if(rate >= 0.9){
cout << "半径 = " << r[i] << endl;
a = r[i];
break;
}
}
r.clear();
e.clear();
R[n] = a;
cout << " 終了回数 : " << n << endl;
}
char ofname[100];
sprintf (ofname, "m-%s.dat",name);
ofstream ofs(ofname);
for(int i=0; i<=l; i++){
ofs << R[i] << endl;
}
ofs.close();
}
54

Download Report