%Montecarlo, Eficiencia Tragaldabas.

%Rs Radio de la zona activa de la muestra
%Esp Espesor del cilindro que contiene la muestra
%Rd radio de la ventana del geiger
%n numero de rayos gamma detectados.
%c numero de rayos gamma que atraviesan ambos detectores.
clear
close all

d=1.2; %Distancia entre ambos detectores.
n=10000; %Número de rayos gamma detectados.
c=0; %Número de detecciones en el detector  7de abajo
a=0; %Contador.
lx=1.2; %Fondo del detector (ambos son iguales)
ly=1.5; %Ancho del detector (ambos son iguales)

%Punto aleatorio del detector. Rayos incidentes con ángulo aleatorio.


for i=1:n

    x0(i)=rand*lx; %Genero un punto aleatorio de detección en el primer detector (eje x)
    y0(i)=rand*ly; %Genero un punto aleatorio de detección en el primer detector (eje y)
    theta(i)=pi*rand; %Ángulo de incidencia colatitud aleatorio.
    phi(i)=2*pi*rand; %Ángulo de incidencia azimutal aleatorio.
    l(i)=d*tan(theta(i)); %Distancia absoluta horizontal entre el punto de detección en el primer detector y en el segundo.

    x1(i)=x0(i)+l(i).*cos(phi(i)); %Posición de la detección en el segundo detector (eje x)
    y1(i)=y0(i)+l(i).*sin(phi(i)); %Posición de la detección en el segundo detector (eje y)

    %Está x1,y1 dentro del segundo detector?

    if (x1(i)>0) & (x1(i)<lx) & (y1(i)>0) & (y1(i)<ly) %Si se cumple, hay detección!
        c=c+1;
    end


end
    eficiencia=c/n

%Distribución de las detecciones en el primer detector.
figure(1)
plot(x0,y0,'.')
axis([0,lx,0,ly])

%Distribución de las detecciones en el segundo detector.
figure(2)
plot(x1,y1,'.')
axis([0,lx,0,ly])

%Comparación de las que entran en el detector y las que caen fuera
figure(3)
if (x1>=0) & (x1<=lx) & (y1>=0) & (y1<=ly)
plot(x1,y1,'.')
axis([-3*lx,lx+3*lx,-3*ly,ly+3*ly])
else
plot(x1,y1,'.')
axis([-3*lx,lx+3*lx,-3*ly,ly+3*ly])
end