Commit 5692e2ec authored by Yoann Sallaz Damaz's avatar Yoann Sallaz Damaz
Browse files

new version

parent d8dcfb7d
Pipeline #42670 failed with stages
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, exp
from BM07RaytracingClass import BM07Raytracing
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)
......@@ -19,303 +14,32 @@ parser.add_argument('--outfile', metavar="outfile", type=str, help="graph png ou
args = parser.parse_args()
#################################################################
# DRAW TO DRAW THE BEAM
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 = 1000
# APPERTURE OF SL0
APERTURE_SL0 = 0.4 #mm
#################################################################
#################################################################
# ERROR ALIGNEMENT
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
#################################################################
bm07 = BM07Raytracing()
bm07.APERTURE_SL0 = 0.4
bm07.SHIFT_ALT_SOURCE = args.beam_alt
bm07.SHIFT_ANG_SOURCE = args.beam_ang
bm07.SHIFT_ALT_SL0 = args.SL0
bm07.SHIFT_ALT_M1 = args.M1
bm07.SHIFT_ALT_MOVEH = args.MOVEH
bm07.NB_BEAM = 10
NOT_INTER_DISTANCE = 123456789
dataX = [-4, -3.5, -3, -2.5, -2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4]
P_Source, photon_distribution, angles_distribution = bm07.creation_source()
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)
return Point(Xend,Yend)
def angle_reflechi(Psource,Pcontact,Psurf1,Psurf2):
angleSource = atan2(Pcontact.y-Psource.y,Pcontact.x-Psource.x)
angleSurface = atan2(Psurf2.y-Psurf1.y,Psurf2.x-Psurf1.x)
#print("angleSource=",degrees(angleSource))
#print("angleSurface=",degrees(angleSurface))
return 2*angleSurface-angleSource
class objet_optique():
def __init__(self, name, liste_coord, centre_rot, angle, shift=0, reflect=False):
self.name = name
self.geo = LineString(liste_coord)
self.centre = centre_rot
self.geo=rotate(self.geo, angle, origin=centre_rot, use_radians=False)
self.reflect = reflect
def intercept(self, begin, angle):
end = extendline_from_direction(begin, angle, 100000)
line = LineString([begin, end])
intersection = line.intersection(self.geo)
if intersection.is_empty:
return None, NOT_INTER_DISTANCE
if not isinstance(intersection, Point):
intersection = list(intersection)[0]
return intersection, LineString([begin, intersection]).length
def get_surface_facette(self, intersection):
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)
positionList = []
diodeList = []
if DRAW:
scanlist = [0]
SCAN=False
else:
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 = []
#############################
# M1
#############################
courbeM1 = 14.31E6 #mm
lengthM1 = 1280 #mm
heightM1 = 90 #mm
P_M1 = Point(29601, shift+ERROR_ALT_M1) #mm
#M1 partie reflechissante
demiAngleArc = atan(0.5*lengthM1/courbeM1)
listePoint=[]
altituderef = cos(demiAngleArc)*courbeM1
for angle in np.arange(-demiAngleArc, demiAngleArc+demiAngleArc*2/100, demiAngleArc*2/100):
listePoint.append(Point(P_M1.x+sin(angle)*courbeM1, P_M1.y+altituderef-cos(angle)*courbeM1))
AllObject.append(objet_optique("M1_refl", listePoint, P_M1, 0.237,reflect=True))
#M1 corp non reflechissant
first=listePoint[0]
end = listePoint[-1]
listePoint=[first, Point(first.x, first.y-heightM1),
Point(end.x,end.y-heightM1), end]
AllObject.append(objet_optique("M1_corps", listePoint, P_M1, 0.237,reflect=False))
#############################
# MASQUE M1
#############################
P_MasqueM1 = Point(28882,shift+ERROR_ALT_M1)
appertureMasqueM1 = 6.8 #mm
listePoint = [Point(P_MasqueM1.x, 0.5*appertureMasqueM1), Point(P_MasqueM1.x, 0.5*appertureMasqueM1+100)]
AllObject.append(objet_optique("Masque_M1_haut", listePoint, P_MasqueM1, 0, reflect=False))
listePoint = [Point(P_MasqueM1.x, -0.5*appertureMasqueM1), Point(P_MasqueM1.x, -0.5*appertureMasqueM1-100)]
AllObject.append(objet_optique("Masque_M1_bas", listePoint, P_MasqueM1, 0, reflect=False))
#############################
# C1
#############################
P_C1 = Point(32336, 22.62+shift+ERROR_ALT_MOVEH) #mm
lengthC1 = 50 #mm
heightC1 = 50 #mm
#C1 partie reflechissante
listePoint = [Point(P_C1.x-0.5*lengthC1, P_C1.y), Point(P_C1.x+0.5*lengthC1, P_C1.y)]
AllObject.append(objet_optique("C1_refl", listePoint, P_C1, -8.514, reflect=True))
#C1 corp non reflechissant
listePoint = [Point(P_C1.x-0.5*lengthC1, P_C1.y), Point(P_C1.x-0.5*lengthC1, P_C1.y+heightC1),
Point(P_C1.x+0.5*lengthC1, P_C1.y+heightC1), Point(P_C1.x+0.5*lengthC1, P_C1.y)]
AllObject.append(objet_optique("C1_corps", listePoint, P_C1, -8.514, reflect=False))
#############################
# C2
#############################
P_C2 = Point(P_C1.x+208.70, P_C1.y-33.01) #mm
lengthC2 = 50 #mm
heightC2 = 50 #mm
#C2 partie reflechissante
listePoint = [Point(P_C2.x-0.5*lengthC2, P_C2.y), Point(P_C2.x+0.5*lengthC2, P_C2.y)]
AllObject.append(objet_optique("C2_refl", listePoint, P_C1, -8.514, reflect=True))
#C2 corp non reflechissant
listePoint = [Point(P_C2.x-0.5*lengthC2, P_C2.y), Point(P_C2.x-0.5*lengthC2, P_C2.y-heightC2),
Point(P_C2.x+0.5*lengthC2, P_C2.y-heightC2), Point(P_C2.x+0.5*lengthC2, P_C2.y)]
AllObject.append(objet_optique("C2_corps", listePoint, P_C1, -8.514, reflect=False))
#############################
# DIODE
#############################
P_Diode = Point(33550,0)
listePoint = [Point(P_Diode.x, P_Diode.y-50), Point(P_Diode.x,P_Diode.y+0)]
AllObject.append(objet_optique("Diode", listePoint, P_Diode, 0, reflect=False))
#############################
# MASQUE FE
#############################
P_MasqueFE = Point(23437,0)
appertureMasqueFE = 8 #mm
listePoint = [Point(P_MasqueFE.x, P_MasqueFE.y+0.5*appertureMasqueFE),
Point(P_MasqueFE.x, P_MasqueFE.y+0.5*appertureMasqueFE+100)]
AllObject.append(objet_optique("Masque_FE_haut", listePoint, P_MasqueFE, 0, reflect=False))
listePoint = [Point(P_MasqueFE.x, P_MasqueFE.y-0.5*appertureMasqueFE),
Point(P_MasqueFE.x, P_MasqueFE.y-0.5*appertureMasqueFE-100)]
AllObject.append(objet_optique("Masque_FE_bas", listePoint, P_MasqueFE, 0, reflect=False))
#############################
# MASQUE BE
#############################
P_MasqueBE = Point(30577,0)
appertureMasqueBE = 20 #mm A VERIF
listePoint = [Point(P_MasqueBE.x, P_MasqueBE.y+appertureMasqueBE),
Point(P_MasqueBE.x, P_MasqueBE.y+appertureMasqueBE+100)]
AllObject.append(objet_optique("Masque_BE_haut", listePoint, P_MasqueBE, 0, reflect=False))
listePoint = [Point(P_MasqueBE.x, P_MasqueBE.y),
Point(P_MasqueBE.x, P_MasqueBE.y-100)]
AllObject.append(objet_optique("Masque_BE_bas", listePoint, P_MasqueBE, 0, reflect=False))
#############################
# SL0
#############################
P_SL0 = Point(23437+1000,shift+ERROR_ALT_SL0) #AVERIF
appertureSL0 = APERTURE_SL0 #mm
listePoint = [Point(P_SL0.x, P_SL0.y+0.5*appertureSL0), Point(P_SL0.x, P_SL0.y+0.5*appertureSL0+100)]
AllObject.append(objet_optique("SL0_haut", listePoint, P_SL0, 0, reflect=False))
listePoint = [Point(P_SL0.x, P_SL0.y-0.5*appertureSL0), Point(P_SL0.x, P_SL0.y-0.5*appertureSL0-100)]
AllObject.append(objet_optique("SL0_bas", listePoint, P_SL0, 0, reflect=False))
#draw All object
if DRAW :
for obj in AllObject:
x,y = obj.geo.xy
if obj.reflect:
col = "b"
else:
col = "k"
plt.plot(x,y,col)
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")
for shift in dataX:
AllObject = bm07.create_geometry()
nbHit, nbHitWeighted, rayons, rayons_w = bm07.propagation(AllObject, P_Source, photon_distribution, angles_distribution)
print("SHIFT, NB HIT IN DIODE = ", shift, nbHit)
positionList.append(shift)
diodeList.append(nbHit)
diodeList.append(nbHitWeighted)
if DRAW :
plt.axis('equal')
plt.plot(dataX,diodeList, marker = 'o')
plt.xlabel("position du scan (mm)")
plt.ylabel("Nb hit diode")
if args.outfile!="":
plt.savefig(args.outfile)
else:
plt.show()
if SCAN:
plt.plot(positionList,diodeList, marker = 'o')
plt.xlabel("position du scan (mm)")
plt.ylabel("Nb hit diode")
if args.outfile!="":
plt.savefig(args.outfile)
else:
plt.show()
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, exp
import numpy as np
from scipy import stats
NOT_INTER_DISTANCE = 123456789 # distance de nom interception
class objet_optique():
def __init__(self, name, liste_coord, centre_rot, angle, shift=0, reflect=False):
self.name = name
self.geo = LineString(liste_coord)
self.centre = centre_rot
self.geo=rotate(self.geo, angle, origin=centre_rot, use_radians=False)
self.reflect = reflect
def extendline_from_direction(self, orig, angle, dist):
Xend = orig.x+dist*cos(angle)
Yend = orig.y+dist*sin(angle)
return Point(Xend,Yend)
def intercept(self, begin, angle):
end = self.extendline_from_direction(begin, angle, 100000)
line = LineString([begin, end])
intersection = line.intersection(self.geo)
if intersection.is_empty:
return None, NOT_INTER_DISTANCE
if not isinstance(intersection, Point):
intersection = list(intersection)[0]
return intersection, LineString([begin, intersection]).length
def get_surface_facette(self, intersection):
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])
class BM07Raytracing:
def __init__(self):
self.SHIFT_ALT_SOURCE=0
self.SHIFT_ANG_SOURCE=0
self.SHIFT_ALT_SL0=0
self.SHIFT_ALT_M1=0
self.SHIFT_ALT_MOVEH=0
self.APERTURE_SL0=6
self.NB_BEAM=1000
self.NRJ = 12.6 #keV
def divVertFWHM(self, 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(self, 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(self, FWHM):
return FWHM/2.355
def extendline_from_direction(self, orig, angle, dist):
Xend = orig.x+dist*cos(angle)
Yend = orig.y+dist*sin(angle)
return Point(Xend,Yend)
def angle_reflechi(self, Psource,Pcontact,Psurf1,Psurf2):
angleSource = atan2(Pcontact.y-Psource.y,Pcontact.x-Psource.x)
angleSurface = atan2(Psurf2.y-Psurf1.y,Psurf2.x-Psurf1.x)
#print("angleSource=",degrees(angleSource))
#print("angleSurface=",degrees(angleSurface))
return 2*angleSurface-angleSource
def create_geometry(self, shiftOptics=0.0):
AllObject = []
#############################
# M1
#############################
courbeM1 = 14.31E6 #mm
lengthM1 = 1280 #mm
heightM1 = 90 #mm
P_M1 = Point(29601, shiftOptics+self.SHIFT_ALT_M1) #mm
#M1 partie reflechissante
demiAngleArc = atan(0.5*lengthM1/courbeM1)
listePoint=[]
altituderef = cos(demiAngleArc)*courbeM1
for angle in np.arange(-demiAngleArc, demiAngleArc+demiAngleArc*2/100, demiAngleArc*2/100):
listePoint.append(Point(P_M1.x+sin(angle)*courbeM1, P_M1.y+altituderef-cos(angle)*courbeM1))
AllObject.append(objet_optique("M1_refl", listePoint, P_M1, 0.237,reflect=True))
#M1 corp non reflechissant
first=listePoint[0]
end = listePoint[-1]
listePoint=[first, Point(first.x, first.y-heightM1),
Point(end.x,end.y-heightM1), end]
AllObject.append(objet_optique("M1_corps", listePoint, P_M1, 0.237,reflect=False))
#############################
# MASQUE M1
#############################
P_MasqueM1 = Point(28882, shiftOptics+self.SHIFT_ALT_M1)
appertureMasqueM1 = 6.8 #mm
listePoint = [Point(P_MasqueM1.x, appertureMasqueM1), Point(P_MasqueM1.x, appertureMasqueM1+100)]
AllObject.append(objet_optique("Masque_M1_haut", listePoint, P_M1, 0.237, reflect=False))
listePoint = [Point(P_MasqueM1.x, 0), Point(P_MasqueM1.x, -100)]
AllObject.append(objet_optique("Masque_M1_bas", listePoint, P_M1, 0.237, reflect=False))
#############################
# TUNGSTEN
#############################
P_TUNGSTEN = Point(31780, 4.9) #mm
listePoint = [Point(P_TUNGSTEN.x, P_TUNGSTEN.y), Point(P_TUNGSTEN.x, P_TUNGSTEN.y-100)]
AllObject.append(objet_optique("Bloc_tungsten", listePoint, P_TUNGSTEN, 0, reflect=False))
#############################
# C1
#############################
P_C1 = Point(32336, 22.62+shiftOptics+self.SHIFT_ALT_MOVEH) #mm
lengthC1 = 50 #mm
heightC1 = 50 #mm
#C1 partie reflechissante
listePoint = [Point(P_C1.x-0.5*lengthC1, P_C1.y), Point(P_C1.x+0.5*lengthC1, P_C1.y)]
AllObject.append(objet_optique("C1_refl", listePoint, P_C1, -8.514, reflect=True))
#C1 corp non reflechissant
listePoint = [Point(P_C1.x-0.5*lengthC1, P_C1.y), Point(P_C1.x-0.5*lengthC1, P_C1.y+heightC1),
Point(P_C1.x+0.5*lengthC1, P_C1.y+heightC1), Point(P_C1.x+0.5*lengthC1, P_C1.y)]
AllObject.append(objet_optique("C1_corps", listePoint, P_C1, -8.514, reflect=False))
#############################
# C2
#############################
P_C2 = Point(P_C1.x+208.70, P_C1.y-33.01) #mm
lengthC2 = 95 #mm
heightC2 = 50 #mm
#C2 partie reflechissante
listePoint = [Point(P_C2.x-0.5*lengthC2, P_C2.y), Point(P_C2.x+0.5*lengthC2, P_C2.y)]
AllObject.append(objet_optique("C2_refl", listePoint, P_C1, -8.514, reflect=True))
#C2 corp non reflechissant
listePoint = [Point(P_C2.x-0.5*lengthC2, P_C2.y), Point(P_C2.x-0.5*lengthC2, P_C2.y-heightC2),
Point(P_C2.x+0.5*lengthC2, P_C2.y-heightC2), Point(P_C2.x+0.5*lengthC2, P_C2.y)]
AllObject.append(objet_optique("C2_corps", listePoint, P_C1, -8.514, reflect=False))
#############################
# MASQUE FE
#############################
P_MasqueFE = Point(23437,0)
appertureMasqueFE = 50 #mm normalement 8mm
listePoint = [Point(P_MasqueFE.x, P_MasqueFE.y+0.5*appertureMasqueFE),
Point(P_MasqueFE.x, P_MasqueFE.y+0.5*appertureMasqueFE+100)]
AllObject.append(objet_optique("Masque_FE_haut", listePoint, P_MasqueFE, 0, reflect=False))
listePoint = [Point(P_MasqueFE.x, P_MasqueFE.y-0.5*appertureMasqueFE),
Point(P_MasqueFE.x, P_MasqueFE.y-0.5*appertureMasqueFE-100)]
AllObject.append(objet_optique("Masque_FE_bas", listePoint, P_MasqueFE, 0, reflect=False))
#############################
# MASQUE BE
#############################
P_MasqueBE = Point(30577,-7.1)
appertureMasqueBE = 30
listePoint = [Point(P_MasqueBE.x, P_MasqueBE.y+appertureMasqueBE),
Point(P_MasqueBE.x, P_MasqueBE.y+appertureMasqueBE+100)]
AllObject.append(objet_optique("Masque_BE_haut", listePoint, P_MasqueBE, 0, reflect=False))
listePoint = [Point(P_MasqueBE.x, P_MasqueBE.y),
Point(P_MasqueBE.x, P_MasqueBE.y-100)]
AllObject.append(objet_optique("Masque_BE_bas", listePoint, P_MasqueBE, 0, reflect=False))
#############################
# SL0
#############################
P_SL0 = Point(26958,shiftOptics+self.SHIFT_ALT_SL0)
appertureSL0 = self.APERTURE_SL0 #mm
listePoint = [Point(P_SL0.x, P_SL0.y+0.5*appertureSL0), Point(P_SL0.x, P_SL0.y+0.5*appertureSL0+100)]
AllObject.append(objet_optique("SL0_haut", listePoint, P_SL0, 0, reflect=False))
listePoint = [Point(P_SL0.x, P_SL0.y-0.5*appertureSL0), Point(P_SL0.x, P_SL0.y-0.5*appertureSL0-100)]
AllObject.append(objet_optique("SL0_bas", listePoint, P_SL0, 0, reflect=False))
return AllObject
def creation_source(self):
# SOURCE
P_Source = Point(0, self.SHIFT_ALT_SOURCE)
beam_sigma = self.FWHM_to_sigma(self.divVertFWHM(self.NRJ))
angles_distribution = np.linspace(-2*beam_sigma, 2*beam_sigma, self.NB_BEAM)
photon_distribution = stats.norm.pdf(angles_distribution, 0, beam_sigma)
angles_distribution+=self.SHIFT_ANG_SOURCE
return P_Source, photon_distribution, angles_distribution
def propagation(self, geometry, P_Source, photon_distribution, angles_distribution):
nbHit=0
nbHitWeighted=0
rayons = []
rayons_w = []
for idxPhoton in range(0, len(angles_distribution)):
angle = angles_distribution[idxPhoton]
poids = photon_distribution[idxPhoton]
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(self.extendline_from_direction(rayon[-1], angle, 33200))
else:
idx = liste_distance.index(distance_min)
obj = geometryCopy[idx]
intersection = liste_intersection[idx]
#print("\trayon hit",obj.name,"at",intersection)
if obj.reflect: