Index: /BPM/print_datas.py
===================================================================
--- /BPM/print_datas.py	(revision 797)
+++ /BPM/print_datas.py	(revision 798)
@@ -11,4 +11,33 @@
 basename='bpm_name_bpm_number_x_motor_step_x_motor_mm_y_motor_step_y_motor_mm_Va_Vb_Vc_Vd_Sum_x_libera_mm_y_libera_mm_'
 
+#Big dataset 1 (along X)
+#graph_fr(basedir+basename+"BPM0-E_BPM1-impr_BPM2-C_100pts-horizontal_3pts-vertical_20180801_1.txt", "x_motor_mm", "x_libera_mm",lin_reg = "yes", residu = "yes", rms = "yes",legende = "yes" )#, centrage = "yes", y_minimal = -20., y_maximal = -10.)
+#plt.show()
+#graph_fr(basedir+basename+"BPM0-E_BPM1-impr_BPM2-C_100pts-horizontal_3pts-vertical_20180801_1.txt", "x_motor_mm", "x_libera_mm",lin_reg = "yes", residu = "yes", rms = "yes",legende = "yes" , x_minimal = -50000, x_maximal = 50000) #, centrage = "yes"
+#plt.show()
+#graph_fr(basedir+basename+"BPM0-E_BPM1-impr_BPM2-C_100pts-horizontal_3pts-vertical_20180801_1.txt", "x_motor_mm", "x_libera_mm",lin_reg = "yes", residu = "yes", rms = "yes",legende = "yes" , x_minimal = -50000, x_maximal = 50000) #, centrage = "yes"
+
+
+#Decrochage
+graph_fr(basedir+basename+"BPM0-E_BPM1-impr_BPM2-C_100pts-horizontal_3pts-vertical_20180801_1.txt", "x_motor_mm", "x_libera_mm",lin_reg = "yes", residu = "yes", rms = "yes",legende = "yes", x_minimal = -200000., x_maximal = -150000. )#, centrage = "yes", y_minimal = -20., y_maximal = -10.)
+graph_fr(basedir+basename+"BPM0-E_BPM1-impr_BPM2-C_100pts-horizontal_3pts-vertical_20180801_1.txt", "x_motor_mm", "Va",lin_reg = "yes", residu = "yes", rms = "yes",legende = "yes", x_minimal = -200000., x_maximal = -150000. )#, centrage = "yes", y_minimal = -20., y_maximal = -10.)
+graph_fr(basedir+basename+"BPM0-E_BPM1-impr_BPM2-C_100pts-horizontal_3pts-vertical_20180801_1.txt", "x_motor_mm", "Vb",lin_reg = "yes", residu = "yes", rms = "yes",legende = "yes", x_minimal = -200000., x_maximal = -150000. )#, centrage = "yes", y_minimal = -20., y_maximal = -10.)
+graph_fr(basedir+basename+"BPM0-E_BPM1-impr_BPM2-C_100pts-horizontal_3pts-vertical_20180801_1.txt", "x_motor_mm", "Vc",lin_reg = "yes", residu = "yes", rms = "yes",legende = "yes", x_minimal = -200000., x_maximal = -150000. )#, centrage = "yes", y_minimal = -20., y_maximal = -10.)
+graph_fr(basedir+basename+"BPM0-E_BPM1-impr_BPM2-C_100pts-horizontal_3pts-vertical_20180801_1.txt", "x_motor_mm", "Vd",lin_reg = "yes", residu = "yes", rms = "yes",legende = "yes", x_minimal = -200000., x_maximal = -150000. )#, centrage = "yes", y_minimal = -20., y_maximal = -10.)
+plt.show()
+
+#Big data set 2 (along y)
+#graph_fr(basedir+basename+"BPM0-E_BPM1-impr_BPM2-C_3pts-horizontal_100pts-vertical_20180802_0.txt", "y_motor_step", "y_libera_mm",lin_reg = "yes", residu = "yes", rms = "yes",legende = "yes" , save_figure = 'y_20180802_0')#, centrage = "yes", y_minimal = -20., y_maximal = -10.)
+#plt.show()
+#graph_fr(basedir+basename+"BPM0-E_BPM1-impr_BPM2-C_3pts-horizontal_100pts-vertical_20180802_0.txt", "y_motor_mm", "y_libera_mm",lin_reg = "yes", residu = "yes", rms = "yes",legende = "yes" , centrage = "yes", y_minimal = -20., y_maximal = -10., save_figure = 'y_zoom_20180802_0' )
+#plt.show()
+
+
+#Ultra Big data set 3 (along x)
+#graph_fr(basedir+basename+"BPM0-E_BPM1-impr_BPM2-C_100pts-horizontal_5pts-vertical_20180803_0.txt", "x_motor_step", "x_libera_mm",lin_reg = "yes", residu = "yes", rms = "yes",legende = "yes", save_figure = 'x_20180803_0' )#, centrage = "yes", y_minimal = -20., y_maximal = -10.)
+#plt.show()
+
+
+#tests de decrochage
 #graph_fr(basedir+basename+"test_decrochage_BPM0-E_BPM1-impr_BPM2-C_High_Stat_20180803_0.txt", "y_motor_step", "y_libera_mm",lin_reg = "yes", residu = "yes", rms = "yes",legende = None )#, centrage = "yes", y_minimal = -20., y_maximal = -10.)
 #graph_fr(basedir+basename+"test_decrochage_BPM0-E_BPM1-impr_BPM2-C_High_Stat_20180803_1.txt", "y_motor_step", "y_libera_mm",lin_reg = "yes", residu = "yes", rms = "yes",legende = None )#, centrage = "yes", y_minimal = -20., y_maximal = -10.)
