Geant4 による カロリメータのシミュレーション 信州大学理学部物理科学科高エネルギー研究室 02S2030E 堀越健作

Geant4 による

カロリメータのシミュレーション

信州大学理学部

物理科学科

高エネルギー研究室

02S2030E

堀越

健作

目次

第一章 ～はじめに～ --- １ 1, 目的 第二章 カロリメータ内の現象 --- ２ 1, 電磁シャワー 第三章 シミュレーション --- ５ 1, 概要 ～Geant4 とは～ 2, Geant4 の特徴 3, exampleN03 のカロリメータ 4, カロリメータに空気を入れる（Gap の数を増やす） 5, 測定器 6, 入射粒子 7, 磁場 8, シミュレーションの流れ 第四章 シミュレーションの結果 --- １２ 1, Pb8mm.sci2mm 2, Pb8mm.sci2mm.Air1mm 3, Pb8mm.sci2mm.Air2mm 4, Pb8mm.Air1mm.sci2mm.Air1mm 5, Pb8mm.Air2mm.sci2mm 第五章 結果に基づいたデータ処理・考察 --- １９ 1, カロリメータ 1、2 について 2, カロリメータ 3、4、5 について 第六章 まとめ --- ２６ 第七章 付録 --- ２７

１

第一章 ～はじめに～

目的

カロリメータに粒子を入れ、その振舞いを観測する。またカロリメータに磁場をかけ、今 後実際の実験に備えてシミュレートする。理由として現実にメートル単位で磁場をかける のは容易ではないことが挙げられる。そのためシミュレーションをすることが必要である。 そしてこの研究結果を今後の測定器開発に役立てる。 Geant4 の一般的なカロリメータについてシミュレーションする。Geant4 については三章で 触れる。

２

第二章 カロリメータ内の現象

電磁シャワー

物質に高エネルギーの電子や陽電子が入射し原子核のそばを通ったとき、原子核の電場に よりγ線を放出する(制動放射)。放出されたγ線は、電子と陽電子を作る(電子・陽電子対 生成) 。 このように制動放射と電子対生成の繰り返しで、電子と陽電子がネズミ算式に増えていく 現象を電磁シャワーという。シャワーというのは粒子の飛散の様子から呼ばれている。 このシャワー現象は電子らのエネルギーが以下で表されるＥc（臨界エネルギー）になるま で続く。 （Z は物質の原子量） ●制動放射 高エネルギー荷電粒子が物質中の電場の影響を受け減速して 1 個の光子を生成する。 100MeV 以上の電子、陽電子ではこの反応が主である。 ●対生成 原子核による強い電場の影響で光子が消滅して電子、陽電子を生成する。電子と陽電子の 質量は 1.022MeV に相当するため対生成には 1.022MeV 以上のエネルギーが必要となる。

対生成

Z + 1.2

800

MeV

Ｅc ＝

制動放射

電子

光子（γ線）

陽電子

３ また、陽電子の速度が十分遅くなると周りの電子に引き寄せられ反応し、消滅して 2 個の 光子を生成するという対消滅反応が起こりやすくなる。 シャワー粒子数が最大となる層の厚さは ln(E0/Ec) その粒子数は 0.3E0/Ec*(E0/Ec)**-1/2 シャワーの横方向の拡がり x は次のように表される。 <x**2>**1/2～0.8Es/Ec E0 入射エネルギー、Ec 臨界エネルギー、Es=21MeV である。 後述するが今回のシミュレーションの場合、 シャワーを起こすターゲットは鉛（原子量 207.2）、入射する電子は E0 = 4GeV であるから Ec = 800/(207.2 + 1.2) = 3.83877159MeV およそ 3.8MeV 以上の電子がシャワーを起こす。 またシャワー粒子数が最大となる層の厚さは ln(E0/Ec)=ln(4000 / 3.83877159) = 6.94889722 mm 粒子数は 0.3E0/Ec*(E0/Ec)**-1/2 =((0.3 * 4 000) / 3.83877159) * sqr(3.83877159 / 4 000) = 9.68400744 シャワーの横方向の拡がり x は<x**2>**1/2～0.8Es/Ec=(0.8 * 21) / 3.83877159 = 4.3764 で表される。

