Commit 6281df41 authored by Yoann Sallaz Damaz's avatar Yoann Sallaz Damaz
Browse files

multrithreading and good beam distribution

parent 526ce3d3
Pipeline #41798 canceled with stages
{
"python.pythonPath": "/usr/bin/python3.8"
}
\ No newline at end of file
from multiprocessing import Pool
from shapely.geometry import Point, LineString, MultiPoint
from shapely.affinity import rotate
from shapely.ops import nearest_points, substring
from math import cos, sin, atan, atan2, degrees, radians, pi
from math import cos, sin, atan, atan2, degrees, radians, pi, exp
import matplotlib.pyplot as plt
import numpy as np
import itertools
import argparse
from scipy import stats
parser = argparse.ArgumentParser()
parser.add_argument('--SL0', metavar="SL0", type=float, help="shift in mm", default=0.0)
parser.add_argument('--M1', metavar="M1", type=float, help="shift in mm", default=0.0)
parser.add_argument('--MOVEH', metavar="MOVEH", type=float, help="shift in mm", default=0.0)
parser.add_argument('--beam_alt', metavar="beam_alt", type=float, help="shift in mm", default=0.0)
parser.add_argument('--beam_ang', metavar="beam_ang", type=float, help="shift in rad", default=0.0)
parser.add_argument('--outfile', metavar="outfile", type=str, help="graph png output name", default="")
args = parser.parse_args()
#################################################################
# DRAW TO DRAW THE BEAM
DRAW = True
DRAW = False
# SCAN TO SCAN AND PLOT DIODE COUNT DURING SCAN
SCAN = True
# Beam ENERGY
NRJ = 12.6 #keV
# NB beam in simu (2 to see beam size, 200 is good for scan)
NB_BEAM = 2
NB_BEAM = 1000
# APPERTURE OF SL0
APERTURE_SL0 = 5 #mm
APERTURE_SL0 = 0.4 #mm
#################################################################
#################################################################
# ERROR ALIGNEMENT
ERROR_ALT_SL0 = 0 #mm
ERROR_ALT_M1 = 0 #mm
ERROR_ALT_MOVEH = 0 #mm
ERROR_ALT_SOURCE = 0 #mm
ERROR_ANG_SOURCE = 0 #rad
ERROR_ALT_SL0 = args.SL0 #mm
ERROR_ALT_M1 = args.M1 #mm
ERROR_ALT_MOVEH = args.MOVEH #mm
ERROR_ALT_SOURCE = args.beam_alt #mm
ERROR_ANG_SOURCE = args.beam_ang #rad
#################################################################
NOT_INTER_DISTANCE = 123456789
def divVertFWHM(E): #nfrom ils
a1=164 #µrad
a2=245 #µrad
b1=0.01375 #kev^-1
b2=0.14418 #kev^-1
return (a1*exp(-b1*E)+a2*exp(-b2*E))*1e-6 #rad FWHM
def divVert(E): #from nils
a1=164 #µrad
a2=245 #µrad
b1=0.01375 #kev^-1
b2=0.14418 #kev^-1
return (a1*exp(-b1*E)+a2*exp(-b2*E))*1e-6*0.849*2 #rad FWHM
def FWHM_to_sigma(FWHM):
return FWHM/2.355
def extendline_from_direction(orig, angle, dist):
Xend = orig.x+dist*cos(angle)
Yend = orig.y+dist*sin(angle)
......@@ -63,10 +96,43 @@ class objet_optique():
ratio_intersection = (intersection.x-self.geo.bounds[0])/(self.geo.bounds[-2]-self.geo.bounds[0])
nearest_geoms = list(substring(self.geo, ratio_intersection-0.002, ratio_intersection+0.002, normalized=True).coords)
return Point(nearest_geoms[0]), Point(nearest_geoms[-1])
def propagation(angle, geometry):
nbHit=0
geometryCopy = geometry.copy()
rayon = [P_Source]
dead = False
while(not dead):
liste_intersection=[]
liste_distance=[]
for obj in geometryCopy:
intersection, distance = obj.intercept(rayon[-1], angle)
liste_intersection.append(intersection)
liste_distance.append(distance)
distance_min = min(liste_distance)
if distance_min == NOT_INTER_DISTANCE:
dead=True
rayon.append(extendline_from_direction(rayon[-1], angle, 1000))
else:
idx = liste_distance.index(distance_min)
obj = geometryCopy[idx]
intersection = liste_intersection[idx]
#print("\trayon hit",obj.name,"at",intersection)
if obj.reflect:
objSurf_P1, objSurf_P2 = obj.get_surface_facette(intersection)
new_angle = angle_reflechi(rayon[-1],intersection,objSurf_P1,objSurf_P2)
angle=new_angle
else:
dead=True
if obj.name=="Diode":
nbHit+=1
rayon.append(intersection)
del(geometryCopy[idx])
return nbHit, rayon
# SOURCE
P_Source = Point(0, ERROR_ALT_SOURCE)
Angle_Source = ERROR_ANG_SOURCE
positionList = []
diodeList = []
......@@ -75,11 +141,14 @@ if DRAW:
scanlist = [0]
SCAN=False
else:
scanlist = list(range(-40,42,2))
for shift in scanlist:
shift/=10.
scanlist = np.linspace(-4, 4, 40)
beam_sigma = FWHM_to_sigma(divVertFWHM(NRJ))
angles_distribution = np.linspace(-2*beam_sigma, 2*beam_sigma, NB_BEAM)
photon_distribution = stats.norm.pdf(angles_distribution, ERROR_ANG_SOURCE, beam_sigma)
for shift in scanlist:
AllObject = []
#############################
......@@ -221,41 +290,15 @@ for shift in scanlist:
col = "k"
plt.plot(x,y,col)
vert_div_source = 0.174E-3 # radian
step = vert_div_source/(NB_BEAM-1)
nbHit = 0
for angle in np.arange(-0.5*vert_div_source+Angle_Source, 0.5*vert_div_source+Angle_Source+step, step):
AllObjectCopy = AllObject.copy()
rayon = [P_Source]
dead = False
while(not dead):
liste_intersection=[]
liste_distance=[]
for obj in AllObjectCopy:
intersection, distance = obj.intercept(rayon[-1], angle)
liste_intersection.append(intersection)
liste_distance.append(distance)
distance_min = min(liste_distance)
if distance_min == NOT_INTER_DISTANCE:
dead=True
rayon.append(extendline_from_direction(rayon[-1], angle, 1000))
else:
idx = liste_distance.index(distance_min)
obj = AllObjectCopy[idx]
intersection = liste_intersection[idx]
#print("\trayon hit",obj.name,"at",intersection)
if obj.reflect:
objSurf_P1, objSurf_P2 = obj.get_surface_facette(intersection)
new_angle = angle_reflechi(rayon[-1],intersection,objSurf_P1,objSurf_P2)
angle=new_angle
else:
dead=True
if obj.name=="Diode":
nbHit+=1
rayon.append(intersection)
del(AllObjectCopy[idx])
if DRAW :
listarg = list(zip(angles_distribution, itertools.repeat(AllObject)))
with Pool() as p:
nbHits, rayons =zip(*p.starmap(propagation, listarg))
nbHit=sum([a*b for a,b in zip(photon_distribution,nbHits)])
if DRAW:
for rayon in rayons:
full_ray = LineString(rayon)
x,y = full_ray.xy
plt.plot(x,y,"g")
......@@ -272,4 +315,7 @@ if SCAN:
plt.plot(positionList,diodeList, marker = 'o')
plt.xlabel("position du scan (mm)")
plt.ylabel("Nb hit diode")
plt.show()
if args.outfile!="":
plt.savefig(args.outfile)
else:
plt.show()
Supports Markdown
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