@@ -20,5 +49,5 @@
     print(ifile)
     graph_fr(basedir+basename+"test_decrochage_BPM0-E_BPM1-impr_BPM2-C_20180803_"+str(ifile)+".txt", "y_motor_step", "y_libera_mm",lin_reg = "yes", residu = "yes", rms = "yes",legende = None )#, centrage = "yes", y_minimal = -20., y_maximal = -10.)
-""""
+"""
 
 #plt.show(block = False)
Index: /BPM/print_datas_lib.py
===================================================================
--- /BPM/print_datas_lib.py	(revision 797)
+++ /BPM/print_datas_lib.py	(revision 798)
@@ -8,5 +8,5 @@
 global_index_figure = 0
 
-figure_size_in_hundreds_pixels=(8,6)
+figure_size_in_hundreds_pixels=(12,8)
 
 
@@ -262,11 +262,6 @@
     if y_data_name == "x_motor_step" or y_data_name == "y_motor_step":
         Y = [i[y_data_number]/1000000. for i in data]
-    if legende is None:
+    if not(lin_reg is None):
         plt.plot(X,Y, color+".")
-        plt.errorbar(X,Y,yerr = rms, ecolor = color)
-    else:
-        plt.plot(X,Y, color+".", label = legende)
-        plt.errorbar(X,Y,yerr = rms, ecolor = color)
-    if not(lin_reg is None):
         slope, intercept, r_value, p_value, std_err = stats.linregress(X,Y) # return tuple (pente,ordonnee a l'origine, coef de correlation, p-value, erreur standard de l'estimation)
         if legende is None:
@@ -274,6 +269,14 @@
             plt.plot(X,[slope*i + intercept for i in X], color+"--")
         else:
-            plt.plot(X,[slope*i + intercept for i in X], color+"--", label = legende + " : regression lineaire, \nerreur = "+str(round(std_err,3))+",\ncoefficient de correlation  = "+str(round(r_value,3))+ "\npente de la regression lineaire = " + str(round(slope,3)))            
-            plt.legend(fontsize=15)
+            plt.plot(X,[slope*i + intercept for i in X], color+"--", label = legende + " : reg. lin, \nerr = "+str(round(std_err,3))+",\ncoef correl = "+str(round(r_value,3))+ "\npente = " + str(round(slope,3)))            
+            plt.legend(fontsize=15,bbox_to_anchor=(1.11, 1.11))            
+    else:
+        if legende is None:
+            plt.plot(X,Y, color+".")
+            plt.errorbar(X,Y,yerr = rms, ecolor = color)
+        else:
+            plt.plot(X,Y, color+".", label = legende)
+            plt.errorbar(X,Y,yerr = rms, ecolor = color)
+
     x_label = plot_legend_fr(x_data_name)
     y_label = plot_legend_fr(y_data_name)