全断面積

これらの過程の起こりやすさを考える。全断面積 σ とは、入射粒子が衝突する物質の標的、 ターゲットの面積に相当する。この値 σ と物質の密度により反応が起こる確率が求まり、 粒子の平均自由行程も求めることができる。 例えば粒子 B からなる板に粒子 A を衝突させるとする。 板の粒子数密度を n、厚さを d、よって単位面積あたり nd 個の粒子がある B の拡がりをσとして衝突粒子 A から板を見たとき、進路上に B がある割合は P = σnd この σ = P / nd を全断面積という。ここでは粒子 B は鉛であり、密度は 11.34*g/cm3 粒子 A は電子である。

4

平均自由行程

粒子が物質中で衝突せずに進む距離 x の平均値を平均自由行程 l といい、 l = 1 / nσ である。（n は粒子の数密度） 鉛の数密度 n は密度 11340kg/m3 原子量 207.2 から n = (11340 / 0.2072) * 6.02 * (10 ** 23) = 3.29472973 × 1028_m-3 となる。

５

第三章 シミュレーション

概要

～Geant4 とは～

Geant4 は C++言語を用いた高エネルギー大規模測定器シミュレータで、陽子・電子・中性 子・γ線・π粒子・μ粒子といった粒子と物質の複雑な相互作用をシミュレーションでき る。多くの分野で役に立ち、実際に宇宙開発・放射線医療などでも本格的な使用が開始さ れている。

Geant4 の特徴

・実験・観測結果との比較 ・複雑な物理シミュレーションを行える ・目的に応じてユーザーがソースを作成、修正できる

exampleN03 のカロリメータ

exampleN03 のカロリメータは、正方形の層(layer)がサンドイッチ状に並んだ構造をしてお り、鉛や鉄などを使う吸収層(Absorber)と、シンチレータを主とする検出層(Gap)がある。 カロリメータに入射した粒子は吸収層で前述の電磁シャワーを起こす。そして発生した荷 電粒子が検出層で出す光により、粒子の入射時刻やその種類・エネルギー値についての情 報を得ることが可能となる。 1layer （層）

鉛層

シンチレータ層

６

カロリメータに空気を入れる（Gap の数を増やす）

基本的にカロリメータは吸収層(Absorber)と検出層(Gap)の二つから成る。しかし現実にカ ロリメータを作ろうとする時にはどうしても空気が入ることがある。そこで本研究ではこ れを想定し、空気を含んだカロリメータも扱う。 そのために検出層(Gap)の数を増やすようソースを修正・作成する。具体的に述べると以下 のファイルである。 src/ExN03DetectorConstruction.cc src/ExN03DetectorMessenger.cc include/ExN03DetectorConstruction.hh include/ExN03DetectorMessenger.hh src/ExN03SteppingAction.cc

これらの中にある『Gap』と書かれた部分を、『Gap1』、『Gap2』として Gap1 をシンチレータ、 Gap2 を空気と定義する（逆の定義も可能）。

更に Gap3、Gap4・・・と加えれば必要に応じて複雑なカロリメータの作成が可能になる。

