Eficiencia_de_sensores_Monte_Carlo.py

import numpy as np 
import matplotlib.pyplot as plt
from matplotlib import animation
from mpl_toolkits.mplot3d import Axes3D
from scipy.interpolate import interp1d
import pandas as pd
# ================================
#     PARÁMETROS GLOBALES
# ================================

plt.rcParams['font.family'] = 'serif'

# Tamaños de fuente (los que usaste en tu ejemplo)
plt.rcParams['axes.labelsize'] = 30     # xlabel, ylabel
plt.rcParams['axes.titlesize'] = 30     # title
plt.rcParams['xtick.labelsize'] = 24
plt.rcParams['ytick.labelsize'] = 24
plt.rcParams['legend.fontsize'] = 26
plt.rcParams['font.size'] = 28

# Grosor de bordes del gráfico ("spines")
plt.rcParams['axes.linewidth'] = 2.5

# Grosor y longitud de los ticks
plt.rcParams['xtick.major.width'] = 3
plt.rcParams['ytick.major.width'] = 3
plt.rcParams['xtick.major.size'] = 14
plt.rcParams['ytick.major.size'] = 14

plt.rcParams['xtick.minor.width'] = 1.75
plt.rcParams['ytick.minor.width'] = 1.75
plt.rcParams['xtick.minor.size'] = 10
plt.rcParams['ytick.minor.size'] = 10

# Grilla por defecto
plt.rcParams['grid.linestyle'] = '-'
plt.rcParams['grid.linewidth'] = 1.0
plt.rcParams['grid.alpha'] = 0.7

plt.rcParams['figure.figsize'] = (16, 12)

# ================================
#  FIN DE PARÁMETROS GLOBALES
# ================================

"""-----------------Dimensiones del sensor------------------"""
Espesor_skipper_micrones = 675  # Espesor del silicio en micrones
Altura_skipper_micrones = 92160 # Altura del silicio en micrones
Ancho_skipper_micrones = 15060 # Ancho del silicio en micrones
Energias_mev = np.linspace(0.1, 10, 100)  # Energías de los fotones en MeV
"""----------------------------Creación de fotones------------------------"""


def crear_fotones(num_fotones):
    """
    Crea fotones con posición aleatoria dentro del volumen del sensor
    y energía entre 0.1 y 10 MeV.
    """

    Lista_fotones = []
    for _ in range(num_fotones):
        x = np.random.uniform(0, Ancho_skipper_micrones)
        y = np.random.uniform(0, Altura_skipper_micrones)
        z = np.random.uniform(0, Espesor_skipper_micrones)
        x = Ancho_skipper_micrones/2
        y = Altura_skipper_micrones/2
        z = 0
        energia = np.random.choice(Energias_mev)  # Usar valores discretos
        Lista_fotones.append(np.array([x, y, z, energia]))
    return Lista_fotones

"""----------------Dirección aleatoria esférica--------------------------"""
def direccion_de_propagacion():
    """
    Devuelve una dirección (dx, dy, dz) aleatoria sobre una esfera unitaria.
    """
    u = np.random.uniform(0, 1)
    v = np.random.uniform(0, 1)
    theta = np.arccos(1 - 2 * u)  # theta en [0, pi]
    phi = 2 * np.pi * v  # phi en [0, 2
 

    dx = np.sin(theta) * np.cos(phi)
    dy = np.sin(theta) * np.sin(phi)
    dz = np.cos(theta)
    return np.array([dx, dy, dz])

"""------------Fin de Pruebas 1/Pruebas 2--------------------------"""
# Generar muchas direcciones
N = 5000
direcciones = np.array([direccion_de_propagacion() for _ in range(N)])
dx, dy, dz = direcciones.T

# Graficar
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(dx, dy, dz, s=2, alpha=0.5)

# Esfera de referencia
ax.set_xlim([-1, 1])
ax.set_ylim([-1, 1])
ax.set_zlim([-1, 1])
ax.set_xlabel('dx')
ax.set_ylabel('dy')
ax.set_zlabel('dz')
ax.set_title("Distribución de direcciones sobre la esfera unitaria")
#plt.show()



cos_thetas = []
phis = []

for _ in range(N):
    u= np.random.uniform(0, 1)
    v = np.random.uniform(0, 1)

    theta = np.arccos(1 - 2 * u)  # theta en [0, pi]
    phi = 2 * np.pi * v  # phi en [0, 2
    cos_theta= np.cos(theta)  # coseno del ángulo polar
    cos_thetas.append(cos_theta)
    phis.append(phi)

plt.figure(figsize=(12, 4))

plt.subplot(1, 2, 1)
plt.hist(cos_thetas, bins=50, density=True, color='skyblue', edgecolor='k')
plt.title("Distribución de cos(θ)")
plt.xlabel("cos(θ)"); plt.ylabel("Densidad")

