source: ETALON/pepper-pot/pepper-pot.py

Last change on this file was 798, checked in by delerue, 6 years ago

Added Codesys directory
File size: 14.1 KB
Line 
1#Beam tracking and pepper-pot simulations
2
3# -*- coding: utf-8 -*-
4#module import
5import numpy as np
6import matplotlib.pyplot as plt
7import time as tm
8
9#gaussian distribution's settings 
10N_particles=10000                                   #no units
11offset_particlesx0p_mm = 0                           #mm (mean/center 1)
12offset_particlesy0p_mm = 0                           #mm (mean/center 2)
13offset_particlesx0d_mrad = 0                         #mrad (mean/center 3)
14offset_particlesy0d_mrad = 0                         #mrad (mean/center 4)
15width_particlesx0p_mm = 1                            #mm (standard deviation 1)
16width_particlesy0p_mm = 1                            #mm (standard deviation 2)
17width_particlesx0d_mrad = 2                          #mrad (standard deviation 3)
18width_particlesy0d_mrad = 2                          #mrad (standard deviation 4)
19slit_width_mm = 0.2                                  #mm
20#factor = 0.5
21n_slits = 15                                         #no units
22#n_mean_slits = 9
23lower_value_mm = -3                                  #mm
24upper_value_mm = 3                                   #mm
25y_slit_positions=np.linspace(lower_value_mm,upper_value_mm,n_slits).tolist() 
26visibility=np.linspace(1,1,N_particles)              #no units
27z_slits_mm = -1                                      #mm
28i_fig = 1                                            #no units
29coeff_mrad_mm=(1./1000)                               #no units 
30#coeff_u=1000                                         #no units
31
32             
33#settings of the initial beam
34def generate_beam():
35    global x0position_mm
36    global y0position_mm
37    global x0divergence_mrad
38    global y0divergence_mrad
39    global visibility
40    x0position_mm=offset_particlesx0p_mm+width_particlesx0p_mm*np.random.randn(N_particles)   #mm
41    y0position_mm=offset_particlesy0p_mm+width_particlesy0p_mm*np.random.randn(N_particles)   #mm   
42    x0divergence_mrad=offset_particlesx0d_mrad+width_particlesx0d_mrad*(np.random.randn(N_particles))          #mrad
43    y0divergence_mrad=offset_particlesy0d_mrad+width_particlesy0d_mrad*(np.random.randn(N_particles))          #mrad
44    visibility=np.linspace(1,1,N_particles)
45
46
47#create slits and return the evolution of visibility     
48def slits(zpos_slits):
49#    MID=[]
50    min_edge_slits=[]
51    max_edge_slits=[]
52    global z_slits_mm                                 #mm
53    global visibility
54    global coeff_mrad_mm
55    global edge   
56#    global middle 
57
58    z_slits_mm=zpos_slits   
59    for i in range(0,len(y_slit_positions)):
60#        print(y_slit_positions[i])                                                   #mm
61#        print(slit_width_mm)                          #mm           
62        visibility=visibility*(1-(abs((y0position_mm+(y0divergence_mrad*z_slits_mm*coeff_mrad_mm))-y_slit_positions[i])<(slit_width_mm/2)*1))
63#        print(visibility)
64#        print(np.sum(visibility))
65#        MID.append(y_slit_positions[i])
66        min_edge_slits.append(y_slit_positions[i]-(slit_width_mm/2))
67        max_edge_slits.append(y_slit_positions[i]+(slit_width_mm/2))
68        edge=min_edge_slits+max_edge_slits
69#        middle=MID
70    return
71
72#plot electron beam (a function of z (the distance between the source and the screen))  #um     
73def plot_beam(z_screen_mm):
74   
75    global i_fig
76    global slit_width_mm
77    global x1position_mm
78    global y1position_mm
79    global x1positionvisible_mm
80    global y1positionvisible_mm
81    global x1divergencevisible_mrad
82    global y1divergencevisible_mrad
83    global density_beamlet_i
84
85    x1position_mm=x0position_mm+(x0divergence_mrad*z_screen_mm*coeff_mrad_mm)                         
86    y1position_mm=y0position_mm+(y0divergence_mrad*z_screen_mm*coeff_mrad_mm)                         
87    x1positionvisible_mm=[x1position_mm[irange] for irange in range(0,N_particles) if (visibility[irange]==1)]         
88    y1positionvisible_mm=[y1position_mm[irange] for irange in range(0,N_particles) if (visibility[irange]==1)] 
89    x1divergencevisible_mrad=[x0divergence_mrad[irange] for irange in range(0,N_particles) if (visibility[irange]==1)] 
90    y1divergencevisible_mrad=[y0divergence_mrad[irange] for irange in range(0,N_particles) if (visibility[irange]==1)] 
91       
92    emitX= round(emitXYopt_mm_mrad(x1positionvisible_mm,x1divergencevisible_mrad),2)
93    xsq= round(moyennesquare(x1positionvisible_mm),2)
94    xprimesq= round(moyenneprimesquare(x1divergencevisible_mrad),2)
95    correlationx=round(correlation_opt(x1positionvisible_mm,x1divergencevisible_mrad),2)
96   
97    emitY= round(emitXYopt_mm_mrad(y1positionvisible_mm,y1divergencevisible_mrad),2)
98    ysq= round(moyennesquare(y1positionvisible_mm),2)
99    yprimesq= round(moyenneprimesquare(y1divergencevisible_mrad),2)
100    correlationy=round(correlation_opt(y1positionvisible_mm,y1divergencevisible_mrad),2)
101   
102#    X=[2*lower_value_mm,2*upper_value_mm]
103#    print(middle)
104#    y_slit_positions_mean=[factor*(middle[0]+middle[1]),factor*(middle[1]+middle[2]),factor*(middle[2]+middle[3]),factor*(middle[3]+middle[4]),factor*(middle[4]+middle[5]),factor*(middle[5]+middle[6]),factor*(middle[6]+middle[7]),factor*(middle[7]+middle[8]),factor*(middle[8]+middle[9])]
105#    print(y_slit_positions_mean)
106
107    #plot electron beam when particles have visibility=1 after passing through n slits in (x,y) plan
108    plt.figure(i_fig,figsize = (10, 10))
109    #plt.title(" Electron Beam After Passing Through N Slits \n  Zslits = " + str(z_slits_mm) + "mm \n " + " Wslits = "  + str(slit_width_mm) + "mm \n Zscreen = " + str(z_screen_mm) + "mm\n" + "EmitX= " + str(emitX) + " moyxsquare = " + str(xsq) + " moyxprimesquare = " + str(xprimesq) + "\n" + "correlationx = " + str(correlationx) + "\n" + "EmitY= " + str(emitY) +  "\n" + " moyysquare = " + str(ysq) + " moyyprimesquare = " + str(yprimesq) + " correlationy = " + str(correlationy) )
110    plt.subplot(221)
111    plt.xlabel("X (en mm)")
112    plt.ylabel("Y (en mm)")
113    plt.title("Plan XY")
114    plt.scatter(x1positionvisible_mm,y1positionvisible_mm,s=1,c='red')
115    #plt.tight_layout()
116    #plt.show(block=False)
117    #i_fig=i_fig+1
118
119       
120    #plot electron beam when particles have visibility = 1 after passing through n slits in (x,x') plan
121    plt.subplot(223)
122    plt.title("Plan XX'")
123    plt.xlabel("X (en mm)")
124    plt.ylabel("X' (en mrad)")
125    #plt.title(" Electron Beam After Passing Through N Slits \n  Zslits = " + str(z_slits_mm) + "mm \n " + " Wslits = "  + str(slit_width_mm) + "mm \n Zscreen = " + str(z_screen_mm) + "mm\n" + "EmitX= " + str(emitX) + " moyxsquare = " + str(xsq) + " moyxprimesquare = " + str(xprimesq) + "\n" + "correlationx = " + str(correlationx) + "\n" + "EmitY= " + str(emitY) +  "\n" + " moyysquare = " + str(ysq) + " moyyprimesquare = " + str(yprimesq) + " correlationy = " + str(correlationy) )
126    plt.scatter(x1positionvisible_mm,x1divergencevisible_mrad,s=1,c='grey')
127#    plt.tight_layout()
128 
129    #plot electron beam when particles have visibility = 1 after passing through n slits in (y,y') plan
130    #plt.figure(i_fig)
131    #plt.title(" Electron Beam After Passing Through N Slits \n  Zslits = " + str(z_slits_mm) + "mm \n " + " Wslits = "  + str(slit_width_mm) + "mm \n Zscreen = " + str(z_screen_mm) + "mm\n" + "EmitX= " + str(emitX) + " moyxsquare = " + str(xsq) + " moyxprimesquare = " + str(xprimesq) + "\n" + "correlationx = " + str(correlationx) + "\n" + "EmitY= " + str(emitY) +  "\n" + " moyysquare = " + str(ysq) + " moyyprimesquare = " + str(yprimesq) + " correlationy = " + str(correlationy) )
132    plt.subplot(222)
133    plt.title("Plan YY'")
134    plt.ylabel("Y (en mm)")
135    plt.xlabel("Y' (en mrad)")
136    plt.scatter(y1divergencevisible_mrad,y1positionvisible_mm,s=1,c='blue')
137    #plt.tight_layout()
138       
139    #plt.figure(i_fig)
140    ax = plt.subplot(224)
141    plt.xticks([])
142    plt.yticks([])
143    #ax.set_xticklabels([])
144    #ax.set_yticklabels([])
145    plot_text=" Electron Beam After Passing Through N Slits \n"
146    plot_text+= "width_particlesx0p_mm = " + str(width_particlesx0p_mm) + "\n"
147    plot_text+= "width_particlesx0d_mrad = " + str(width_particlesx0d_mrad) + "\n"
148    plot_text+= "width_particlesy0p_mm = " + str(width_particlesy0p_mm) + "\n"
149    plot_text+= "width_particlesy0d_mrad = " + str(width_particlesy0d_mrad) + "\n"
150    plot_text+= "Zslits = " + str(z_slits_mm) + "mm \n " + " Wslits = "  + str(slit_width_mm) + "mm \n"
151    plot_text+= "Zscreen = " + str(z_screen_mm) + "mm\n" + "EmitX= " + str(emitX) + " moyxsquare = " + str(xsq) + " \n moyxprimesquare = " + str(xprimesq) + "\n" + "correlationx = " + str(correlationx) + "\n" + "EmitY= " + str(emitY) +  "\n" + " moyysquare = " + str(ysq) + "\n moyyprimesquare = " + str(yprimesq) + " correlationy = " + str(correlationy) 
152    plt.text(0.05,0.0,plot_text)
153#    plt.show(block=False)
154    i_fig=i_fig+1
155   
156
157    #plt.title(" Electron Beam After Passing Through N Slits \n  Zslits = " + str(z_slits_mm) + "mm \n " + " Wslits = "  + str(slit_width_mm) + "mm \n Zscreen = " + str(z_screen_mm) + "mm\n" + "EmitX= " + str(emitX) + " moyxsquare = " + str(xsq) + " moyxprimesquare = " + str(xprimesq) + "\n" + "correlationx = " + str(correlationx) + "\n" + "EmitY= " + str(emitY) +  "\n" + " moyysquare = " + str(ysq) + " moyyprimesquare = " + str(yprimesq) + " correlationy = " + str(correlationy) )
158    '''
159    for y_value_edge_slits in edge:
160        plt.plot(X,[y_value_edge_slits, y_value_edge_slits],'k',lw=2)
161    for y_value_middle_slits in middle:
162        plt.plot(X,[y_value_middle_slits, y_value_middle_slits],'b',lw=2)
163    for y_value_middle_slits in y_slit_positions_mean :
164        plt.plot(X,[y_value_middle_slits, y_value_middle_slits],'r',lw=2)
165    plt.scatter(x1positionvisible_mm,y1positionvisible_mm,s=1,c='red')
166    plt.show()
167    i_fig=i_fig+1
168    '''
169   
170    plt.figure(i_fig,figsize = (10, 10))
171    plt.subplot(223)
172    H=plt.hist2d(x1positionvisible_mm,y1positionvisible_mm,100)
173    plt.xlabel("x1positionvisible_mm")
174    plt.ylabel("y1positionvisible_mm")
175    plt.title("HIST 2D XY")
176    plt.set_cmap("plasma")
177   
178    plt.subplot(221)
179    I=plt.hist(x1positionvisible_mm,100)
180    plt.xlabel("x1positionvisible_mm")
181    plt.ylabel("density")
182    plt.title("HIST 1Dx")
183   
184    plt.subplot(222)
185    histo=plt.hist(y1positionvisible_mm,100)
186    plt.xlabel("y1positionvisible_mm")
187    plt.ylabel("density")
188    plt.title("HIST 1Dy")
189#    plt.show()
190   
191    '''
192    density_beamlet_e=np.linspace(0,0,len(middle))
193    for y in y1positionvisible_mm:
194        density_beamlet_e[number_slits(y)]+=1
195    print(density_beamlet_e)
196   
197   
198
199    matrix_histo=np.array(histo)
200    list_histo=matrix_histo.tolist()
201    y_sj=list_histo[1]
202    print(y_sj) 
203    print(len(y_sj))
204    print(len(density_beamlet_e))
205    y_sj=y_sj.tolist()
206    del y_sj[-1]
207    print(len(y_sj))
208    y_sj_slit=y_slit_positions_mean
209    print(y_sj_slit)
210    y_slit_positions_mean_r=[factor*(middle[1]+middle[2]),factor*(middle[2]+middle[3]),factor*(middle[3]+middle[4]),factor*(middle[4]+middle[5]),factor*(middle[5]+middle[6]),factor*(middle[6]+middle[7]),factor*(middle[7]+middle[8]),factor*(middle[8]+middle[9])]
211    print(y_slit_positions_mean_r)
212    y_sj_slit_r=y_slit_positions_mean_r
213    D=[x*y for x, y in zip(density_beamlet_e,y_sj_slit_r)]
214    print(D)
215    E=sum(D)
216    print((1/len(y1positionvisible_mm))*E)
217   
218   
219    #L=[0]
220    #for j in D:
221        #if j!=0:
222            #L[-1]+=j
223        #elif L[-1]!=0:
224            #L+=[0]
225    #x_mean_j=[z/a for z, a in zip(L,density_beamlet_e)]
226    #print(x_mean_j)
227   
228   
229    #plt.figure(i_fig)
230    plt.subplot(224)
231    plt.xticks([])
232    plt.yticks([])
233    plt.text(0.1,0.2," Electron Beam After Passing Through N Slits \n  Zslits = " + str(z_slits_mm) + "mm \n " + " Wslits = "  + str(slit_width_mm) + "mm \n Zscreen = " + str(z_screen_mm) + "mm\n" + "EmitX= " + str(emitX) + " moyxsquare = " + str(xsq) + " moyxprimesquare = " + str(xprimesq) + "\n" + "correlationx = " + str(correlationx) + "\n" + "EmitY= " + str(emitY) +  "\n" + " moyysquare = " + str(ysq) + " moyyprimesquare = " + str(yprimesq) + " correlationy = " + str(correlationy) )
234    i_fig=i_fig+1
235
236    #return density_beamlet_i
237    '''
238
239# /END OF PROCEDURE   
240
241
242#OPTIMIZATION BEAM WITH POSITION
243
244   
245def squarelist(L):   
246    return L**2 
247   
248     
249def absroot(L):
250    return (abs(L)**(1/2))
251
252 
253def moyennesquare(L):
254    #L=x1positionvisible_mm
255    #L=y1positionvisible_mm
256    global x0position_mm
257    global x0divergence_mrad
258    global coeff_mrad_mm
259    global x1position_mm
260    return (1/len(L))*sum(squarelist(L-np.mean(L)))
261
262
263def moyenneprimesquare(M):
264    #L=x1divergencevisible_mrad
265    #L=y1divergencevisible_mrad
266    global N_particles
267    global x0divergence_mrad
268    return sum((1/len(M))*squarelist(M-np.mean(M)))
269
270   
271def correlation_opt(N,O):
272    #M=x1positionvisible_mm
273    #M=y1positionvisible_mm
274    #N=x1divergencevisible_mrad
275    #N=y1divergencevisible_mrad
276    global x0position_mm
277    global x0divergence_mrad
278    global coeff_mrad_mm
279    global visibility             
280    return (1/len(N))*sum((N-np.mean(N))*(O-np.mean(O)))
281
282
283def emitXYopt_mm_mrad(L,M):
284    return absroot(moyennesquare(L)*moyenneprimesquare(M)-squarelist(correlation_opt(L,M)))
285
286
287def execute(z_slits_mm,z2_screen_mm):
288    generate_beam()
289    slits(z_slits_mm)
290    plot_beam(z2_screen_mm)
291
292
293#TESTS MIN ZHANG
294
295
296def number_slits(L):
297    i=0
298    for y in middle:
299        if L>y:
300            i+=1
301    return i
302
303
304def mean_x_square(L):
305    #L=x1positionvisible_mm
306    #L=y1positionvisible_mm
307    #note julien brossard distance 230
308    #mean_x_beamlet
309    #position moyenne de chacune des fentes
310   
311    global x0position_mm
312    global x0divergence_mrad
313    global coeff_mrad_mm
314    global x1position_mm
315   
316    for i in density_beamlet_i:
317        for j in range(0,len(density_beamlet_i)):
318            return [density_beamlet_i[0],density_beamlet_i[j]]
319    #return (1/len(L))*sum(len(density_beamlet_i)*L)
320   
321    #objective - tuesday : be able to choose j_density_beamlet
322    #execute(200,250)
323   
324
325#def mean_x'_square(L)
326
327# MAIN CODE
328generate_beam()
329plot_beam(100)
330slits(100)
331plot_beam(100)
332plot_beam(500)
333plt.show
334
335
Note: See TracBrowser for help on using the repository browser.