しかし src/ExN03SteppingAction.cc では以下に注意する。 if ((edep!=0.)) {

if (volume == detector->GetAbsorber()) eventaction->AddAbs(edep,stepl); if (volume == detector->GetGap1()) eventaction->AddGap(edep,stepl); if (volume == detector->GetGap2()) eventaction->AddGap(edep,stepl); ・・・以下略

この三行目と四行目

if (volume == detector->GetGap1()) eventaction->AddGap(edep,stepl); if (volume == detector->GetGap2()) eventaction->AddGap(edep,stepl);

はどちらかだけを使う。例えば Gap1（シンチレータ）の結果が欲しいのなら Gap2(空気)の 行を無効にしないと値が得られない。 もう一つ必要な作業として、位置(position)の設定がある。 例として src/ExN03DetectorConstruction.cc 内の Gap1 を挙げる。 ここでは 1layer は （鉛）Absorber/（シンチレータ）Gap1/（空気）Gap2 で構成されているとする。

７

// // Gap1 //

solidGap1=0; logicGap1=0; physiGap1=0;

if (Gap1Thickness > 0.)

{ solidGap1 = new G4Box("Gap1",Gap1Thickness/2,CalorSizeYZ/2,CalorSizeYZ/2); logicGap1 = new G4LogicalVolume(solidGap1,

Gap1Material,

Gap1Material->GetName());

physiGap1 = new G4PVPlacement(0, //no rotation G4ThreeVector((AbsorberThickness-Gap2Thickness)/2,0.,0.), //its position

logicGap1, //its logical volume

Gap1Material->GetName(), //its name logicLayer, //its mother false, //no boulean operat

0); //copy number }

この//its position と書かれた行に着目する。 1layer は Absorber/Gap1/Gap2 の順で並んでいる。

Gap1 の x 軸の position が(AbsorberThickness-Gap2Thickness)/2 と表されるのは、Gap1 よ りも前にあるものを正、後ろにあるものを負とすることで幾何学的に正しく層が配置され るためである。従って他の場合、 Absorber/Gap1/Gap2/Gap3/Gap4 であれば、Gap2 の x-position は (AbsorberThickness+Gap1Thickness-Gap3Thickness-Gap4Thickness)/2 となる。 これにより 1layer 内の構造が正しく設定され、必要とするカロリメータが用意できる。

８

測定器

カロリメータの断面は、一辺 100 cm の大きさをしている。 layer 数は、100layer である。これらは src/ExN03DetectorConstruction.cc で変更可能。

入射粒子

エネルギー4GeV の電子を入射する。 N03/run1.mac で決める。

磁場

粒子の入射方向に対して、垂直に 0 tesla から 10 tesla の磁場をかける。 そのために、exampleN04 のプログラムをコピーする。 １、N03/src 内に以下を ExN04Field.cc というファイル名で保存する。 #include "ExN04Field.hh" ExN04Field::ExN04Field() { Bz = 3.0*tesla; rmax_sq = sqr(50.*cm); zmax = 100.*cm; } ExN04Field::~ExN04Field() {;}

void ExN04Field::GetFieldValue(const double Point[3],double *Bfield) const {

Bfield[0] = 0.; Bfield[1] = 0.;

if(abs(Point[2])<zmax && (sqr(Point[0])+sqr(Point[1]))<rmax_sq) { Bfield[2] = Bz; }

100cm

100cm

磁場

100layer（層） 電子が入射する

９ else { Bfield[2] = 0.; } } ２、N03/include 内に以下を ExN04Field.hh というファイル名で保存する。 #ifndef ExN04Field_H #define ExN04Field_H 1 #include "globals.hh" #include "G4MagneticField.hh"

class ExN04Field : public G4MagneticField {

public:

ExN04Field(); ~ExN04Field();

void GetFieldValue( const double Point[3],

double *Bfield ) const;

private: G4double Bz; G4double rmax_sq; G4double zmax; }; #endif ３、以下を src/ExN03DetectorConstruction.cc に追加する。 #include "ExN04Field.hh" #include "G4TransportationManager.hh"

１０ #include "G4MaterialTable.hh" #include "G4Element.hh" #include "G4ElementTable.hh" #include "G4Tubs.hh" #include "G4ThreeVector.hh" #include "G4Transform3D.hh" #include "G4RotationMatrix.hh" #include "G4FieldManager.hh" 以下を G4VPhysicalVolume* ExN03DetectorConstruction::Construct()に追加する。 //---// Magnetic field

static G4bool fieldIsInitialized = false; if(!fieldIsInitialized)

{

ExN04Field* myField = new ExN04Field; G4FieldManager* fieldMgr = G4TransportationManager::GetTransportationManager() ->GetFieldManager(); fieldMgr->SetDetectorField(myField); fieldMgr->CreateChordFinder(myField); fieldIsInitialized = true; } これで ExN04Field.cc を使って磁場をかけることができる。

１１

シミュレーションの流れ

１、カロリメータに電子を一つ入れて、 シンチレータ 100layer に吸収されたエネルギーを合計する ２、電子 1000 個分についてエネルギー合計のヒストグラムを作り、 ３、それぞれ磁場 0 ～ 10 tesla でシミュレーションする ４、これをカロリメータ １～５ でそれぞれ行なう 理由 ・ランダム要素の大きい現象のため、1000event 分のデータを必要とする ・10tesla という値はかなり大きいが、どんな影響が出るか不明なので扱っておく。 ・前述した通り、空気の浸入を考慮したカロリメータ ２～５ を用意した。 色々な位置に、厚さの違う空気が入っている。2mm という値はかなり大きいが、磁場と同 じく扱ってみることで、影響を観察する。 カロリメータ １～５ は以下のように 1layer を構成する。 １、鉛(Pb)8mm ┃ シンチレータ(sci)2mm ２、Pb8mm ┃ sci2mm ┃ 空気(Air)1mm ３、Pb8mm ┃ sci2mm ┃ Air2mm

４、Pb8mm ┃ Air1mm ┃ sci2mm ┃ Air1mm ５、Pb8mm ┃ Air2mm ┃ sci2mm

１２

第四章 シミュレーションの結果

１、Pb8mm ┃ sci2mm のシンチレータ 100layer に吸収されたエネルギーの合計値

(電子 1000event)磁場の範囲は 0～10tesla 横軸【エネルギー合計値(MeV)】、縦軸【頻度】

■磁場 0tesla ■磁場 1tesla

■磁場 2tesla ■磁場 3tesla

１３ ■磁場 6tesla ■磁場 7tesla ■磁場 8tesla ■磁場 9tesla ■磁場 10tesla これらガウス曲線の頂点の値を 表・グラフにしたものを次に示す。

１４ １、Pb8mm ┃ sci2mm の磁場 0～10tesla の結果 磁場(tesla) 0 1 2 3 4 5 6 7 8 9 10 結果(MeV) 95.97 97.99 98.66 100.1 102.2 103.1 104.5 105.4 104.8 104.5 104.6 Sigma 0.34 0.33 0.34 0.3 0.3 0.3 0.3 0.4 0.3 0.4 0.4 ■横軸【磁場(tesla)】、縦軸【結果(MeV)】■

１５

２、Pb8mm ┃ sci2mm ┃ 空気(Air)1mm のシンチレータ 100layer に吸収されたエネルギー の合計値 (電子 1000event) 磁場の範囲は 0～10tesla ■磁場 0tesla（1～10tesla のグラフは省略）横軸【エネルギー合計値(MeV)】、縦軸【頻度】 磁場(tesla) 0 1 2 3 5 6 7 8 9 10 結果(MeV) 96.95 98.15 101.6 105.1 110.8 112 113.7 112.6 114.3 112.2 Sigma 0.33 0.34 0.3 0.3 0.4 0.4 0.4 0.4 0.4 0.4 ■横軸【磁場(tesla)】、縦軸【結果(MeV)】■

１６

３、Pb8mm ┃ sci2mm ┃ Air2mm のシンチレータ 100layer に吸収された エネルギーの合計値 (電子 1000event) 磁場の範囲は 0～10tesla ■磁場 0tesla（1～10tesla のグラフは省略）横軸【エネルギー合計値(MeV)】、縦軸【頻度】 磁場(tesla) 0 1 2 3 5 6 7 8 9 10 結果(MeV) 99.52 100.1 104.7 109.6 116.9 117.9 118.4 118.1 117.3 114.8 Sigma 0.33 0.3 0.4 0.4 0.4 0.4 0.4 0.5 0.5 0.5 ■横軸【磁場(tesla)】、縦軸【結果(MeV)】■

１７

４、Pb8mm ┃ Air1mm ┃ sci2mm ┃ Air1mm のシンチレータ 100layer に吸収された エネルギーの合計値 (電子 1000event) 磁場の範囲は 0～10tesla ■磁場 0tesla（1～10tesla のグラフは省略）横軸【エネルギー合計値(MeV)】、縦軸【頻度】 磁場(tesla) 0 1 2 3 5 6 7 8 9 10 結果(MeV) 96.92 97.57 100.5 103.3 106.9 106.9 106.3 104.5 102.6 100.1 Sigma 0.33 0.34 0.4 0.4 0.4 0.4 0.4 0.5 0.4 0.5 ■横軸【磁場(tesla)】、縦軸【結果(MeV)】■

１８

５、Pb8mm ┃ Air2mm ┃ sci2mm のシンチレータ 100layer に吸収された エネルギーの合計値 (電子 1000event) 磁場の範囲は 0～10tesla ■磁場 0tesla（1～10tesla のグラフは省略）横軸【エネルギー合計値(MeV)】、縦軸【頻度】 磁場(tesla) 0 1 2 3 4 5 6 7 8 9 10 結果(MeV) 96.96 95.05 95.07 95.08 94.94 93.02 91.73 89.81 88.29 85.71 83.2 Sigma 0.33 0.33 0.33 0.37 0.38 0.35 0.39 0.39 0.40 0.42 0.4 ■横軸【磁場(tesla)】、縦軸【結果(MeV)】■

１９

第五章 結果に基づいたデータ処理・考察

■ １、Pb8mm ┃ sci2mm

２、Pb8mm ┃ sci2mm ┃ Air1mm

について

横軸【磁場】、縦軸【結果】 磁場の増加に伴ってエネルギーの値が増えていき、磁場がある程度大きさになると、あま り変化は見られなくなった。 カロリメータに磁場をかけると荷電粒子の進行方向が曲がる。従って磁場がない時よりも カロリメータ内を走る距離が長くなるためエネルギーが多く吸収された（下図参照）。 6 tesla 以降エネルギーが増えないのは、電子の走る距離に大きな変化がないものと推測 される。これらの裏付けとして、次にカロリメータ１を使って粒子の移動距離のデータを とった。

Pb

sci

0 tesla

2 tesla

磁場

１ ２

２０ src/ExN03SteppingAction.cc 内を修正して を出力させた。 ■カロリメータ１(Pb8mm ┃ sci2mm）内を粒子の通過総距離 磁場の範囲は 0～10tesla ■磁場 0tesla（1～10tesla のグラフは省略）横軸【粒子の総距離(mm)】、縦軸【頻度】■ 磁場(tesla) 0 1 2 3 4 5 6 7 8 9 10 総距離(mm) 3523 3533 3536 3537 3553 3557 3561 3563 3561 3561 3557 Sigma 1.7 1.5 1.7 1.6 1.6 1.7 1.7 1.7 1.8 1.8 1.8 ■横軸【磁場(tesla)】、縦軸【総距離(mm)】■

２１ この結果から、電子がカロリメータを通過する距離が、磁場の増加に比例していることが 確かめられる。 6 tesla 以降距離に大きな変化がない理由は、磁場によりＵターンしている電子の数が増 えてきているためと考えられる。 電子のＵターンについては次の考察で触れる。

２２

■３、Pb8mm ┃ sci2mm ┃ Air2mm

■４、Pb8mm ┃ Air1mm ┃ sci2mm ┃ Air1mm

■５、Pb8mm ┃ Air2mm ┃ sci2mm

について

横軸【磁場】、縦軸【結果】 ４と５は磁場を強くすると値が大きく減少する結果を得た。また空気の位置が違うだけの ３と５の結果の差は注目したい。 ３と５の違いの理由の一つとして、磁場を上げることで鉛を通過した電子がシンチレータ に到達することなく空気層でＵターンし、再び鉛に入ってしまったことが考えられる（次 の図を参照）。 ３ ４ ５

２３ 電子がＵターンするモデル図 運動量、円運動半径、磁場の間には次の関係がある。これを使って暫定ではあるがこの現 象を考察する。 P = 0.3 R B P 運動量(MeV/c)、R 半径(mm)、B 磁場(tesla) 例 R = 2mm（Air の厚さ）、B = 1tesla とする。 このとき P = 0.3 x 2 x 1 = 0.6 MeV/c となり、0.6MeV/c の電子はシンチレータ層に到達することができない。 この仮説の裏付けとして、第一 layer 目の鉛層に吸収されたエネルギーを測定した。 仮説が正しければ磁場でＵターンした電子により、磁場に比例して値が増加すると予想さ れる。 これに必要なシミュレーションを次に示す。

磁場

2mm

Pb

Air sci

電子

２４

layer 毎の値をとるためには次のことをする。 第一 layer 目の鉛層を Abs1 と呼ぶことにする。

１、src/ExN03EventAction.cc では、

// initialisation per event と書かれたところに次の行を加える。 EnergyAbs1 = 0.;TrackLAbs1 = 0.;

２、src/ExN03EventAction.hh では Abs1 の定義を加える。

３、src/ExN03SteppingAction.cc では、以下のようにする。 if ((edep!=0.)) {

if (volume == detector->GetAbsorber()) eventaction->AddAbs(edep,stepl); // if (volume == detector->GetGap1())

eventaction->AddGap(edep,stepl);

if (volume == detector->GetGap2()) eventaction->AddGap(edep,stepl);

if(volume == detector->GetAbsorber() && -592 > track->GetPosition().x() && -600 < track->GetPosition().x()){

eventaction->AddAbs1(edep,stepl); }

２５ これによって、Abs1 の領域を定義してそこに吸収されたエネルギーが得られる。 一つの磁場に対して電子 1000event の結果を以下に示す。 ■Abs1 に吸収されたエネルギー値 磁場の範囲は 0～10tesla 磁場(tesla) 0 1 2 3 5 6 7 8 9 10 結果(MeV) 51.67 57.89 56.28 57.92 60.5 67.98 67.75 76.25 77.63 82.28 Sigma 1.27 1.18 1.26 1.09 1.3 1.18 1.50 1.44 1.27 1.40 ■横軸【磁場(tesla)】、縦軸【結果(MeV)】■ この結果から、磁場の増加に伴って第一 layer 目の鉛層に吸収されるエネルギーが増加す ることが確かめられた。先程の式 P = 0.3 R B P 運動量(MeV/c)、R 半径(mm)、B 磁場(tesla) から一定の運動量の時、磁場と回転半径は B = P / 0.3 R = K / R （K は定数）となり、 反比例している。磁場が増加すれば半径は小さくなり、Ｕターンしやすい。 従って磁場の影響からＵターンする電子が多くあると言える。 なお、このＵターン現象は空気層の有無に関わらず、物資中でも起こることが、前述の隙 間のないカロリメータ１、２の考察から言える。空気層を入れたことで鉛とシンチレータ の間隔が大きくなったため、Ｕターン現象が顕著に現れるようになったのである。

２６

第六章 まとめ

１、基本的に磁場をかけることで得るエネルギーが増加するが、空気の侵入と強すぎる磁 場はエネルギーの減少に繋がることが確認された。 ２、鉛とシンチレータの間には、できる限り空気を入れないようにするべきである。 ３、今回は扱った空気層の厚さは 1mm と 2mm だけだったが、今後 0.1mm、0.01mm オーダー でシミュレーションを行えば、磁場により値が減少する臨界の厚さが求まるだろう。 ４、これまでこのカロリメータに磁場をかけた実験は行なわれていない為、今回のシミュ レーションを今後の測定器開発に役立てたい。 最後に、本研究に助言を下さった竹下先生、同研究室の小島君、中山さんに感謝します。

２７

第七章 付録

値をファイルに落とす方法

src/ExN03EventAction.cc に以下を加える ofstream fo; fo.open("filename",std::ios::app); if(!fo) {

G4cout << G4endl << " can't !!!! " << G4endl << G4endl; exit(1); } if (evtNb%printModulo == 0) { fo << std::setw(7) << EnergyAbs << " " << std::setw(7) << EnergyGap << G4endl; fo.close(); }

これで EnergyAbs と EnergyGap の結果が filename に出力される。