Code indexing in gitaly is broken and leads to code not being visible to the user. We work on the issue with highest priority.

Skip to content
Snippets Groups Projects
Commit 3d7c9c7e authored by snuverink_j's avatar snuverink_j
Browse files

add size and angle from static beam

parent 861ec10c
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 129 additions and 56 deletions
beamtheta.py
+ 129
− 56
View file @ 3d7c9c7e
......@@ -30,16 +30,15 @@ def makeGaussian2(x_center=0, y_center=0, theta=0, sigma_x = 10, sigma_y=10, x_s
framesize=30
theta=0
static_std = 8
rot_std = 6
rot_radius = 9
static_radius = 0
z=np.zeros((2*framesize,2*framesize))
z1=np.zeros((2*framesize,2*framesize))
rot_radius = 13
rot_std = 6
static = np.zeros((2*framesize,2*framesize))
rotation = np.zeros((2*framesize,2*framesize))
# +
pi = np.arccos(-1)
......@@ -47,21 +46,21 @@ pi = np.arccos(-1)
for angle in np.arange(0,720,1):
x0=np.cos(pi*angle/180)*rot_radius
y0=np.sin(pi*angle/180)*rot_radius
z += makeGaussian2(x0, y0, angle, rot_std, rot_std, framesize, framesize)
rotation += makeGaussian2(x0, y0, angle, rot_std, rot_std, framesize, framesize)
# -
for angle in np.arange(0,360,1):
x0=np.cos(pi*angle/180)*static_radius
y0=np.sin(pi*angle/180)*static_radius
z1 += makeGaussian2(x0, y0, angle, static_std, static_std, framesize, framesize)
static += makeGaussian2(x0, y0, angle, static_std, static_std, framesize, framesize)
bins = np.arange(-framesize+0.5,framesize+0.5,1)
values_rot = norm(z[framesize,:])
values_stat = norm(z1[framesize,:])
values_rot = norm(rotation[framesize,:])
values_stat = norm(static [framesize,:])
plt.plot(bins, values_rot, color='green', linestyle='--', linewidth=3, label='Rotating beam', marker="x", markerfacecolor='red')
plt.plot(bins, values_stat, color='red', linestyle='--', linewidth=3, label='Static beam', marker="x", markerfacecolor='black')
plt.plot(bins, values_stat, color='red', linestyle='--', linewidth=3, label='Static beam', marker="x", markerfacecolor='black')
plt.ylim(bottom=0)
plt.legend()
plt.xlabel('${radius\ [mm]}$',size=11)
......@@ -94,32 +93,28 @@ rows=2
columns=2
#relative intensities per sq.mm
int_rot = z[30,30]/np.sum(z)
int_static = z1[30,30]/np.sum(z1)
int_rot = rotation[30,30]/np.sum(rotation)
int_static = static[30,30]/np.sum(static)
print(int_rot)
print(int_static)
print("Reduction of peak is", int_rot/int_static*100,"%" )
#Plot 2D Gaussians for unrotated and rotated beams
fig.add_subplot(rows,columns,1)
plt.imshow(z1)
plt.title('2D Gaussian beam distribution z1(x,y) - static beam')
plt.imshow(static)
plt.title('2D Gaussian beam distribution (x,y) - static beam')
plt.clim(0,100)
plt.colorbar()
plt.colorbar();
fig.add_subplot(rows,columns,4)
plt.imshow(z)
plt.title('2D Gaussian beam distribution z(x,y) - rotating beam')
plt.imshow(rotation)
plt.title('2D Gaussian beam distribution (x,y) - rotating beam')
plt.clim(0,100)
plt.colorbar()
plt.xlabel('${x[mm]}$',size=11)
plt.ylabel('${y[mm]}$',size=11)
plt.show()
plt.colorbar();
# # Particle creation
# ## With Gaussian
mu, sigma = 9, 6 # mean and standard deviation
s = np.random.default_rng().normal(mu, sigma, 10000)
......@@ -130,58 +125,136 @@ s
plt.hist(s)
# ?np.zeros
rot_radius = 9
static_radius = 0
static_std = 8
rot_std = 6
np.sin(4.5e-3)
# ## One by one
# Distribution constants
length_wobble_target = 8.35 # [m] length from wobbler to target, to be checked
rot_radius = 9 # [mm]
wobble_angle = rot_radius / length_wobble_target * 1e-3 # before 4.5e-3 # [rad]
static_std = 6 # [mm] symmetric in both planes
length_triplet_target = length_wobble_target + 1.1 # [m] length from triplet to target
beamsize_x_triplet = 0.317 # [mm] TO BE CONFIRMED
beamsize_y_triplet = 2*beamsize_x_triplet
x_angle = static_std / length_triplet_target / beamsize_x_triplet * 1e-3 # 2e-3 # [rad]
y_angle = static_std / length_triplet_target / beamsize_y_triplet * 1e-3
print("wobbled_angle", wobble_angle)
print("x_angle", x_angle)
print("y_angle", y_angle)
# Creation constants
nrPoints = 100
angleStep = 1
particles = np.zeros((nrPoints * 360//angleStep,6))
# +
wobble_angle = 4.5e-3
pi = np.arccos(-1)
particles = np.zeros((nrPoints * int(360/angleStep),6))
for i,angle in enumerate(np.arange(0,360,angleStep)):
x0=np.cos(pi*angle/180)*rot_radius
y0=np.sin(pi*angle/180)*rot_radius
x = np.random.default_rng().normal(x0, sigma, nrPoints)
y = np.random.default_rng().normal(y0, sigma, nrPoints)
xp = x0 / rot_radius * wobble_angle
yp = y0 / rot_radius * wobble_angle
# wobbled angle
xp = np.cos(np.pi*angle/180) * wobble_angle
yp = np.sin(np.pi*angle/180) * wobble_angle
# central position from wobbling
x0 = [xp * rot_radius / wobble_angle if wobble_angle > 0 else xp]
y0 = [yp * rot_radius / wobble_angle if wobble_angle > 0 else yp]
# add static size
x = np.random.default_rng().normal(x0, static_std, nrPoints)
y = np.random.default_rng().normal(y0, static_std, nrPoints)
# correct for static angle
xp += (x-x0) / static_std * x_angle
yp += (y-y0) / static_std * y_angle
# directional cosines
xp = np.sin(xp)
yp = np.sin(yp)
zp = np.cos(wobble_angle)
zp = np.sqrt(1- xp*xp - yp*yp) #np.cos(wobble_angle)
particles[i*nrPoints:(i+1)*nrPoints,:] = np.array([x,np.zeros(nrPoints),y,np.ones(nrPoints)*xp,np.ones(nrPoints)*zp, np.ones(nrPoints)*yp]).T
# -
particles[0]
np.sqrt(particles[:,3]* particles[:,3] + particles[:,4] * particles[:,4] + particles[:,5] * particles[:,5])
plt.title("sin(xp)")
plt.hist(particles[:,3]);
plt.xlabel("x")
plt.ylabel("y")
plt.hist2d(particles[:,0], particles[:,2], bins = 50);
plt.xlabel("x")
plt.ylabel("xp")
plt.hist2d(particles[:,0], particles[:,3], bins = 50);
plt.xlabel("y")
plt.ylabel("yp")
plt.hist2d(particles[:,2], particles[:,5], bins = 50);
plt.title("rotation")
plt.xlabel("sin(xp)")
plt.ylabel("sin(yp)")
plt.hist2d(particles[:,3], particles[:,5], bins = 50);
plt.hist(particles[:,0]);
plt.xlabel("x [mm]")
plt.hist(particles[:,0], bins = 50);
plt.plot(np.arange(0,2*framesize,1),norm(z[framesize,:]), color='green', linestyle='--', linewidth=3, label='Rotating beam', marker="x", markerfacecolor='red')
plt.plot(np.arange(0,2*framesize,1),norm(z1[framesize,:]), color='red', linestyle='--', linewidth=3, label='Static beam', marker="x", markerfacecolor='black')
plt.ylim(bottom=0)
plt.legend()
plt.xlabel('${radius\ [mm]}$',size=11)
plt.ylabel('$norm.\ beam\ intensity$',size=11)
plt.title('Sliced beam distribution')
plt.xlabel("y [mm]")
plt.hist(particles[:,2], bins = 50);
plt.hist(np.sqrt(particles[:,0]*particles[:,0] + particles[:,2]*particles[:,2]), bins = 50);
particles
np.savetxt("flattened.csv", particles, delimiter=",")
# # MCNP format
# +
# Simple in mm
#np.savetxt("flattened.csv", particles, delimiter=",")
# -
import pandas as pd
def MCNPsave():
'''
Save in MCP format
beam travelling in z but y in MCNP
unit cm
first positions then directional cosines
'''
np.savetxt("flattened.csv", particles, delimiter=" ")
df = pd.read_csv('flattened.csv', sep=" ", header=None)
df.columns = ["x", "y", "z", "xp", "yp", "zp"]
# convert to cm
df['x'] = df['x'].div(10)
df['y'] = df['y'].div(10)
df['z'] = df['z'].div(10)
df2 = df[["xp", "yp","zp"]]
df = df.drop('xp', 1)
df = df.drop('yp', 1)
df = df.drop('zp', 1)
df=df.append(df2)
df = df.round(11)
df.insert(0, 'BlS1', '')
df.insert(1, 'BlS2', '')
df.insert(2, 'BlS3', '')
df.insert(3, 'BlS4', '')
df.insert(4, 'BlS5', '')
df.to_csv('flat.txt', sep=' ', float_format='%.11e',index=False, header=None)
# Add header
lines = ['SDEF PAR=h POS=D1 VEC=FPOS D2 DIR=-1 ERG=590 wgt=1', 'SP1 1 35999r','SI1 L ']
with open("flat.txt", 'r+') as file:
readcontent = file.read()
all_lines=file.readlines()
file.seek(0, 0)
for i in lines:
file.write(str(i) + "\n")
file.write(readcontent)
exp = 36003
with open('flat.txt','r') as b:
lines = b.readlines()
with open('flat.txt','w') as b:
for i,line in enumerate(lines):
# Add header for directions
if i == exp:
b.write('DS2 L'+"\n")
b.write(line)
MCNPsave()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment