433 lines
No EOL
14 KiB
Python
433 lines
No EOL
14 KiB
Python
#!/usr/bin/env python3
|
||
# -*- coding: utf-8 -*-
|
||
|
||
"""
|
||
Bethe-Bloch energy loss study using RANGE method with fine grids and thresholds.
|
||
|
||
Setups:
|
||
- CHAOS (with/without Al)
|
||
- Multi-BGO (no silicon)
|
||
- Silicon-only penetration thresholds
|
||
"""
|
||
|
||
import numpy as np
|
||
from scipy.constants import *
|
||
|
||
# =========================================================
|
||
# GLOBAL ENERGY GRID SETTING
|
||
# =========================================================
|
||
|
||
E_MIN = -1 # log10(0.1 MeV)
|
||
E_MAX = 4
|
||
|
||
# =========================================================
|
||
# CONSTANTS
|
||
# =========================================================
|
||
|
||
e0 = epsilon_0
|
||
me = m_e
|
||
|
||
def J_to_MeV(J):
|
||
return J / (e * 1e6)
|
||
|
||
def MeV_to_J(MeV):
|
||
return MeV * (e * 1e6)
|
||
|
||
# =========================================================
|
||
# PARTICLES
|
||
# =========================================================
|
||
|
||
particles = {
|
||
'Proton': {'m0': m_p, 'z': 1},
|
||
'Helium': {'m0': 6.644657e-27, 'z': 2},
|
||
'Muon': {'m0': 1.883531627e-28, 'z': 1}
|
||
}
|
||
|
||
# =========================================================
|
||
# MATERIALS
|
||
# =========================================================
|
||
|
||
materials = {
|
||
'Si': {
|
||
'ne': (2.336e3 / (28.085e-3)) * 14.0 * N_A,
|
||
'I': 173 * e
|
||
},
|
||
'BGO': {
|
||
'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),
|
||
'I': (4 * 250 * e + 3 * 290 * e + 12 * 150 * e) / (4 + 3 + 12)
|
||
},
|
||
'Al': {
|
||
'ne': (2.70e3 / (26.98e-3)) * 13 * N_A,
|
||
'I': 166 * e
|
||
},
|
||
'BC430': {
|
||
'ne': 3.36e29 / 1e3,
|
||
'I': 60. * e # typisch für H/C‑Plastik
|
||
}
|
||
}
|
||
|
||
# =========================================================
|
||
# RELATIVITY
|
||
# =========================================================
|
||
|
||
def rel_speed(Ekin, m0):
|
||
gamma_ = Ekin / (m0 * c**2) + 1
|
||
beta2 = 1 - 1 / gamma_**2
|
||
beta2 = np.clip(beta2, 1e-12, 1.0)
|
||
return c * np.sqrt(beta2)
|
||
|
||
def gamma(v):
|
||
beta2 = v**2 / c**2
|
||
beta2 = np.clip(beta2, 0, 1 - 1e-12)
|
||
return 1 / np.sqrt(1 - beta2)
|
||
|
||
# =========================================================
|
||
# BETHE-BLOCH
|
||
# =========================================================
|
||
|
||
def bethebloch(Ekin, m0, z, ne, I):
|
||
Ekin = np.maximum(Ekin, 1e3 * e)
|
||
v = rel_speed(Ekin, m0)
|
||
v = np.maximum(v, 1e-6 * c)
|
||
g = gamma(v)
|
||
beta = v**2 / c**2
|
||
C = ne * z**2 * e**4 / (4 * np.pi * me * v**2 * e0**2)
|
||
arg = 2 * g**2 * me * v**2 / I
|
||
arg = np.maximum(arg, 1.0001)
|
||
D = np.log(arg)
|
||
return C * (D - beta)
|
||
|
||
# =========================================================
|
||
# RANGE METHOD
|
||
# =========================================================
|
||
|
||
def build_range_table(particle, material):
|
||
p = particles[particle]
|
||
m = materials[material]
|
||
|
||
E_grid = np.logspace(E_MIN, E_MAX, 10000)
|
||
E_J = MeV_to_J(E_grid)
|
||
dEdx = bethebloch(E_J, p['m0'], p['z'], m['ne'], m['I'])
|
||
|
||
R = np.zeros_like(E_grid)
|
||
for i in range(1, len(E_grid)):
|
||
dE = E_J[i] - E_J[i-1]
|
||
R[i] = R[i-1] + dE / dEdx[i]
|
||
|
||
return E_grid, R
|
||
|
||
def propagate_range(E0, layers, tables):
|
||
E = E0
|
||
for mat, thickness in layers:
|
||
if E <= 0:
|
||
return 0
|
||
E_grid, R = tables[mat]
|
||
R0 = np.interp(E, E_grid, R)
|
||
if R0 < thickness:
|
||
return 0
|
||
R_rest = R0 - thickness
|
||
E = np.interp(R_rest, R, E_grid)
|
||
return E
|
||
|
||
# =========================================================
|
||
# THRESHOLDS (BGO)
|
||
# =========================================================
|
||
|
||
def find_thresholds(particle, layers):
|
||
tables = {mat: build_range_table(particle, mat) for mat, _ in layers}
|
||
E0_grid = np.logspace(E_MIN, E_MAX, 10000)
|
||
|
||
Ec_list = []
|
||
Ebgo_list = []
|
||
|
||
for E0 in E0_grid:
|
||
E = E0
|
||
for mat, d in layers[:-1]:
|
||
E = propagate_range(E, [(mat, d)], tables)
|
||
Ec = E
|
||
Ebgo = propagate_range(Ec, [layers[-1]], tables)
|
||
|
||
Ec_list.append(Ec)
|
||
Ebgo_list.append(Ebgo)
|
||
|
||
Ec_list = np.array(Ec_list)
|
||
Ebgo_list = np.array(Ebgo_list)
|
||
|
||
E_reach0 = E0_grid[np.argmax(Ec_list > 0)]
|
||
E_reach4 = E0_grid[np.argmax(Ec_list > 4)]
|
||
E_transition = E0_grid[np.argmax(Ebgo_list > 0)]
|
||
|
||
return E_reach0, E_reach4, E_transition
|
||
|
||
# =========================================================
|
||
# THRESHOLDS (Detector C)
|
||
# =========================================================
|
||
|
||
def find_thresholds_C(particle, layers):
|
||
tables = {mat: build_range_table(particle, mat) for mat, _ in layers}
|
||
|
||
# Nur bis zur Al-Schicht vor Detector C
|
||
layers_to_C = layers[:2] # Si (A) + Al
|
||
|
||
E0_grid = np.logspace(E_MIN, E_MAX, 10000)
|
||
Ec_list = []
|
||
|
||
for E0 in E0_grid:
|
||
E = propagate_range(E0, layers_to_C, tables)
|
||
Ec_list.append(E)
|
||
|
||
Ec_list = np.array(Ec_list)
|
||
E_reach0 = E0_grid[np.argmax(Ec_list > 0)]
|
||
E_reach34keV = E0_grid[np.argmax(Ec_list > 0.034)]
|
||
|
||
return E_reach0, E_reach34keV
|
||
|
||
# =========================================================
|
||
# SETUPS
|
||
# =========================================================
|
||
|
||
CHAOS = [
|
||
('Si', 300e-6), # Detector A
|
||
('Al', 2e-3), # Housing
|
||
('Si', 300e-6), # Detector C
|
||
('BGO', 2e-2)
|
||
]
|
||
|
||
CHAOS_Al = [
|
||
('Al', 5.5e-3), # Pressure housing
|
||
('Si', 300e-6), # Detector A
|
||
('Al', 2e-3), # Housing
|
||
('Si', 300e-6), # Detector C
|
||
('BGO', 2e-2)
|
||
]
|
||
|
||
# # =========================================================
|
||
# # MAIN CHAOS COMPARISON
|
||
# # =========================================================
|
||
|
||
# for setup_name, layers in [
|
||
# ("CHAOS with Cherenkov housing (2mm), without pressure housing", CHAOS),
|
||
# ("CHAOS including Cherenkov housing (2mm) and pressure housing (5.5 mm)", CHAOS_Al)
|
||
# ]:
|
||
|
||
# print(f"\n==============================")
|
||
# print(setup_name)
|
||
# print("==============================")
|
||
|
||
# for p in particles:
|
||
# # Threshold nach Al-Schicht für Detector C
|
||
# Ec0, Ec34 = find_thresholds_C(p, layers)
|
||
# # Thresholds für BGO und Transition
|
||
# E0_0, E0_4, E_transition = find_thresholds(p, layers)
|
||
|
||
# print(f"\n--- {p} ---")
|
||
# print(f"Minimum E0 to reach Detector C : {Ec0:.3f} MeV (E_rest>0), {Ec34:.3f} MeV (E_rest>34 keV)")
|
||
# print(f"Minimum E0 to reach BGO : {E0_0:.3f} MeV (E_rest>0), {E0_4:.3f} MeV (E_rest>4 MeV)")
|
||
# print(f"Transition (stop → penetrate) : {E_transition:.3f} MeV")
|
||
|
||
|
||
# =========================================================
|
||
# MULTI-BGO (NO SILICON)
|
||
# =========================================================
|
||
|
||
print("\n========================================")
|
||
print("MULTI-BGO (no silicon)")
|
||
print("========================================")
|
||
|
||
BGO_thicknesses = [2e-2, 4e-2, 6e-2]
|
||
|
||
for p in particles:
|
||
print(f"\nParticle: {p}")
|
||
tables = {'BGO': build_range_table(p, 'BGO')}
|
||
for d in BGO_thicknesses:
|
||
E0_grid = np.logspace(E_MIN, E_MAX, 5000)
|
||
for E0 in E0_grid:
|
||
E_rest = propagate_range(E0, [('BGO', d)], tables)
|
||
if E_rest > 0:
|
||
print(f"BGO thickness = {d*100:.0f} cm → Transition E0 = {E0:.3f} MeV")
|
||
break
|
||
|
||
# =========================================================
|
||
# SILICON ONLY
|
||
# =========================================================
|
||
|
||
print("\n========================================")
|
||
print("SILICON ONLY")
|
||
print("========================================")
|
||
|
||
Si_configs = {
|
||
"300 um Si + 300 um Si": 300e-6,
|
||
"500 um Si + 500 um Si": 500e-6
|
||
}
|
||
|
||
for label, thickness in Si_configs.items():
|
||
print(f"\n--- {label} ---")
|
||
for p in particles:
|
||
tables = {'Si': build_range_table(p, 'Si')}
|
||
layers = [('Si', thickness), ('Si', thickness)]
|
||
E0_grid = np.logspace(E_MIN, E_MAX, 5000)
|
||
for E0 in E0_grid:
|
||
E_rest = propagate_range(E0, layers, tables)
|
||
if E_rest > 0:
|
||
print(f"{p}: Minimum E0 to penetrate Si = {E0:.3f} MeV")
|
||
break
|
||
|
||
|
||
# =========================================================
|
||
# MAIN ProHEPaM ANALYSIS
|
||
# =========================================================
|
||
|
||
# FEINERES GRID NUR FÜR ProHEPaM (nicht alle anderen Setups)
|
||
E0_grid = np.logspace(E_MIN, E_MAX, 50000) # statt 5000
|
||
|
||
print("\n========================================")
|
||
print("ProHEPaM setup")
|
||
print("========================================")
|
||
|
||
|
||
ProHEPaM = [
|
||
('Si', 2*500e-6), # Detector S1
|
||
('Al', 2e-3), # Housing
|
||
('Si', 2*500e-6), # Detector S2
|
||
('BGO', 2e-2), # Detector B1
|
||
('Si', 2*500e-6), # Detector S3
|
||
('BGO', 2e-2), # Detector B2
|
||
('Si', 2*500e-6), # Detector S4
|
||
('BGO', 2e-2), # Detector B3
|
||
('Si', 2*500e-6), # Detector S5
|
||
('BC430', 2e-2) # Detector P
|
||
]
|
||
|
||
|
||
def find_E_to_reach(particle, layers, target_index, tables, E0_grid=None):
|
||
"""Return minimum E0 (MeV) so that res at entrance of layer[target_index] > 0."""
|
||
if E0_grid is None:
|
||
E0_grid = np.logspace(E_MIN, E_MAX, 5000)
|
||
|
||
for E0 in E0_grid:
|
||
E = E0
|
||
for i, (mat, d) in enumerate(layers):
|
||
if i > target_index:
|
||
break
|
||
E = propagate_range(E, [(mat, d)], tables)
|
||
if E <= 0:
|
||
break
|
||
if E > 0:
|
||
return E0
|
||
return 0.0
|
||
|
||
|
||
def find_transition_E_to_BGO(particle, layers, bgo_index, tables, E0_grid=None):
|
||
"""Return E0 such that the particle just penetrates the BGO block at bgo_index."""
|
||
if E0_grid is None:
|
||
E0_grid = np.logspace(E_MIN, E_MAX, 5000)
|
||
|
||
for i, E0 in enumerate(E0_grid):
|
||
E = E0
|
||
for i_mat, (mat, d) in enumerate(layers):
|
||
if i_mat >= bgo_index:
|
||
break
|
||
E = propagate_range(E, [(mat, d)], tables)
|
||
if E <= 0:
|
||
break
|
||
if E <= 0:
|
||
continue
|
||
# Energie am Eingang des BGO‑Blocks
|
||
E_in_BGO = E
|
||
# Durchqueren des BGO‑Blocks
|
||
E_out_BGO = propagate_range(E_in_BGO, [('BGO', layers[bgo_index][1])], tables)
|
||
if E_out_BGO > 0:
|
||
return E0
|
||
return 0.0
|
||
|
||
def find_thresholds_BX(particle, layers, s_before_index, bgo_index, tables, E0_grid):
|
||
"""
|
||
Für ProHEPaM: B1, B2, B3.
|
||
- s_before_index: Index des Si‑Detektors direkt vor dem BGO (z.B. S2 vor B1).
|
||
- bgo_index: Index des BGO‑Blocks (B1/B2/B3).
|
||
"""
|
||
E0_reach_BX = 0.0 # Min E0 so dass BGO erreicht wird (E nach S > 0)
|
||
E0_detect_BX = 0.0 # Min E0 so dass E nach S > 4 MeV
|
||
E0_trans_BX = 0.0 # Min E0 so dass BGO durchdrungen wird (E_out > 0)
|
||
|
||
for E0 in E0_grid:
|
||
E = E0
|
||
for i_mat, (mat, d) in enumerate(layers):
|
||
if i_mat > bgo_index:
|
||
break
|
||
E = propagate_range(E, [(mat, d)], tables)
|
||
if E <= 0:
|
||
break
|
||
|
||
if E <= 0:
|
||
continue
|
||
# Energie nach S vor dem BGO (bei s_before_index)
|
||
E_after_S = E
|
||
# Durchqueren des BGO‑Blocks
|
||
E_out_BGO = propagate_range(E_after_S, [('BGO', layers[bgo_index][1])], tables)
|
||
|
||
# Reach: BGO überhaupt noch erreicht
|
||
if E0_reach_BX == 0.0 and E_after_S > 0:
|
||
E0_reach_BX = E0
|
||
|
||
# Detection‑threshold in BX (z.B. E_rest > 4 MeV nach S)
|
||
if E0_detect_BX == 0.0 and E_after_S > 4.0:
|
||
E0_detect_BX = E0
|
||
|
||
# Transition: BGO durchdrungen
|
||
if E0_trans_BX == 0.0 and E_out_BGO > 0:
|
||
E0_trans_BX = E0
|
||
|
||
return E0_reach_BX, E0_detect_BX, E0_trans_BX
|
||
|
||
|
||
for p in particles:
|
||
print(f"\n--- {p} ---")
|
||
tables = {
|
||
'Si': build_range_table(p, 'Si'),
|
||
'Al': build_range_table(p, 'Al'),
|
||
'BGO': build_range_table(p, 'BGO'),
|
||
'BC430': build_range_table(p, 'BC430')
|
||
}
|
||
|
||
E0_grid = np.logspace(E_MIN, E_MAX, 20000)
|
||
|
||
|
||
# ==============================================
|
||
# 1. Liste: Min. E0 to reach S1, S2, B1, S3, B2, S4, B3, S5, P
|
||
# ==============================================
|
||
|
||
E0_S2 = find_E_to_reach(p, ProHEPaM, 2, tables, E0_grid) # S2 index = 2
|
||
E0_B1 = find_E_to_reach(p, ProHEPaM, 3, tables, E0_grid) # B1 index = 3
|
||
E0_S3 = find_E_to_reach(p, ProHEPaM, 4, tables, E0_grid) # S3 index = 4
|
||
E0_B2 = find_E_to_reach(p, ProHEPaM, 5, tables, E0_grid) # B2 index = 5
|
||
E0_S4 = find_E_to_reach(p, ProHEPaM, 6, tables, E0_grid) # S4 index = 6
|
||
E0_B3 = find_E_to_reach(p, ProHEPaM, 7, tables, E0_grid) # B3 index = 7
|
||
E0_S5 = find_E_to_reach(p, ProHEPaM, 8, tables, E0_grid) # S5 index = 8
|
||
E0_P = find_E_to_reach(p, ProHEPaM, 9, tables, E0_grid) # P index = 9
|
||
|
||
print("List of E0 thresholds (reach):")
|
||
print(f" Min. E0 to reach S2 (2*500 um Si + Al) : {E0_S2:.3f} MeV")
|
||
print(f" Min. E0 to reach B1 (E_after_S2>0) : {E0_B1:.3f} MeV")
|
||
print(f" Min. E0 to reach S3 (incl. B1) : {E0_S3:.3f} MeV")
|
||
print(f" Min. E0 to reach B2 (E_after_S3>0) : {E0_B2:.3f} MeV")
|
||
print(f" Min. E0 to reach S4 (incl. B2) : {E0_S4:.3f} MeV")
|
||
print(f" Min. E0 to reach B3 (E_after_S4>0) : {E0_B3:.3f} MeV")
|
||
print(f" Min. E0 to reach S5 (incl. B3) : {E0_S5:.3f} MeV")
|
||
print(f" Min. E0 to reach P (BC430) : {E0_P:.3f} MeV")
|
||
|
||
|
||
# ==============================================
|
||
# 2. Drei Transition‑Energien für B1, B2, B3
|
||
# ==============================================
|
||
|
||
_, _, E0_t_B1 = find_thresholds_BX(p, ProHEPaM, 2, 3, tables, E0_grid)
|
||
_, _, E0_t_B2 = find_thresholds_BX(p, ProHEPaM, 4, 5, tables, E0_grid)
|
||
_, _, E0_t_B3 = find_thresholds_BX(p, ProHEPaM, 6, 7, tables, E0_grid)
|
||
|
||
print("Transition energies (stop → penetrate):")
|
||
print(f" Transition E0 for B1 (stop→penetrate) : {E0_t_B1:.3f} MeV")
|
||
print(f" Transition E0 for B2 (stop→penetrate) : {E0_t_B2:.3f} MeV")
|
||
print(f" Transition E0 for B3 (stop→penetrate) : {E0_t_B3:.3f} MeV") |