@@ -512,5 +515,5 @@
 
 
-def graph_fr(fichier, x_data_name, y_data_name, lin_reg = None, residu = None, rms = None, x_begin_mm = None, x_end_mm = None, y_begin_mm = None, y_end_mm = None, x_minimal = None, x_maximal = None, y_minimal = None,  y_maximal = None, centrage = None, one_bpm = None, legende = "yes"):
+def graph_fr(fichier, x_data_name, y_data_name, lin_reg = None, residu = None, rms = None, x_begin_mm = None, x_end_mm = None, y_begin_mm = None, y_end_mm = None, x_minimal = None, x_maximal = None, y_minimal = None,  y_maximal = None, centrage = None, one_bpm = None, legende = "yes", save_figure = None):
     (list_data_name,data) = read_data(fichier)
     if not(x_data_name in list_data_name):
@@ -526,8 +529,22 @@
     linear_fonction_on_data(list_data_name, data, x_data_name, begin_mm = x_begin_mm, end_mm = x_end_mm)
     linear_fonction_on_data(list_data_name, data, y_data_name, begin_mm = y_begin_mm, end_mm = y_end_mm)
-    data = cut_data(list_data_name, data, x_data_name, inf = x_maximal, equal = x_maximal)
-    data = cut_data(list_data_name, data, x_data_name, sup = x_minimal, equal = x_minimal)
-    data = cut_data(list_data_name, data, y_data_name, inf = y_maximal, equal = y_maximal)
-    data = cut_data(list_data_name, data, y_data_name, sup = y_minimal, equal = y_minimal)
+    print("Data size")
+    print(len(data))
+    if not(x_maximal is None):
+        data = cut_data(list_data_name, data, x_data_name, inf = x_maximal, equal = x_maximal)
+        print("x_maximal "+str(x_maximal))
+        print(len(data))
+    if not(x_minimal is None):
+        data = cut_data(list_data_name, data, x_data_name, sup = x_minimal, equal = x_minimal)
+        print("x_minimal "+str(x_minimal))
+        print(len(data))
+    if not(y_maximal is None):
+        print("y_maximal "+str(y_maximal))
+        data = cut_data(list_data_name, data, y_data_name, inf = y_maximal, equal = y_maximal)
+        print(len(data))
+    if not(y_minimal is None):
+        print("y_minimal "+str(y_minimal))
+        data = cut_data(list_data_name, data, y_data_name, sup = y_minimal, equal = y_minimal)
+        print(len(data))
     if data == []:
         print("cut to high")
@@ -551,24 +568,39 @@
                 if not(legende is None):
                     the_legend=all_bpm_name[bpm_number]
+                    print(the_legend)
                 else:
                     the_legend=None
                 print_graph_fr(list_data_name, cut_data(list_data_name, data, "bpm_name", equal =  all_bpm_name[bpm_number]), x_data_name, y_data_name, global_index_figure, color = color_list[bpm_number], legende = the_legend, lin_reg = lin_reg, rms = rms, centrage = centrage)
+                if not(save_figure is None):
+                    plt.savefig(save_figure+'_'+str(global_index_figure)+'.png')
+                    print("Figure saved")
+                    print(save_figure+'_'+str(global_index_figure)+'.png')
                 if not(residu is None):
                     print_residu_fr(list_data_name,  cut_data(list_data_name, data, "bpm_name", equal = all_bpm_name[bpm_number]), x_data_name, y_data_name, global_index_figure + 1, color_list[bpm_number], legende = the_legend, rms = rms, centrage = centrage)
+                    if not(save_figure is None):
+                        plt.savefig(save_figure+'_residu_'+str(global_index_figure)+'.png')
             global_index_figure += 2
         else:
             all_oposit_y_value = all_value_of_data(list_data_name, data,oposit_y_data_name)
-            if not(legende is None):
-                the_legend=all_bpm_name[bpm_number]
-            else:
-                the_legend=None
             for index_value in range(len(all_oposit_y_value)):
                 data_cut =  cut_data(list_data_name, data, oposit_y_data_name, equal = all_oposit_y_value[index_value])
                 for bpm_number in range(len(all_bpm_name)):
+                    if not(legende is None):
+                        the_legend=all_bpm_name[bpm_number]
+                    else:
+                        the_legend=None
                     #print(cut_data(list_data_name, data_cut, "bpm_name", equal =  all_bpm_name[bpm_number]))
                     print_graph_fr(list_data_name, cut_data(list_data_name, data_cut, "bpm_name", equal =  all_bpm_name[bpm_number]), x_data_name, y_data_name, global_index_figure + 2*index_value, color = color_list[bpm_number], legende = the_legend, lin_reg = lin_reg, rms = rms, title = oposit_y_data_name + " = " + str(all_oposit_y_value[index_value]), centrage = centrage)
                     if not(residu is None):
-                        print_residu_fr(list_data_name, cut_data(list_data_name, data_cut, "bpm_name", equal = all_bpm_name[bpm_number]), x_data_name, y_data_name, global_index_figure + 2*index_value + 1, color_list[bpm_number], legende = the_legend, rms = rms, title = oposit_y_data_name + " = " + str(all_oposit_y_value[index_value]), centrage = centrage)
-            global_index_figure += 2*len(all_oposit_y_value)
+                        print_residu_fr(list_data_name, cut_data(list_data_name, data_cut, "bpm_name", equal = all_bpm_name[bpm_number]), x_data_name, y_data_name, global_index_figure + 2*index_value + 1, color_list[bpm_number], legende = the_legend, rms = rms, title = oposit_y_data_name + " = " + str(all_oposit_y_value[index_value]), centrage = centrage)                        
+                if not(save_figure is None):
+                    plt.figure(global_index_figure + 2*index_value)
+                    plt.savefig(save_figure+'_'+str(global_index_figure + 2*index_value)+'.png')
+                    print("Figure saved")
+                    print(save_figure+'_'+str(global_index_figure + 2*index_value)+'.png')
+                    if not(residu is None):
+                        plt.figure(global_index_figure + 2*index_value+1)
+                        plt.savefig(save_figure+'_residu_'+str(global_index_figure + 2*index_value)+'.png')
+                global_index_figure += 2*len(all_oposit_y_value)
     else: # need to be done, same as befor with "bpm_number" to replace "bpm_name" 
         print("ERROR TYPE OF DATA")
Index: /papers/2016_HDR_ND/Advanced_diags/smithpurcell.tex
===================================================================
--- /papers/2016_HDR_ND/Advanced_diags/smithpurcell.tex	(revision 797)
+++ /papers/2016_HDR_ND/Advanced_diags/smithpurcell.tex	(revision 798)
@@ -145,6 +145,6 @@
 
 \newgeometry{width=200mm,height=280mm,top=5mm,bottom=10mm} 
+\begin{landscape}
 \pagestyle{empty}
-\begin{landscape}
 \hspace{-1cm}
 \begin{table}[htbp]
@@ -212,6 +212,7 @@
 
 \end{landscape}
-\restoregeometry 
-\pagestyle{plain}
+%\restoregeometry 
+%\pagestyle{plain}
+\newgeometry{a4paper,total={145mm,217mm},  top=45mm,headsep=18mm}    
 
 
Index: /papers/2016_HDR_ND/header.tex
===================================================================
--- /papers/2016_HDR_ND/header.tex	(revision 797)
+++ /papers/2016_HDR_ND/header.tex	(revision 798)
@@ -8,5 +8,6 @@
 \documentclass[11pt]{book}
 \usepackage{geometry}        
-\geometry{a4paper}    
+\geometry{a4paper,total={145mm,217mm},  top=65mm,headsep=18mm}    
+%\geometry{a4paper}    
 \usepackage[parfill]{parskip}  
 \usepackage{graphicx}
@@ -23,5 +24,4 @@
             citecolor=blue, urlcolor=blue]{hyperref}
 \usepackage[latin1]{inputenc} 
-\usepackage{geometry}
 %\usepackage[utf8x]{inputenc} 
 
