CDR plot; generation of beam distribution for Aleksandar

beamtheta.py

@@ -48,6 +48,7 @@ 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)
+# -
 
 
 for angle in np.arange(0,360,1):
@@ -70,6 +71,21 @@ plt.axvline(x= 23,color='black',linestyle='--')
 plt.axvline(x=-23,color='black',linestyle='--')
 plt.text( 23,0.025,'target',ha='center',va='center',rotation='vertical',backgroundcolor='white',bbox={'facecolor':'white','pad':5})
 plt.text(-23,0.025,'target',ha='center',va='center',rotation='vertical',backgroundcolor='white',bbox={'facecolor':'white','pad':5})
+plt.savefig("15.11beamprofile.png",dpi=400)
+
+values_stat
+
+mean_rot  = np.average(bins, weights=values_rot)
+mean_stat = np.average(bins, weights=values_stat)
+print(mean_rot, mean_stat)
+
+var_rot  = np.average((bins - mean_rot) **2, weights=values_rot)
+var_stat = np.average((bins - mean_stat)**2, weights=values_stat)
+print(var_rot, var_stat)
+
+std_rot  = np.sqrt(var_rot)
+std_stat = np.sqrt(var_stat)
+print(std_rot, std_stat)
 
 
 
@@ -78,22 +94,94 @@ rows=2
 columns=2
 
 #relative intensities per sq.mm
-print(z[30,30]/np.sum(z))
-print(z1[30,30]/np.sum(z1))
+int_rot    = z[30,30]/np.sum(z)
+int_static = z1[30,30]/np.sum(z1)
+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.clim(0,100)
-#plt.colorbar()
-
-#fig.add_subplot(rows,columns,4)
-#plt.imshow(z)
-#plt.title('2D Gaussian beam distribution z(x,y) - rotating beam')
-#plt.clim(0,100)
-#plt.colorbar()
-
-#plt.xlabel('${x[mm]}$',size=11)
-#plt.ylabel('${y[mm]}$',size=11)
-#plt.show()
+fig.add_subplot(rows,columns,1)
+plt.imshow(z1)
+plt.title('2D Gaussian beam distribution z1(x,y) - static beam')
+plt.clim(0,100)
+plt.colorbar()
+
+fig.add_subplot(rows,columns,4)
+plt.imshow(z)
+plt.title('2D Gaussian beam distribution z(x,y) - rotating beam')
+plt.clim(0,100)
+plt.colorbar()
+
+plt.xlabel('${x[mm]}$',size=11)
+plt.ylabel('${y[mm]}$',size=11)
+plt.show()
+
+
+
+
+
+mu, sigma = 9, 6 # mean and standard deviation
+s = np.random.default_rng().normal(mu, sigma, 10000)
+
+abs(mu - np.mean(s))
+
+s
+
+plt.hist(s)
+
+# ?np.zeros
+
+rot_radius = 9
+static_radius = 0
+static_std = 8
+rot_std = 6
+
+np.sin(4.5e-3)
+
+nrPoints = 100
+angleStep = 1
+particles = np.zeros((nrPoints * 360//angleStep,6))
+
+# +
+wobble_angle = 4.5e-3 
+pi = np.arccos(-1)
+
+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
+    xp = np.sin(xp)
+    yp = np.sin(yp)
+    zp = 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.hist(particles[:,3]);
+
+plt.hist2d(particles[:,0], particles[:,2], bins = 50);
+
+plt.hist2d(particles[:,3], particles[:,5], bins = 50);
+
+plt.hist(particles[:,0]);
+
+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')
+
+particles
+
+np.savetxt("flattened.csv", particles, delimiter=",")
+
+