318 lines
7.6 KiB
Python
318 lines
7.6 KiB
Python
|
#!/usr/bin/env python3
|
|||
|
|
# -*- coding: utf-8 -*-
|
|||
|
|
"""
|
|||
|
|
Created on Wed Mar 4 09:30:07 2026
|
|||
|
|
@author: ava
|
|||
|
|
"""
|
|||
|
|
|
|||
|
|
import numpy as np
|
|||
|
|
import matplotlib.pyplot as plt
|
|||
|
|
from scipy.constants import *
|
|||
|
|
import pandas as pd
|
|||
|
|
from pathlib import Path
|
|||
|
|
import sys
|
|||
|
|
|
|||
|
|
sys.path.append(str(Path(__file__).resolve().parents[1]))
|
|||
|
|
import plotstyle # Dein plotstyle!
|
|||
|
|
|
|||
|
|
MODE = "paper"
|
|||
|
|
plotstyle.set_style(MODE)
|
|||
|
|
|
|||
|
|
e0 = epsilon_0
|
|||
|
|
me = m_e
|
|||
|
|
|
|||
|
|
def J_to_MeV(J):
|
|||
|
|
return J / (e * 1e6)
|
|||
|
|
|
|||
|
|
def MeV_to_J(MeV):
|
|||
|
|
return MeV * (e * 1e6)
|
|||
|
|
|
|||
|
|
Na = 6.022140857e23
|
|||
|
|
u = 1e-3 / Na
|
|||
|
|
|
|||
|
|
|
|||
|
|
particles_BB = {
|
|||
|
|
'Proton': {
|
|||
|
|
'E0_min': 30, # MeV
|
|||
|
|
'E0_max': 10000, # MeV
|
|||
|
|
'm0': m_p,
|
|||
|
|
'z': 1
|
|||
|
|
},
|
|||
|
|
'Helium': {
|
|||
|
|
'E0_min': 30, # MeV
|
|||
|
|
'E0_max': 10000, # MeV
|
|||
|
|
'm0': 6.644657e-27,
|
|||
|
|
'z': 2
|
|||
|
|
},
|
|||
|
|
'Muon': {
|
|||
|
|
'E0_min': 30, # MeV
|
|||
|
|
'E0_max': 10000, # MeV
|
|||
|
|
'm0': 1.883531627e-28,
|
|||
|
|
'z': 1
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
materials = {
|
|||
|
|
'BGO': {
|
|||
|
|
# Elektronendichte in BGO (gewichteter Durchschnitt)
|
|||
|
|
'ne': (4 * (9.78e3 / (208.98e-3)) * 83 * N_A + 3 * (5.323e3 / (72.63e-3)) * 32 * N_A + 12 * (1.429e3 / (16.00e-3)) * 8 * N_A) / (4 + 3 + 12),
|
|||
|
|
# Ionisierungsenergie in BGO (gewichteter Durchschnitt)
|
|||
|
|
'I': (4 * 250 * e + 3 * 290 * e + 12 * 150 * e) / (4 + 3 + 12),
|
|||
|
|
# Dicke in m
|
|||
|
|
'd': 2e-2
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
def rel_speed(Ekin, m0):
|
|||
|
|
return c * np.sqrt(1 - 1 / ((Ekin / m0 / c / c + 1) ** 2))
|
|||
|
|
|
|||
|
|
def gamma(v):
|
|||
|
|
return 1 / np.sqrt(1 - v * v / c / c)
|
|||
|
|
|
|||
|
|
def bethebloch(Ekin, m0, z, ne, I):
|
|||
|
|
v = rel_speed(J_to_MeV(Ekin) if isinstance(Ekin, np.ndarray) else Ekin, m0)
|
|||
|
|
lorentz_gamma = gamma(v)
|
|||
|
|
C = ne * z**2 * e**4 / (4 * np.pi * me * v**2 * e0**2)
|
|||
|
|
D = np.log(2 * lorentz_gamma**2 * me * v**2 / I)
|
|||
|
|
beta = v**2 / c**2
|
|||
|
|
return C * (D - beta)
|
|||
|
|
|
|||
|
|
def Eloss(E0, m0, z, ne, I, dx, distance):
|
|||
|
|
steps = round(distance / dx)
|
|||
|
|
E = MeV_to_J(E0)
|
|||
|
|
for i in range(steps):
|
|||
|
|
dE = bethebloch(E, m0, z, ne, I) * dx
|
|||
|
|
E -= dE
|
|||
|
|
Ediff = MeV_to_J(E0) - E
|
|||
|
|
return J_to_MeV(Ediff)
|
|||
|
|
|
|||
|
|
|
|||
|
|
# -----------------------------
|
|||
|
|
# Geant4‑Konfiguration (ohne Elektron)
|
|||
|
|
# -----------------------------
|
|||
|
|
|
|||
|
|
# Nur Proton, Helium, Mu‑ (kein e-)
|
|||
|
|
PARTICLES_BASE_G4 = {
|
|||
|
|
'proton': 'G4outfiles/proton_70energies',
|
|||
|
|
'alpha': 'G4outfiles/helium_70energies',
|
|||
|
|
'mu-': 'G4outfiles/muon_70energies'
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
colors_G4 = {
|
|||
|
|
'proton': 'red',
|
|||
|
|
'alpha': 'green',
|
|||
|
|
'mu-': (0.5, 0.3, 0.8)
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
LABELS_G4 = {
|
|||
|
|
'proton': 'Proton',
|
|||
|
|
'alpha': 'Helium',
|
|||
|
|
'mu-': 'Muon'
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
LINSTYLES = {
|
|||
|
|
2: '-', # 2 cm
|
|||
|
|
4: '--', # 4 cm
|
|||
|
|
6: ':', # 6 cm
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
THICKNESSES = [2, 4, 6]
|
|||
|
|
|
|||
|
|
# -----------------------------
|
|||
|
|
# Funktion: GESAMTER Energieverlust (ALLE Zeilen!) – nur für G4
|
|||
|
|
# -----------------------------
|
|||
|
|
|
|||
|
|
def load_particle_data(filename, part_name, label):
|
|||
|
|
try:
|
|||
|
|
df = pd.read_csv(filename, sep='\t')
|
|||
|
|
df['z_cm'] = df['z'] / 10.0
|
|||
|
|
df_part = df[df['part'] == part_name].copy()
|
|||
|
|
|
|||
|
|
if len(df_part) == 0:
|
|||
|
|
print(f" → {label}: KEINE {part_name}-Hits gefunden!")
|
|||
|
|
return None
|
|||
|
|
|
|||
|
|
total_eloss = df_part.groupby(['primaryE', 'event'])['edep'].sum().reset_index()
|
|||
|
|
total_eloss.columns = ['E0', 'event', 'dE_total']
|
|||
|
|
|
|||
|
|
means = total_eloss.groupby('E0')['dE_total'].agg(['mean','std']).reset_index()
|
|||
|
|
means.columns = ['E0', 'dE_mean', 'dE_std']
|
|||
|
|
return means
|
|||
|
|
|
|||
|
|
except FileNotFoundError:
|
|||
|
|
print(f" → Datei '{filename}' nicht gefunden!")
|
|||
|
|
return None
|
|||
|
|
except Exception as e:
|
|||
|
|
print(f" → Fehler bei {label}: {e}")
|
|||
|
|
return None
|
|||
|
|
|
|||
|
|
|
|||
|
|
# -----------------------------
|
|||
|
|
# 9 Dateien laden: Proton/Helium/Mu‑, 2/4/6 cm (ohne e-)
|
|||
|
|
# -----------------------------
|
|||
|
|
|
|||
|
|
all_data = {} # (part, thickness) -> DataFrame
|
|||
|
|
|
|||
|
|
for part_g4 in PARTICLES_BASE_G4.keys():
|
|||
|
|
label = LABELS_G4[part_g4]
|
|||
|
|
base = PARTICLES_BASE_G4[part_g4]
|
|||
|
|
|
|||
|
|
for thickness in THICKNESSES:
|
|||
|
|
thickness_label = f"{thickness} cm"
|
|||
|
|
filename = f"{base}_{thickness}cm_0.hits"
|
|||
|
|
print(f"\n=== {label} ({part_g4}), {thickness_label} ===")
|
|||
|
|
print(f" Datei: {filename}")
|
|||
|
|
|
|||
|
|
df = load_particle_data(filename, part_g4, label)
|
|||
|
|
if df is None:
|
|||
|
|
print(f" → {label} / {thickness_label}: KEINE DATEN")
|
|||
|
|
else:
|
|||
|
|
print(f" → {label} / {thickness_label}: {len(df)} Energiepunkte")
|
|||
|
|
all_data[(part_g4, thickness)] = df
|
|||
|
|
|
|||
|
|
|
|||
|
|
# -----------------------------
|
|||
|
|
# PLOT: Bethe‑Bloch + Geant4 in einem Plot
|
|||
|
|
# -----------------------------
|
|||
|
|
|
|||
|
|
plt.figure()
|
|||
|
|
|
|||
|
|
material = materials['BGO']
|
|||
|
|
dx = 1e-5
|
|||
|
|
|
|||
|
|
distances = {
|
|||
|
|
'2 cm': 0.02,
|
|||
|
|
'4 cm': 0.04,
|
|||
|
|
'6 cm': 0.06
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
# Farben und Linestyles für Bethe‑Bloch (nur Proton/Helium/Mu‑, wie in BB‑Skript)
|
|||
|
|
colors_BB = {
|
|||
|
|
'Proton': 'red',
|
|||
|
|
'Helium': 'green',
|
|||
|
|
'Muon': (0.5, 0.3, 0.8)
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
linestyles_BB = {
|
|||
|
|
'2 cm': '-',
|
|||
|
|
'4 cm': '--',
|
|||
|
|
'6 cm': ':'
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
# 1) Bethe‑Bloch‑Kurven
|
|||
|
|
for name, particle_params in particles_BB.items():
|
|||
|
|
E_values = np.logspace(
|
|||
|
|
np.log10(particle_params['E0_min']),
|
|||
|
|
np.log10(particle_params['E0_max']),
|
|||
|
|
700
|
|||
|
|
)
|
|||
|
|
|
|||
|
|
for label, distance in distances.items():
|
|||
|
|
Eloss_values = np.array([
|
|||
|
|
Eloss(
|
|||
|
|
E0,
|
|||
|
|
particle_params['m0'],
|
|||
|
|
particle_params['z'],
|
|||
|
|
material['ne'],
|
|||
|
|
material['I'],
|
|||
|
|
dx,
|
|||
|
|
distance
|
|||
|
|
) for E0 in E_values
|
|||
|
|
])
|
|||
|
|
|
|||
|
|
plt.plot(
|
|||
|
|
E_values,
|
|||
|
|
Eloss_values,
|
|||
|
|
color=colors_BB[name],
|
|||
|
|
linestyle=linestyles_BB[label],
|
|||
|
|
linewidth=2,
|
|||
|
|
label=f"{name} (Bethe‑Bloch, {label})"
|
|||
|
|
)
|
|||
|
|
|
|||
|
|
|
|||
|
|
# 2) Geant4‑Kurven (Errorbars)
|
|||
|
|
for (part_g4, thickness), df in all_data.items():
|
|||
|
|
if df is None or len(df) == 0:
|
|||
|
|
continue
|
|||
|
|
|
|||
|
|
mask = df['dE_mean'] > 0
|
|||
|
|
if mask.sum() == 0:
|
|||
|
|
continue
|
|||
|
|
|
|||
|
|
color = colors_G4[part_g4]
|
|||
|
|
linestyle = LINSTYLES[thickness]
|
|||
|
|
|
|||
|
|
x = df.loc[mask, 'E0']
|
|||
|
|
y = df.loc[mask, 'dE_mean']
|
|||
|
|
yerr = df.loc[mask, 'dE_std']
|
|||
|
|
|
|||
|
|
# Mittlere Linie für Geant4
|
|||
|
|
plt.plot(
|
|||
|
|
x, y,
|
|||
|
|
color=color,
|
|||
|
|
linestyle=linestyle,
|
|||
|
|
linewidth=1,
|
|||
|
|
label=f"{LABELS_G4[part_g4]} (G4, {thickness} cm)"
|
|||
|
|
)
|
|||
|
|
|
|||
|
|
# Fehlerschlauch (50% sichtbar, leicht durchsichtig)
|
|||
|
|
plt.fill_between(
|
|||
|
|
x,
|
|||
|
|
y - yerr,
|
|||
|
|
y + yerr,
|
|||
|
|
color=color,
|
|||
|
|
linestyle=linestyle,
|
|||
|
|
alpha=0.4
|
|||
|
|
)
|
|||
|
|
|
|||
|
|
# Referenzlinie: IDEAL ΔE = E0
|
|||
|
|
x = np.linspace(30, 10000, 100)
|
|||
|
|
plt.plot(x, x, color="black", linewidth=1, linestyle='-', label=r"$\Delta E = E_{\text{kin}}$")
|
|||
|
|
|
|||
|
|
|
|||
|
|
plt.xscale('log')
|
|||
|
|
plt.yscale('log')
|
|||
|
|
plt.ylim(10, 1000)
|
|||
|
|
plt.xlim(29, 10000)
|
|||
|
|
plt.xlabel('Kinetic energy $E_{kin}$ in MeV')
|
|||
|
|
plt.ylabel('Energy loss $\\Delta E$ in MeV')
|
|||
|
|
plt.grid(True, which="both", ls="--", lw=0.5)
|
|||
|
|
|
|||
|
|
|
|||
|
|
# -----------------------------
|
|||
|
|
# ZWEI LEGENDEN
|
|||
|
|
# 1) Particles (mit BB‑ und G4‑Labels, aber nur 3 Teilchen)
|
|||
|
|
from matplotlib.lines import Line2D
|
|||
|
|
|
|||
|
|
legend_particle_handles = []
|
|||
|
|
|
|||
|
|
for name in ['Proton', 'Helium', 'Muon']:
|
|||
|
|
color = colors_BB[name]
|
|||
|
|
legend_particle_handles.append(
|
|||
|
|
Line2D([0], [0], color=color, lw=2, label=name)
|
|||
|
|
)
|
|||
|
|
|
|||
|
|
# 2) BGO thickness (2/4/6 cm, Linestyles, schwarz)
|
|||
|
|
legend_thickness_handles = []
|
|||
|
|
for label in distances.keys():
|
|||
|
|
legend_thickness_handles.append(
|
|||
|
|
Line2D([0], [0], color='black', linestyle=linestyles_BB[label], lw=2, label=label)
|
|||
|
|
)
|
|||
|
|
|
|||
|
|
|
|||
|
|
legend_particles = plt.legend(
|
|||
|
|
handles=legend_particle_handles,
|
|||
|
|
loc='upper right',
|
|||
|
|
title="Particles"
|
|||
|
|
)
|
|||
|
|
legend_thickness = plt.legend(
|
|||
|
|
handles=legend_thickness_handles,
|
|||
|
|
loc='upper left',
|
|||
|
|
title="BGO thickness"
|
|||
|
|
)
|
|||
|
|
plt.gca().add_artist(legend_particles)
|
|||
|
|
|
|||
|
|
plt.tight_layout()
|
|||
|
|
plotstyle.savefig("BGO_Eloss_BetheBloch_vs_G4_246cm", category="BB")
|
|||
|
|
plt.show()
|