Update glitches/deviation_intensity/files/big_mono_scan.zip,...

Update glitches/deviation_intensity/files/big_mono_scan.zip, glitches/deviation_intensity/files/white_source_beginning_very_interesting.zip, glitches/deviation_intensity/files/plot/mccode_large_source_beginning.dat, glitches/deviation_intensity/files/plot/remove_spikes.py, glitches/deviation_intensity/files/plot/remove_spikes_simplified_large_source.py, glitches/deviation_intensity/files/plot/super_simple.py, glitches/deviation_intensity/files/plot/digitalized_data.txt, glitches/deviation_intensity/files/plot/theory_itsty_firstxtal.dat, glitches/deviation_intensity/files/plot/theory_itsty_mantid.dat, glitches/deviation_intensity/files/plot/mccode.dat

glitches/deviation_intensity/files/plot/remove_spikes.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+c=np.loadtxt("digitalized_data.txt", encoding='utf-8-sig',delimiter=";")
+
+sizec=len(c)
+
+Xdeux_ssrlc = np.zeros(sizec)
+Ydeux_ssrlc = np.zeros(sizec)
+
+i=0
+for elt in c:    
+    Xdeux_ssrlc[i]=1000.0*elt[0]
+    Ydeux_ssrlc[i]=elt[1]
+    i+=1
+
+########### Get data from glitches_positions.txt
+#1, screen variable from previous code
+screen=1
+
+#a=np.loadtxt("theory_itsty_firstxtal.dat") #first and second crystal glitches give the same thing, good
+#a=np.loadtxt("theory_itsty_mine.dat") #not good results
+a=np.loadtxt("theory_itsty_mantid.dat") #not that good imo because there are no 111 planes
+sizea = len(a)
+
+Xdeux_theorique = np.zeros(sizea)
+Ydeux_theorique = np.zeros(sizea)
+Ydeux_theorique_deux = np.zeros(sizea)
+
+i=0
+
+for elt in a:
+    Xdeux_theorique[i]=elt[0]    
+    #Ydeux_theorique[i]=-elt[screen]*(30.69/71350.83256) #first, second xtal
+    #Ydeux_theorique[i]=-elt[screen]*(30.69/311) #f2 311.23, f1 262775
+    Ydeux_theorique[i]=-elt[screen]*(30.69/7408) #(6/595.69)
+    i+=1
+ 
+i=0
+
+aaa=np.loadtxt("theory_itsty_firstxtal.dat")
+
+for elt in aaa:
+    #Xdeux_theorique[i]=elt[0]
+    Ydeux_theorique_deux[i]=-elt[screen]*(30.69/71350.83)
+    i+=1
+  
+ 
+########### Get data from McXtrace's mccode.dat
+#1 3 ou 5
+screen=1
+
+a=np.loadtxt("mccode.dat")
+sizea=len(a)
+
+Xdeux = np.zeros(sizea)
+Ydeux = np.zeros(sizea)
+Ydeux_nan = np.zeros(sizea)
+
+#whats above isnt a glitch
+#whats under is
+threshold = 0.999
+
+not_in_glitch = 1
+i=0
+for elt in a:    
+    Xdeux[i]=elt[0]
+    Ydeux[i]=elt[screen]
+    Ydeux_nan[i]=elt[screen]
+    i+=1
+      
+max_Ydeux = max(Ydeux)    
+
+#### Get rid of glitches, detect where they are and nan the portions, then replace the nans by interpolation       
+i=0
+for elt in a:
+    if i!=0 and i!=(sizea-1):
+      y_diff_sign = Ydeux[i]-Ydeux[i-1]
+      if y_diff_sign<0:
+        y_diff_sign = -1
+        y_div = Ydeux[i]/Ydeux[i-1]
+      else:
+        y_diff_sign = 1
+        y_div = Ydeux[i-1]/Ydeux[i]
+        
+      y_diff = np.abs(Ydeux[i]-Ydeux[i-1])
+            
+      if not_in_glitch==1:
+        if y_div<threshold:
+          not_in_glitch = 0          
+      elif not_in_glitch==0: 
+        if y_div>threshold:
+          if(y_diff_sign>=1):
+            if((Ydeux[i]/Ydeux[i+1])>threshold): #also the point after it
+              not_in_glitch=1
+          
+      if not_in_glitch==0:
+        Ydeux_nan[i] = np.nan
+        
+    i+=1    
+    
+not_nan = np.logical_not(np.isnan(Ydeux_nan))
+indices = np.arange(len(Ydeux_nan))
+without_glitches = np.interp(indices, indices[not_nan], Ydeux_nan[not_nan])
+
+max_without_glitches = max(without_glitches)
+
+for i,elt in enumerate(without_glitches):
+    without_glitches[i]=(without_glitches[i]/max_without_glitches)*100
+
+for i,elt in enumerate(Ydeux):
+    Ydeux[i]=(Ydeux[i]/max_Ydeux)*100
+
+
+#moving avg
+def moving_average(x, w):
+    return np.convolve(x, np.ones(w), 'valid') / w
+    
+#window length of the moving avg
+#nn=500 for the captures ecrans, compare
+nn=500
+o = moving_average(without_glitches, nn)
+y_padded = np.pad(o, (nn//2, nn-1-nn//2), mode='edge')
+
+print("shape",y_padded.shape,without_glitches.shape)
+
+#Look at the normal plot and the plot without the glitches
+#Check that it is correct before looking at the plot of Ydeux 
+#plt.plot(Xdeux, Ydeux,marker='x', color="r")
+#plt.plot(Xdeux, y_padded,marker='x', color="g")
+
+#take off the moving avg without the glitches from the original plot with glitches
+for i,elt in enumerate(Ydeux):
+    Ydeux[i]=Ydeux[i]-y_padded[i]
+
+#optional, you can try to bring a glitch to the same intensity as one of the ssrl glitches to see if the other glitches fit
+for i,elt in enumerate(Ydeux):
+    Ydeux[i]=Ydeux[i]*(30.69/95.69)
+
+
+plt.title("Intensity deviation. Si 220 phi 90.")
+plt.xlabel("E")
+plt.ylabel("I deviation")
+
+
+#plt.plot(Xdeux_theorique,Ydeux_theorique,marker='o', color="b", label='theory new')
+
+plt.plot(Xdeux_theorique, Ydeux_theorique_deux, marker='x',color="y", label='theory')
+
+plt.plot(Xdeux_ssrlc,Ydeux_ssrlc, marker='o',color="r", label='main db ssrl-db BL11-2')
+
+plt.plot(Xdeux, Ydeux, marker='.',color="g", label='McXtrace')
+
+
+plt.legend(fontsize=10)
+plt.show()

glitches/deviation_intensity/files/plot/remove_spikes_simplified_large_source.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import matplotlib.pyplot as plt
+ 
+
+
+c=np.loadtxt("digitalized_data.txt", encoding='utf-8-sig',delimiter=";")
+
+sizec=len(c)
+
+Xdeux_ssrlc = np.zeros(sizec)
+Ydeux_ssrlc = np.zeros(sizec)
+
+i=0
+for elt in c:    
+    Xdeux_ssrlc[i]=1000.0*elt[0]
+    Ydeux_ssrlc[i]=elt[1]
+    i+=1
+
+########### Get data from glitches_positions.txt
+#1, screen variable from previous code
+screen=1
+
+a=np.loadtxt("theory_itsty_firstxtal.dat")
+
+sizea = len(a)
+
+Xdeux_theorique = np.zeros(sizea)
+Ydeux_theorique = np.zeros(sizea)
+
+
+i=0
+
+for elt in a:
+    Xdeux_theorique[i]=elt[0]    
+    Ydeux_theorique[i]=-elt[screen]*(30.69/71350.83)
+    i+=1
+ 
+ 
+ 
+########### Get data from McXtrace's mccode.dat
+#1 3 ou 5
+screen=1
+
+a=np.loadtxt("mccode_large_source_beginning.dat")
+sizea=len(a)
+
+Xdeux = np.zeros(sizea)
+Ydeux = np.zeros(sizea)
+Ydeux_nan = np.zeros(sizea)
+
+#whats above isnt a glitch
+#whats under is
+threshold = 0.999
+
+not_in_glitch = 1
+i=0
+for elt in a:    
+    Xdeux[i]=elt[0]
+    Ydeux[i]=elt[screen]
+    Ydeux_nan[i]=elt[screen]
+    i+=1
+      
+max_Ydeux = max(Ydeux)    
+
+for i,elt in enumerate(Ydeux):
+    Ydeux[i]=(Ydeux[i]/max_Ydeux)*100
+
+
+#moving avg
+def moving_average(x, w):
+    return np.convolve(x, np.ones(w), 'valid') / w
+    
+#window length of the moving avg
+nn= 100
+o = moving_average(Ydeux, nn)
+y_padded = np.pad(o, (nn//2, nn-1-nn//2), mode='edge')
+
+
+
+#Look at the normal plot and the plot without the glitches
+#Check that it is correct before looking at the plot of Ydeux 
+#plt.plot(Xdeux, Ydeux,marker='x', color="r")
+#plt.plot(Xdeux, y_padded,marker='x', color="g")
+
+#take off the moving avg without the glitches from the original plot with glitches
+for i,elt in enumerate(Ydeux):
+    Ydeux[i]=Ydeux[i]-y_padded[i]
+
+plt.title("Intensity deviation. Si 220 phi 90.")
+plt.xlabel("E")
+plt.ylabel("I deviation")
+
+
+plt.plot(Xdeux_theorique,Ydeux_theorique,marker='x', color="y", label='theory')
+
+plt.plot(Xdeux_ssrlc,Ydeux_ssrlc, marker='.',color="r", label='main db ssrl-db BL11-2')
+
+plt.plot(Xdeux, Ydeux, marker='.',color="g", label='McXtrace')
+
+plt.legend(fontsize=10)
+plt.show()
+

glitches/deviation_intensity/files/plot/super_simple.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+import numpy as np
+import matplotlib.pyplot as plt
+import math 
+       
+########### Get data from McXtrace's mccode.dat
+#change value here
+screen=3
+
+a=np.loadtxt("mccode.dat")
+
+sizea=len(a)
+
+Xdeux = np.zeros(sizea)
+Ydeux = np.zeros(sizea)
+
+
+i=0
+for elt in a:    
+    Xdeux[i]=elt[0]
+    Ydeux[i]=elt[screen]
+    i+=1
+
+
+plt.title("Intensity")
+plt.xlabel("E")
+plt.ylabel("Intensity")
+plt.plot(Xdeux,Ydeux,marker='o', color="g")
+plt.show()
+