Geant4のGeometryの最もシンプルな記述法

一辺2メートルの立方体の空気(不可視)に、一辺2メートルの立方体のアルゴンガスを詰めた例

Geant4で定義されている物質一覧

G4VPhysicalVolume* Geometry::Construct()
{
    G4NistManager* materi_Man = G4NistManager::Instance();

    G4double leng_X_World = 2.0 * m;
    G4double leng_Y_World = 2.0 * m;
    G4double leng_Z_World = 2.0 * m;

    G4VSolid* solid_World = new G4Box("Solid_World", leng_X_World / 2.0, leng_Y_World / 2.0, leng_Z_World / 2.0);
    G4Material* materi_World = materi_Man->FindOrBuildMaterial("G4_AIR");
    G4LogicalVolume* logVol_World = new G4LogicalVolume(solid_World, materi_World, "LogVol_World");
    logVol_World->SetVisAttributes(G4VisAttributes::GetInvisible());
    G4VPhysicalVolume* physVol_World = new G4PVPlacement(G4Transform3D(), "PhysVol_World", logVol_World, 0, false, 0);

    G4VSolid* solid_Box = new G4Box("Solid_Box", leng_X_World / 2.0, leng_X_World / 2.0, leng_X_World / 2.0);
    G4Material* materi_Box = materi_Man->FindOrBuildMaterial("G4_Ar");
    G4LogicalVolume* logVol_Box = new G4LogicalVolume(solid_Box, materi_Box, "LogVol_Box", 0, 0, 0);
    G4VPhysicalVolume* PhysVol_Box = new G4PVPlacement(G4Transform3D(), "PhysVol_Box", logVol_Box, physVol_World, false, 1000, true);
    //上記のポインタは使わないので下記のように表現されることが多い
    //new G4PVPlacement(G4Transform3D(), "PhysVol_Box", logVol_Box, physVol_World, false, 1000, true); 

    return physVol_World;
}

Geant4で重イオンをPrimaryGeneratorで入射する

Geant4で重イオンを入射したいという欲求に駆られた。Geant4で標準で扱える入射粒子はZが2以上だとヘリウム-4のみである。それ以外の原子核はやや特殊な方法で作る。

class PrimaryGenerator : public G4VUserPrimaryGeneratorAction
{
public:
    PrimaryGenerator();
    ~PrimaryGenerator();

public:
    void GeneratePrimaries(G4Event*);

private:
    G4ParticleGun * fpParticleGun;
};

G4VUserPrimaryGeneratorAction を継承したPrimaryGenerator クラスがあったとして、

PrimaryGenerator::PrimaryGenerator()
    : fpParticleGun(0)
{
    fpParticleGun = new G4ParticleGun();
    G4ParticleDefinition* particle = G4IonTable::GetIonTable()->GetIon(90, 232);
    fpParticleGun->SetParticleDefinition(particle);
    fpParticleGun->SetParticleEnergy(250 * MeV);
}

PrimaryGeneratorクラスのコンストラクタ内で重イオン(Z=90、A=232のトリウム)を作ろうとすると、

-------- WWWW ------- G4Exception-START -------- WWWW -------
*** G4Exception : PART105
      issued by : G4IonTable::CreateIon()
Can not create ions because GenericIon is not ready
*** This is just a warning message. ***
-------- WWWW -------- G4Exception-END --------- WWWW -------

-------- EEEE ------- G4Exception-START -------- EEEE -------
*** G4Exception : Event0101
      issued by : G4ParticleGun::SetParticleDefinition()
Null pointer is given.
*** Fatal Exception *** core dump ***
 **** Track information is not available at this moment
 **** Step information is not available at this moment

-------- EEEE -------- G4Exception-END --------- EEEE -------

以上のようなエラーが出る。後半はNull pointerだよと言ってるだけであまり意味はなくて、前半のGenericIon がないよというのが本命のエラー。何らかの準備が必要なのだろう(注: コメントで /run/initialize の前ではGenericIonができていないので使えませんと教えてもらいました)。解決策は見つからなかったので、

void PrimaryGenerator::GeneratePrimaries(G4Event* anEvent)
{
    G4ParticleDefinition* particle = G4IonTable::GetIonTable()->GetIon(90, 232);
    fpParticleGun->SetParticleDefinition(particle);
    fpParticleGun->SetParticleEnergy(250*A * MeV);
    //方向の設定など

    fpParticleGun->GeneratePrimaryVertex(anEvent);
}

GeneratePrimaries内で作ってやるとエラーは出なかった。