Index: /pepper-pot/pepper-pot.py
===================================================================
--- /pepper-pot/pepper-pot.py	(revision 798)
+++ /pepper-pot/pepper-pot.py	(revision 798)
@@ -0,0 +1,335 @@
+#Beam tracking and pepper-pot simulations
+
+# -*- coding: utf-8 -*-
+#module import
+import numpy as np
+import matplotlib.pyplot as plt
+import time as tm
+
+#gaussian distribution's settings  
+N_particles=10000                                   #no units
+offset_particlesx0p_mm = 0                           #mm (mean/center 1)
+offset_particlesy0p_mm = 0                           #mm (mean/center 2)
+offset_particlesx0d_mrad = 0                         #mrad (mean/center 3)
+offset_particlesy0d_mrad = 0                         #mrad (mean/center 4)
+width_particlesx0p_mm = 1                            #mm (standard deviation 1)
+width_particlesy0p_mm = 1                            #mm (standard deviation 2)
+width_particlesx0d_mrad = 2                          #mrad (standard deviation 3)
+width_particlesy0d_mrad = 2                          #mrad (standard deviation 4)
+slit_width_mm = 0.2                                  #mm 
+#factor = 0.5
+n_slits = 15                                         #no units
+#n_mean_slits = 9
+lower_value_mm = -3                                  #mm 
+upper_value_mm = 3                                   #mm 
+y_slit_positions=np.linspace(lower_value_mm,upper_value_mm,n_slits).tolist() 
+visibility=np.linspace(1,1,N_particles)              #no units
+z_slits_mm = -1                                      #mm
+i_fig = 1                                            #no units
+coeff_mrad_mm=(1./1000)                               #no units  
+#coeff_u=1000                                         #no units
+
+              
+#settings of the initial beam
+def generate_beam():
+    global x0position_mm
+    global y0position_mm
+    global x0divergence_mrad
+    global y0divergence_mrad
+    global visibility
+    x0position_mm=offset_particlesx0p_mm+width_particlesx0p_mm*np.random.randn(N_particles)   #mm
+    y0position_mm=offset_particlesy0p_mm+width_particlesy0p_mm*np.random.randn(N_particles)   #mm   
+    x0divergence_mrad=offset_particlesx0d_mrad+width_particlesx0d_mrad*(np.random.randn(N_particles))          #mrad
+    y0divergence_mrad=offset_particlesy0d_mrad+width_particlesy0d_mrad*(np.random.randn(N_particles))          #mrad
+    visibility=np.linspace(1,1,N_particles)
+
+
+#create slits and return the evolution of visibility      
+def slits(zpos_slits):
+#    MID=[]
+    min_edge_slits=[]
+    max_edge_slits=[]
+    global z_slits_mm                                 #mm
+    global visibility
+    global coeff_mrad_mm
+    global edge    
+#    global middle  
+
+    z_slits_mm=zpos_slits    
+    for i in range(0,len(y_slit_positions)):
+#        print(y_slit_positions[i])                                                   #mm
+#        print(slit_width_mm)                          #mm           
+        visibility=visibility*(1-(abs((y0position_mm+(y0divergence_mrad*z_slits_mm*coeff_mrad_mm))-y_slit_positions[i])<(slit_width_mm/2)*1))
+#        print(visibility)
+#        print(np.sum(visibility))
+#        MID.append(y_slit_positions[i])
+        min_edge_slits.append(y_slit_positions[i]-(slit_width_mm/2))
+        max_edge_slits.append(y_slit_positions[i]+(slit_width_mm/2))
+        edge=min_edge_slits+max_edge_slits
+#        middle=MID
+    return
+
+#plot electron beam (a function of z (the distance between the source and the screen))  #um      
+def plot_beam(z_screen_mm):
+    
+    global i_fig
+    global slit_width_mm 
+    global x1position_mm
+    global y1position_mm
+    global x1positionvisible_mm
+    global y1positionvisible_mm
+    global x1divergencevisible_mrad
+    global y1divergencevisible_mrad
+    global density_beamlet_i
+
+    x1position_mm=x0position_mm+(x0divergence_mrad*z_screen_mm*coeff_mrad_mm)                          
+    y1position_mm=y0position_mm+(y0divergence_mrad*z_screen_mm*coeff_mrad_mm)                          
+    x1positionvisible_mm=[x1position_mm[irange] for irange in range(0,N_particles) if (visibility[irange]==1)]          
+    y1positionvisible_mm=[y1position_mm[irange] for irange in range(0,N_particles) if (visibility[irange]==1)] 
+    x1divergencevisible_mrad=[x0divergence_mrad[irange] for irange in range(0,N_particles) if (visibility[irange]==1)] 
+    y1divergencevisible_mrad=[y0divergence_mrad[irange] for irange in range(0,N_particles) if (visibility[irange]==1)] 
+       
+    emitX= round(emitXYopt_mm_mrad(x1positionvisible_mm,x1divergencevisible_mrad),2)
+    xsq= round(moyennesquare(x1positionvisible_mm),2)
+    xprimesq= round(moyenneprimesquare(x1divergencevisible_mrad),2)
+    correlationx=round(correlation_opt(x1positionvisible_mm,x1divergencevisible_mrad),2)
+    
+    emitY= round(emitXYopt_mm_mrad(y1positionvisible_mm,y1divergencevisible_mrad),2)
+    ysq= round(moyennesquare(y1positionvisible_mm),2)
+    yprimesq= round(moyenneprimesquare(y1divergencevisible_mrad),2)
+    correlationy=round(correlation_opt(y1positionvisible_mm,y1divergencevisible_mrad),2)
+    
+#    X=[2*lower_value_mm,2*upper_value_mm]
+#    print(middle)
+#    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])]
+#    print(y_slit_positions_mean)
+
+    #plot electron beam when particles have visibility=1 after passing through n slits in (x,y) plan
+    plt.figure(i_fig,figsize = (10, 10))
+    #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) )
+    plt.subplot(221)
+    plt.xlabel("X (en mm)")
+    plt.ylabel("Y (en mm)")
+    plt.title("Plan XY")
+    plt.scatter(x1positionvisible_mm,y1positionvisible_mm,s=1,c='red')
+    #plt.tight_layout()
+    #plt.show(block=False)
+    #i_fig=i_fig+1
+
+        
+    #plot electron beam when particles have visibility = 1 after passing through n slits in (x,x') plan
+    plt.subplot(223)
+    plt.title("Plan XX'")
+    plt.xlabel("X (en mm)")
+    plt.ylabel("X' (en mrad)")
+    #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) )
+    plt.scatter(x1positionvisible_mm,x1divergencevisible_mrad,s=1,c='grey')
+#    plt.tight_layout()
+ 
+    #plot electron beam when particles have visibility = 1 after passing through n slits in (y,y') plan
+    #plt.figure(i_fig)
+    #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) )
+    plt.subplot(222)
+    plt.title("Plan YY'")
+    plt.ylabel("Y (en mm)")
+    plt.xlabel("Y' (en mrad)")
+    plt.scatter(y1divergencevisible_mrad,y1positionvisible_mm,s=1,c='blue')
+    #plt.tight_layout()
+        
+    #plt.figure(i_fig)
+    ax = plt.subplot(224)
+    plt.xticks([])
+    plt.yticks([])
+    #ax.set_xticklabels([])
+    #ax.set_yticklabels([])
+    plot_text=" Electron Beam After Passing Through N Slits \n"
+    plot_text+= "width_particlesx0p_mm = " + str(width_particlesx0p_mm) + "\n"
+    plot_text+= "width_particlesx0d_mrad = " + str(width_particlesx0d_mrad) + "\n"
+    plot_text+= "width_particlesy0p_mm = " + str(width_particlesy0p_mm) + "\n"
+    plot_text+= "width_particlesy0d_mrad = " + str(width_particlesy0d_mrad) + "\n"
+    plot_text+= "Zslits = " + str(z_slits_mm) + "mm \n " + " Wslits = "  + str(slit_width_mm) + "mm \n"
+    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) 
+    plt.text(0.05,0.0,plot_text)
+#    plt.show(block=False)
+    i_fig=i_fig+1
+    
+
+    #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) )
+    '''
+    for y_value_edge_slits in edge:
+        plt.plot(X,[y_value_edge_slits, y_value_edge_slits],'k',lw=2)
+    for y_value_middle_slits in middle:
+        plt.plot(X,[y_value_middle_slits, y_value_middle_slits],'b',lw=2)
+    for y_value_middle_slits in y_slit_positions_mean :
+        plt.plot(X,[y_value_middle_slits, y_value_middle_slits],'r',lw=2)
+    plt.scatter(x1positionvisible_mm,y1positionvisible_mm,s=1,c='red')
+    plt.show()
+    i_fig=i_fig+1
+    '''
+    
+    plt.figure(i_fig,figsize = (10, 10))
+    plt.subplot(223)
+    H=plt.hist2d(x1positionvisible_mm,y1positionvisible_mm,100)
+    plt.xlabel("x1positionvisible_mm")
+    plt.ylabel("y1positionvisible_mm")
+    plt.title("HIST 2D XY")
+    plt.set_cmap("plasma")
+    
+    plt.subplot(221)
+    I=plt.hist(x1positionvisible_mm,100)
+    plt.xlabel("x1positionvisible_mm")
+    plt.ylabel("density")
+    plt.title("HIST 1Dx")
+    
+    plt.subplot(222)
+    histo=plt.hist(y1positionvisible_mm,100)
+    plt.xlabel("y1positionvisible_mm")
+    plt.ylabel("density")
+    plt.title("HIST 1Dy")
+#    plt.show()
+    
+    '''
+    density_beamlet_e=np.linspace(0,0,len(middle))
+    for y in y1positionvisible_mm:
+        density_beamlet_e[number_slits(y)]+=1
+    print(density_beamlet_e)
+    
+    
+
+    matrix_histo=np.array(histo)
+    list_histo=matrix_histo.tolist()
+    y_sj=list_histo[1] 
+    print(y_sj)  
+    print(len(y_sj))
+    print(len(density_beamlet_e))
+    y_sj=y_sj.tolist()
+    del y_sj[-1]
+    print(len(y_sj))
+    y_sj_slit=y_slit_positions_mean
+    print(y_sj_slit)
+    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])]
+    print(y_slit_positions_mean_r)
+    y_sj_slit_r=y_slit_positions_mean_r
+    D=[x*y for x, y in zip(density_beamlet_e,y_sj_slit_r)]
+    print(D)
+    E=sum(D)
+    print((1/len(y1positionvisible_mm))*E)
+    
+    
+    #L=[0]
+    #for j in D:
+        #if j!=0:
+            #L[-1]+=j
+        #elif L[-1]!=0:
+            #L+=[0]
+    #x_mean_j=[z/a for z, a in zip(L,density_beamlet_e)]
+    #print(x_mean_j)
+    
+    
+    #plt.figure(i_fig)
+    plt.subplot(224)
+    plt.xticks([])
+    plt.yticks([])
+    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) )
+    i_fig=i_fig+1
+
+    #return density_beamlet_i
+    '''
+
+# /END OF PROCEDURE    
+
+
+#OPTIMIZATION BEAM WITH POSITION
+
+   
+def squarelist(L):    
+    return L**2 
+    
+     
+def absroot(L):
+    return (abs(L)**(1/2))
+
+  
+def moyennesquare(L):
+    #L=x1positionvisible_mm
+    #L=y1positionvisible_mm
+    global x0position_mm
+    global x0divergence_mrad
+    global coeff_mrad_mm
+    global x1position_mm
+    return (1/len(L))*sum(squarelist(L-np.mean(L)))
+
+
+def moyenneprimesquare(M):
+    #L=x1divergencevisible_mrad
+    #L=y1divergencevisible_mrad
+    global N_particles
+    global x0divergence_mrad 
+    return sum((1/len(M))*squarelist(M-np.mean(M)))
+
+   
+def correlation_opt(N,O):
+    #M=x1positionvisible_mm
+    #M=y1positionvisible_mm
+    #N=x1divergencevisible_mrad
+    #N=y1divergencevisible_mrad
+    global x0position_mm
+    global x0divergence_mrad
+    global coeff_mrad_mm
+    global visibility             
+    return (1/len(N))*sum((N-np.mean(N))*(O-np.mean(O)))
+
+
+def emitXYopt_mm_mrad(L,M):
+    return absroot(moyennesquare(L)*moyenneprimesquare(M)-squarelist(correlation_opt(L,M)))
+
+
+def execute(z_slits_mm,z2_screen_mm):
+    generate_beam()
+    slits(z_slits_mm)
+    plot_beam(z2_screen_mm)
+
+
+#TESTS MIN ZHANG
+
+
+def number_slits(L):
+    i=0
+    for y in middle:
+        if L>y:
+            i+=1
+    return i
+
+
+def mean_x_square(L):
+    #L=x1positionvisible_mm
+    #L=y1positionvisible_mm
+    #note julien brossard distance 230
+    #mean_x_beamlet
+    #position moyenne de chacune des fentes
+    
+    global x0position_mm
+    global x0divergence_mrad
+    global coeff_mrad_mm
+    global x1position_mm
+    
+    for i in density_beamlet_i:
+        for j in range(0,len(density_beamlet_i)):
+            return [density_beamlet_i[0],density_beamlet_i[j]]
+    #return (1/len(L))*sum(len(density_beamlet_i)*L)
+    
+    #objective - tuesday : be able to choose j_density_beamlet
+    #execute(200,250)
+    
+
+#def mean_x'_square(L)
+
+# MAIN CODE
+generate_beam()
+plot_beam(100)
+slits(100)
+plot_beam(100)
+plot_beam(500)
+plt.show
+
+
