宇宙線ミューオンしたい

 1PrimaryGeneratorAction::GeneratePrimaries(G4Event *aEvent)
 2{
 3    // イベントごとにランダムにしたいため
 4    // 1回のイベントで生成する粒子の数は1個にする
 5    G4int n_particle = 1;
 6    auto gun = new G4ParticleGun(n_particle);
 7
 8    // 入射粒子を指定
 9    G4ParticleTable *table = new G4ParticleTable::GetParticleTable();
10    G4ParticleDefinition *particle = table->FindParticle("mu-");
11    gun->SetParticleDefinition(particle);
12
13    // 入射位置を指定
14    G4ThreeVector xyz{0., 0., 0.};
15    gun->SetParticlePosition(xyz)
16
17    // 入射エネルギーを指定
18    G4double energy{1. * GeV};
19    gun->SetParticleEnergy(energy);
20
21    // 入射方向を指定
22    G4ThreeVector direction{0., 0., 1.};
23    gun->SetParticleMomentumDirection(direction);
24
25    // イベントにvertexを追加
26    gun->GeneratePrimaryVertex(aEvent);
27}

地表に降り注ぐ宇宙線は、そのほとんどがミュー粒子と考えられます。 Geant4のミューオンは、正ミューオン（mu+）、負ミューオン（mu-）がそれぞれ用意されています。

天頂角分布したい

 1#include "G4UniformRand.hh"
 2#include "G4ThreeVector.hh"
 3#include "G4PhysicalConstants.hh"
 4
 5G4ThreeVector RandomMuonDirection() const
 6{
 7    // 天頂角（zenith angle）
 8    G4double u = G4UniformRand();
 9    G4double theta = std::acos(std::sqrt(u));
10
11    // 方位角（azimuth angle）
12    G4double v = G4UniformRand();
13    G4double phi = CLHEP::twopi * v;
14
15    // 運動方向を計算
16    G4double sinTheta = std::sin(theta);
17    G4double cosTheta = std::cos(theta);
18    G4double sinPhi = std::sin(phi);
19    G4double cosPhi = std::cos(phi);
20
21    G4double x = sinTheta * cosPhi;
22    G4double y = sinTheta * sinPhi;
23    G4double z = cosTheta;
24
25    G4ThreeVector direction{x, y, z};
26    return direction;
27};

宇宙線が降り注ぐ方向は \(\cos^{2}\theta\) に比例することが知られています。

\[\begin{split}R & := \cos^{2}\theta \\ \cos \theta & = \sqrt{R} \\ \theta & = \cos^{-1} R = \arccos R \end{split}\]

天頂角 \(\theta\) と 方位角 \(\phi\) の値を使って、 方向ベクトルは次のように計算できます。

\[\begin{split}x & = \sin\theta \cos\phi \\ y & = \sin\theta \sin\phi \\ z & = \cos\theta \end{split}\]