Physicsが正しいかどうかは別にして、Geant4は重イオンも扱える(ターゲットが重い原子核だと反応して出てくるし)ので、Primaryに出来ないはずはないが、手元のコードではPrimaryGeneratorでは作れなかった。

正しい重イオンビームの作り方を知っている人はこっそり教えて下さい。(注: このやり方で正しいらしいです)

PythonでJPGの2値画像をPNG画像に変換する

元々2値のPNG画像を無理やりJPG画像にすると、アナログノイズ的に輝度値が揺らいでしまう。 また、JPG画像の方が画像サイズが増えてしまうこともある。 JPG画像にしても輝度値は大きく変わらないので、閾値127で2値画像を作り、PNG形式で保存する関数を作った。

# %%
import cv2

def GetBinaryPng(jpgfile, pngfile, showhist=False):
    img = cv2.imread(jpgfile,0)
    hist = cv2.calcHist([img],[0],None,[256],[0,256])
    if showhist:
        import matplotlib.pyplot as plt
        plt.bar([i for i in range(256)], [a[0] for a in hist], width=1.0)
        plt.yscale('log')
        plt.show()

    ret,thresh1 = cv2.threshold(img,127,255,cv2.THRESH_BINARY)
    cv2.imwrite(pngfile, thresh1)

# %%
for i in range(0, 30):
    GetBinaryPng(f"{i}.jpg",f"{i}.png")
# %%

Python + pandas でcsvデータの一部を除去して保存する

pandas で読んだデータの一部を除去して保存する方法。下の例ではdata1という変数が1000以下のものを除去している。

import pandas as pd
df = pd.read_csv('run0001.csv')

index_drop = []
for i, row in enumerate(dfdump.itertuples()):
    if row.data1<1000:
        index_drop.append(i)
        continue

df_valid = df.drop(index_drop) #dfは変化しない
df_valid.to_csv('run0001_valid.csv')

Python (numpy+scipy)でプロファイルヒストグラムを描画する

プロファイルヒストグラム(Profile histograms ROOT: TProfile Class Reference を参照) は、2次元データをあるX軸の区画ごとに、Y軸のデータの平均値と標準偏差などを描画したものである。

matplotlib - Plotting profile hitstograms in python - Stack Overflow のソースコードを参考に、実際に描画してみる

import scipy.stats
import numpy as np
import matplotlib.pyplot as plt

x = np.random.rand(10000)
y = x + scipy.stats.norm(0, 0.2).rvs(10000)

plt.scatter(x, y, marker=".", s=1)
plt.ylim(-0.5, 1.5)
plt.show()

means_result = scipy.stats.binned_statistic(x, [y, y**2], bins=50, range=[0,1], statistic='mean') #binごとの統計を計算

means, means2 = means_result.statistic # yとy**2の平均値
standard_deviations = np.sqrt(means2 - means**2) # 標準偏差を計算

bin_edges = means_result.bin_edges # binの両端
bin_centers = (bin_edges[:-1] + bin_edges[1:])/2.0 # binの中心を計算

plt.errorbar(x=bin_centers, y=means, yerr=standard_deviations, linestyle='none', marker='.')
plt.ylim(-0.5, 1.5)
plt.show()

絵を見るとプロファイルヒストグラムが何かは一目瞭然だろう。

検出効率のような2値を扱う場合も書いてみる。0以上は検出できた=1、0未満は検出できなかった=0とする。

import scipy.stats
import numpy as np
import matplotlib.pyplot as plt

x = np.random.rand(1000)
y = x + scipy.stats.norm(0, 0.2).rvs(1000)
y = np.array([(1 if yi>=0 else 0) for yi in y]) #何かを検出したかどうかを2値で表す
one = np.ones(1000)

plt.scatter(x, y, marker=".", s=1)
plt.ylim(-0.5, 1.5)
plt.show()

means_result = scipy.stats.binned_statistic(x, y, bins=50, range=[0,1], statistic='mean')
means = means_result.statistic #検出効率
ones = scipy.stats.binned_statistic(x, one, bins=50, range=[0,1], statistic='sum').statistic #binに入った数

standard_deviations = np.sqrt(means * (1-means) / ones) # 二項分布の誤差

bin_edges = means_result.bin_edges # binの両端
bin_centers = (bin_edges[:-1] + bin_edges[1:])/2.0 # binの中心を計算

plt.errorbar(x=bin_centers, y=means, yerr=standard_deviations, linestyle='none', marker='.')
plt.ylim(-0.5, 1.5)
plt.show()

できた。