The python code used for testing:
import numpy as np
import matplotlib.pyplot as plt
eV = 1e-6 #geant const
#surface normal
theFacetNormal = np.array([0,0,1])
#value found during debugging for Au->vacuum
energyThreshold = 5 * eV
#energy of electron
ekint = 6 * eV
N = int(1e5)
phi = np.random.sample(N)*2*np.pi
theta = np.random.sample(N)*np.pi/2
in_arr = np.stack((np.cos(phi)*np.sin(theta), np.sin(phi)*np.sin(theta), -np.cos(theta))).T
# in_arr = np.random.sample((3,N)).T*2 - 1
# %%
def orthogonal(vec):
"""
inline Hep3Vector Hep3Vector::orthogonal() const {
double xx = x() < 0.0 ? -x() : x();
double yy = y() < 0.0 ? -y() : y();
double zz = z() < 0.0 ? -z() : z();
if (xx < yy) {
return xx < zz ? Hep3Vector(0,z(),-y()) : Hep3Vector(y(),-x(),0);
}else{
return yy < zz ? Hep3Vector(-z(),0,x()) : Hep3Vector(y(),-x(),0);
}
}
"""
x = vec[0]
y = vec[1]
z = vec[2]
# return np.array((-1,-1,-1))
# res = np.cross(vec, np.random.sample(3))
if abs(x) < abs(y):
if abs(x) < abs(z):
res = np.array([0, z, -y])
else:
return np.array([y, -x, 0])
else:
if abs(y) < abs(z):
res = np.array([-z, 0, x])
else:
res = np.array([y, -x, 0])
return res/np.linalg.norm(res)
# %%
out_arr = np.zeros((len(in_arr), 3))
for i, incoming in enumerate(in_arr):
# incoming[-1] = -1
#momentum of electron - unit vector
incoming = np.array(incoming)
incoming = incoming / np.linalg.norm(incoming)
#supporting method from G4MicroelecSurface
def GetIncidentAngle(oldMomentum):
PdotN = np.dot(oldMomentum, theFacetNormal)
magP = np.linalg.norm(oldMomentum)
magN = np.linalg.norm(theFacetNormal)
incidentangle = np.asin(PdotN/(magP*magN))
# incidentangle = np.acos(PdotN/(magP*magN))
# incidentangle =
return incidentangle
#a ted uz hura na ty vypocty
thetat = GetIncidentAngle(incoming) #angle d incidence
thetat_alt = np.atan(incoming[2]/np.sqrt(incoming[0]**2 + incoming[1]**2))
# print(f"{incoming}")
# print(f"{thetat = } rad")
# print(f"{thetat_alt = } rad\n")
ekinNormalt = ekint * np.cos(thetat)**2
sin_thetaft = np.sqrt(ekint/(ekint+energyThreshold)) * np.sin(thetat)
#Refraction angle
if (1.0-sin_thetaft**2 ) > 0:
cos_thetaft = np.sqrt(1.0-sin_thetaft**2);
else:
cos_thetaft = 0.0;
# print(cos_thetaft)
#finalni vypocet
cos_thetaft = cos_thetaft*np.cos(thetat) + sin_thetaft*np.sin(thetat)
zVerst = incoming
xVerst = orthogonal(zVerst) #zVerst.orthogonal()
yVerst = np.cross(zVerst,xVerst)
xDirt = np.sqrt(1. - cos_thetaft**2);
yDirt = xDirt;
zPrimeVerst=((xDirt*xVerst + yDirt*yVerst + cos_thetaft*zVerst));
#aParticleChange.ProposeMomentumDirection(zPrimeVerst.unit());
#->unit
unitVecResult = zPrimeVerst / np.linalg.norm(zPrimeVerst)
out_arr[i] = unitVecResult
# print(f" in: {incoming} \nout: {unitVecResult} \n")
out_arr = out_arr.T
in_arr = in_arr.T
# %%
fig, ax = plt.subplots()
for x, d in zip(in_arr, ("x_in", "y_in", "z_in")):
# for x, d in zip(in_arr.T, ("x", "y", "z")):
ax. hist(x,
label = d,
bins = 50,
histtype = "step"
)
ax.legend()
fig, ax = plt.subplots()
for x, d in zip(out_arr, ("x", "y", "z")):
# for x, d in zip(in_arr.T, ("x", "y", "z")):
ax. hist(x,
label = d,
bins = 50,
histtype = "step"
)
ax.legend()
diff_arr = out_arr - in_arr
fig, ax = plt.subplots()
for x, d in zip(diff_arr, ("x_diff", "y_diff", "z_diff")):
# for x, d in zip(in_arr.T, ("x", "y", "z")):
ax. hist(x,
label = d,
bins = 50,
histtype = "step"
)
ax.legend()
fig, ax = plt.subplots()
r = np.sqrt(out_arr[0]**2 + out_arr[1]**2)
out_phi = np.atan2(out_arr[1],out_arr[0])
# out_phi = np.asin(out_arr[1]/r)
# out_phi = np.arccos(out_arr[0]/r)
out_theta = np.arctan(out_arr[2]/r)
for d, name in zip((
out_phi,
out_theta
), (
"phi",
"theta"
)):
ax.hist(d,
bins = 70,
label = name,
histtype = "step"
)
ax.legend()