plt.subplot(1, 2, 2)
plt.hist(phis, bins=50, density=True, color='salmon', edgecolor='k')
plt.title("Distribución de φ")
plt.xlabel("φ (radianes)"); plt.ylabel("Densidad")

plt.tight_layout()
plt.show()
"""------------FIN Pruebas 2--------------------------"""

"""------------------ACA VA LA SEECION DE ATENUACION-------------------"""
"""----------------Simulación de trayectoria--------------------------"""
def simular_propagacion_foton(x0, y0, z0, direccion, paso_micrones=10, max_pasos=1000):
    """
    Simula el trayecto del fotón dentro del volumen del sensor.
    Se detiene si el fotón sale del volumen.
    """
    trayectoria = []
    x, y, z = x0, y0, z0
    dx, dy, dz = direccion #/ np.linalg.norm(direccion)  # normalizar

    for _ in range(max_pasos):
        trayectoria.append((x, y, z))
        x += paso_micrones * dx
        y += paso_micrones * dy
        z += paso_micrones * dz

        if not (0 <= x <= Ancho_skipper_micrones and
                0 <= y <= Altura_skipper_micrones and
                0 <= z <= Espesor_skipper_micrones):
            break

    return np.array(trayectoria)

"""----------------Animación 3D--------------------------"""
def animar_trayectorias(lista_trayectorias):
    fig = plt.figure(figsize=(8, 6))
    ax = fig.add_subplot(111, projection='3d')

    ax.set_xlim(0, Ancho_skipper_micrones)
    ax.set_ylim(0, Altura_skipper_micrones)
    ax.set_zlim(0, Espesor_skipper_micrones)
    ax.set_xlabel("x [μm]")
    ax.set_ylabel("y [μm]")
    ax.set_zlabel("z [μm]")
    ax.set_title("Propagación de fotones en el sensor")

    # Crear puntos y líneas para cada fotón
    puntos = []
    caminos = []

    colores = ['red', 'green', 'blue', 'purple', 'orange']  # más si querés más fotones

    for i, trayectoria in enumerate(lista_trayectorias):
        p, = ax.plot([], [], [], 'o', color=colores[i % len(colores)])
        l, = ax.plot([], [], [], '-', color=colores[i % len(colores)], alpha=0.6)
        puntos.append(p)
        caminos.append(l)

    max_len = max(len(t) for t in lista_trayectorias)

    def init():
        for p, l in zip(puntos, caminos):
            p.set_data([], [])
            p.set_3d_properties([])
            l.set_data([], [])
            l.set_3d_properties([])
        return puntos + caminos

    def update(i):
        for idx, trayectoria in enumerate(lista_trayectorias):
            if i < len(trayectoria):
                x, y, z = trayectoria[:i+1].T
                puntos[idx].set_data([x[-1]], [y[-1]])
                puntos[idx].set_3d_properties([z[-1]])
                caminos[idx].set_data(x, y)
                caminos[idx].set_3d_properties(z)
        return puntos + caminos

    ani = animation.FuncAnimation(
        fig, update, frames=max_len,
        init_func=init, blit=True, interval=50, repeat=False
    )

    plt.show()


"""---------------------------Pruebas 3--------------------------"""

# --------------------- Carga de datos de atenuación ---------------------
df = pd.read_csv("Si_attenuation_data.txt", sep="\t")  # Asegurate de que el archivo exista
densidad_silicio = 2.33  # g/cm³


# --------------------- Interpolación de lambda(E) ---------------------
mu_m_over_rho = np.interp(Energias_mev,
                          df["Photon Energy [MeV]"],
                          df["Total Attenuation (With Coherent) [cm^2/g]"])
mu = mu_m_over_rho * densidad_silicio      # [1/cm]
lambda_cm = 1 / mu                         # [cm]
lambda_um = lambda_cm * 1e4                # [micrones]

lambda_interp = interp1d(Energias_mev, lambda_um, kind='linear', fill_value="extrapolate")

def lambda_de_energia(energia_mev):
    return lambda_interp(energia_mev)

def generar_camino_libre(energia_mev):
    return np.random.exponential(scale=lambda_de_energia(energia_mev))

# --------------------- Generador de fotones ---------------------
def crear_fotones(num_fotones):
    Lista_fotones = []
    for _ in range(num_fotones):
        x = np.random.uniform(0, Ancho_skipper_micrones)
        y = np.random.uniform(0, Altura_skipper_micrones)
        z = np.random.uniform(0, Espesor_skipper_micrones)
        x = Ancho_skipper_micrones/2
        y = Altura_skipper_micrones/2
        z = 0
        energia = np.random.choice(Energias_mev)  # Usar valores discretos
        Lista_fotones.append(np.array([x, y, z, energia]))
    return Lista_fotones

# --------------------- Dirección aleatoria en la esfera ---------------------
def direccion_de_propagacion():
    """
    Devuelve una dirección (dx, dy, dz) aleatoria sobre una esfera unitaria.
    """
    u = np.random.uniform(0, 1)
    v = np.random.uniform(0, 1)
    theta = np.arccos(1 - 2 * u)  # theta en [0, pi]
    phi =  np.pi * v  # phi en [0, 2]
 

    dx = np.sin(theta) * np.cos(phi)
    dy = np.sin(theta) * np.sin(phi)
    dz = np.cos(theta)
    #dx, dy, dz = 0, 0, 1
    return np.array([dx, dy, dz])

# --------------------- Simulación de trayectorias ---------------------
def simular_propagacion_con_atenuacion(x0, y0, z0, direccion, energia_mev, paso=10):
    lambda_um = lambda_de_energia(energia_mev)
    distancia_total = np.random.exponential(scale=lambda_um)

    x, y, z = x0, y0, z0
    dx, dy, dz = direccion / np.linalg.norm(direccion)
    distancia = 0

    trayectoria = []
    while distancia < distancia_total:
        if not (0 <= x <= Ancho_skipper_micrones and
                0 <= y <= Altura_skipper_micrones and
                0 <= z <= Espesor_skipper_micrones):
            return np.array(trayectoria), False  # Escapó sin interactuar
        trayectoria.append([x, y, z])
        x += paso * dx
        y += paso * dy
        z += paso * dz
        distancia += paso

    return np.array(trayectoria), True  # Interactuó dentro

# --------------------- Visualización de trayectorias ---------------------
def graficar_trayectorias(lista_trayectorias):
    from mpl_toolkits.mplot3d import Axes3D
    fig = plt.figure(figsize=(8, 6))
    ax = fig.add_subplot(111, projection='3d')
    ax.set_xlim(0, Ancho_skipper_micrones)
    ax.set_ylim(0, Altura_skipper_micrones)
    ax.set_zlim(0, Espesor_skipper_micrones)
    ax.set_xlabel("x [μm]")
    ax.set_ylabel("y [μm]")
    ax.set_zlabel("z [μm]")
    ax.set_title("Trayectorias de fotones en el sensor")

    colores = ['red', 'green', 'blue', 'orange', 'purple']
    for i, trayectoria in enumerate(lista_trayectorias):
        if trayectoria.size == 0:
            continue
        x, y, z = trayectoria.T
        ax.plot(x, y, z, color=colores[i % len(colores)], label=f"Fotón {i+1}")

    ax.legend()
    plt.tight_layout()
    plt.show()


"""----------------Ejecutar simulación--------------------------"""
if __name__ == "__main__":
    # Crear dos fotones
    fotones = crear_fotones(20)
    trayectorias = []

    for foton in fotones:
        x0, y0, z0, E = foton
        direccion = direccion_de_propagacion()
        tray = simular_propagacion_foton(x0, y0, z0, direccion)
        trayectorias.append(tray)

    animar_trayectorias(trayectorias)
# --------------------- Código principal ---------------------
if __name__ == "__main__":

    num_fotones = 1000
    fotones = crear_fotones(num_fotones)
    eficiencias = []

    for i in range ((len(Energias_mev))):
        trayectorias = []
        interacciones = 0
        escapes = 0

        for f in fotones:
            x0, y0, z0, E = f
            direccion = direccion_de_propagacion()
            trayectoria, interactuo = simular_propagacion_con_atenuacion(x0, y0, z0, direccion, Energias_mev[i])
            trayectorias.append(trayectoria)
            if interactuo:
                interacciones += 1
            else:
                escapes += 1
        eficiencia= interacciones / num_fotones
        eficiencias.append(eficiencia)
        print(f"Energia de los fotones: {Energias_mev[i]} MeV")
        print(f"Fotones totales:      {num_fotones}")
        print(f"Interacciones:        {interacciones}")
        print(f"Escaparon sin inter.: {escapes}")
        print("-" * 40)

        # Mostrar solo los primeros trayectos en el gráfico
        
        if i==len(Energias_mev): animar_trayectorias(trayectorias[:5])
    # Graficar eficiencia vs energía
    plt.figure(figsize=(8, 5))
    plt.plot(Energias_mev, np.array(eficiencias), marker='o', color='blue', label="Eficiencia")
    plt.xlabel("Energía del fotón [MeV]")
    plt.ylabel("Eficiencia de interacción")
    plt.title("Eficiencia vs Energía de fotones en silicio")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()