#!/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")