DesignAnalysis/BB-energy-coincidences.py

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