In [1]:
using PyPlot
using Distributions
PyPlot.matplotlib[:rc]("text", usetex=true)
In [2]:
a₁₁ = 11.993
a₂₂ = 32.749
a₁₂ = a₁₁ + a₂₂
a₁ = a₁₁/a₁₂
a₂ = a₂₂/a₁₂
λ₁ = 47.13
λ₂ = 1136.19
println("a₁ = ",a₁)
println("a₂ = ",a₂)
In [3]:
x = linspace(10,400,50)
S(S₀,x)= S₀*((a₁*(exp(-x/λ₁))) + (a₂*(exp(-x/λ₂))))
ylabel(L"$S(1,x)$",fontsize=20,rotation=360, labelpad=20)
xlabel(L"$x$",fontsize=20)
ylim(0.3,1)
plot(x,S(1,x))
grid("on")
In [4]:
k = linspace(0,20,21)
p₁ = Poisson(2)
p₂ = Poisson(5)
p₃ = Poisson(8)
cpois₁ = cdf(p₁,k)
cpois₂ = cdf(p₂,k)
cpois₃ = cdf(p₃,k)
ylim(0,1.3)
plot(k,cpois₁,"r",label=L"$cpois_1$")
plot(k,cpois₂,"b",label=L"$cpois_2$")
plot(k,cpois₃,"g",label=L"$cpois_3$")
legend(loc=4)
xlabel(L"$k$",fontsize=20)
Out[4]:
In [5]:
k = linspace(0,20,21)
p₁ = Poisson(2)
p₂ = Poisson(5)
p₃ = Poisson(8)
dpois₁ = pdf(p₁,k)
dpois₂ = pdf(p₂,k)
dpois₃ = pdf(p₃,k)
plot(k,dpois₁,"r",label=L"$dpois_1$")
plot(k,dpois₂,"b",label=L"$dpois_2$")
plot(k,dpois₃,"g",label=L"$dpois_3$")
legend(loc=4)
xlabel(L"$k$",fontsize=20)
Out[5]:
In [20]:
#Numer of strips
N = 50
#x-direction
i = linspace(1,N,50)
#y-direction
j = linspace(1,N,50)
#Lateral side of the detector
LD = 200
#Optical fiber bending region:
LC = 50
#pixel size
Lpix = LD/N
#Vertical muon deposition in SPEs (2 MeV)
Svm = 7
#Distance from impact point to PMT1 (x-plane)
Ld1(i,j) = LD - Lpix*i
Lc1(i,j) = sqrt(LC^2 + (Lpix*j - LD/2).*(Lpix*j - LD/2))
L1(i,j) = Ld1(i,j) + Lc1(i,j)
#Distance from impact point to PMT2 (y-plane)
Ld2(i,j) = LD - Lpix*j
Lc2(i,j) = sqrt(LC^2 + (Lpix*i - LD/2).*(Lpix*i - LD/2))
L2(i,j) = Ld1(i,j) + Lc1(i,j)
# reshape(L1(i,j),N,N)
matriz1(i,j) = hcat([L1(i,j) for j=1:50]...)
matriz2(i,j) = hcat([L2(i,j) for j=1:50]...)
# println(matriz1(i,j)[1,3])
# println(matriz2(i,j)[1,3])
for i in [1:50]
for j in [1:50]
X(i,j) = Ld1(i,j) + Lc1(i,j)
Y(i,j) = Ld2(i,j) + Lc2(i,j)
Z(i,j) = X(i,j) - Y(i,j)
end
end
#Effective angle for propagation of light inside the optical fiber
ζ = 40
c = 3*10^10*cos(ζ) #cm/sec
for i in [1:50]
for j in [1:50]
Δt(i,j) = (Z(i,j))/(c*10^(-9.0))
#print(Δt(i,j))
end
end
typeof(Δt(i,j))
# for i in [1:50]
# for j in [1:50]
# matshow(Δt(i,j))
# end
# end
#Detection probability for each (i,j) pixel
# Pd1(i,j) = 1-cdf(Poisson(1),S(Svm,L1))
# Pd2(i,j) = 1-cdf(Poisson(1),S(Svm,L2))
# PD(i,j) = Pd1 .* Pd2
# for i in [1:N]
# for j in [1:N]
# #Distance from impact point to PMT1 (x-plane)
# Ld1(i,j) = LD - Lpix*i
# Lc1(i,j) = sqrt(LC^2 + (Lpix*j - LD/2).*(Lpix*j - LD/2))
# L1(i,j) = Ld1(i,j) + Lc1(i,j)
# #Distance from impact point to PMT2 (y-plane)
# Ld2(i,j) = LD - Lpix*j
# Lc2(i,j) = sqrt(LC^2 + (Lpix*i - LD/2).*(Lpix*i - LD/2))
# L2(i,j) = Ld1(i,j) + Lc1(i,j)
# #Detection probability for each (i,j) pixel
# Pd1(i,j) = 1-cdf(Poisson(1),S(Svm,L1))
# Pd2(i,j) = 1-cdf(Poisson(1),S(Svm,L2))
# PD(i,j) = Pd1 .* Pd2
# #Time of flight to PMT1 and PMT2
# #Effective angle for propagation of light inside the optical fiber
# ζ = 40
# c = 3*10^10*cos(ζ) #cm/sec
# Δt(i,j) = (L1(i,j) - L2(i,j))/(c*10^(-9.0))
# end
# end
#plot(Δt(i,j))
#THIS MAKES NO SENSE FOR THE MOMENT
Out[20]:
In [59]:
#Depth of planes [cm]
z = [30,60,150]
#Particle trayectory
T(z,θ,φ,x₀) = z/(c*cos(θ)*10^-9.0)
X(θ,φ,x₀) = z*tan(θ)*cos(φ) + x₀
Y(θ,φ,y₀) = z*tan(θ)*sin(φ) + y₀
I(θ,φ,x₀) = ceil(X(θ,φ,x₀)/Lpix)
J(θ,φ,y₀) = ceil(Y(θ,φ,y₀)/Lpix)
ζ₁(θ,φ,x₀,y₀) = (0 <= X(θ,φ,x₀) <= 200) & (0 <= X(θ,φ,x₀) <= 200) ? 1 : 0
#Arrival time of signal ar each PMT
Tx(plane,θ,φ,x₀,y₀) = L1(I(θ,φ,x₀),J(θ,φ,y₀))/(c*10^-9.0)
Ty(plane,θ,φ,x₀,y₀) = L2(I(θ,φ,x₀),J(θ,φ,y₀))/(c*10^-9.0)
"""
Layer 1 - plane x - PMT1
Layer 1 - plane y - PMT2
Layer 2 - plane x - PMT3
Layer 2 - plane y - PMT4
Layer 3 - plane x - PMT5
Layer 3 - plane y - PMT6
"""
Tp(θ,φ,x₀,y₀,t₀) = [t₀ +T(z[1],θ,φ,x₀) + Tx(1,θ,φ,x₀,y₀),
t₀ + T(z[1],θ,φ,x₀) + Ty(1,θ,φ,x₀,y₀),
t₀ + T(z[2],θ,φ,x₀) + Tx(2,θ,φ,x₀,y₀),
t₀ + T(z[2],θ,φ,x₀) + Ty(1,θ,φ,x₀,y₀),
t₀ + T(z[3],θ,φ,x₀) + Tx(3,θ,φ,x₀,y₀),
t₀ + T(z[3],θ,φ,x₀) + Ty(3,θ,φ,x₀,y₀)]
# Single given μ
# Initial conditions
# Plane #1
t₀₁ = 0
x₀₁ = 80
y₀₁ = 180
θ₁ = 30
φ₁ = -4
Tp(θ₁,φ₁,x₀₁,y₀₁,t₀₁)
#i = linspace(1,18,18)
#plot(i,Tp(θ₁,φ₁,x₀₁,y₀₁,t₀₁))
Out[59]: