ModuleDetector
Provides a class which represents a virtual detector —a plane in 3D space— with methods to analyze the impact points of a ray bundle intercepting it.
Created in Sept 2020
@author: Stefan Haessler + André Kalouguine
1""" 2Provides a class which represents a virtual detector —a plane in 3D space— with methods 3to analyze the impact points of a ray bundle intercepting it. 4 5 6 7 8 9Created in Sept 2020 10 11@author: Stefan Haessler + André Kalouguine 12""" 13# %% Modules 14 15import numpy as np 16import quaternion 17from abc import ABC, abstractmethod 18 19import ARTcore.ModuleProcessing as mp 20import ARTcore.ModuleOpticalRay as mray 21import ARTcore.ModuleGeometry as mgeo 22from ARTcore.ModuleGeometry import Point, Vector, Origin 23 24# This is to allow multiple constructors based on the type of the first argument 25from functools import singledispatchmethod 26import logging 27from copy import copy 28import time 29 30 31logger = logging.getLogger(__name__) 32 33LightSpeed = 299792458000 34 35# %% Abstract class 36class Detector(ABC): 37 """ 38 Abstract class describing a detector and the properties/methods it has to have. 39 It implements the required components for the detector to behave as a 3D object. 40 Non-abstract subclasses of `Detector` should implement the remaining 41 functiuons such as the detection of rays 42 """ 43 type = "Generic Detector (Abstract)" 44 @property 45 def r(self): 46 """ 47 Return the offset of the optical element from its reference frame to the lab reference frame. 48 Not that this is **NOT** the position of the optical element's center point. Rather it is the 49 offset of the reference frame of the optical element from the lab reference frame. 50 """ 51 return self._r 52 53 @r.setter 54 def r(self, NewPosition): 55 """ 56 Set the offset of the optical element from its reference frame to the lab reference frame. 57 """ 58 if isinstance(NewPosition, mgeo.Point) and len(NewPosition) == 3: 59 self._r = NewPosition 60 else: 61 raise TypeError("Position must be a 3D mgeo.Point.") 62 @property 63 def q(self): 64 """ 65 Return the orientation of the optical element. 66 The orientation is stored as a unit quaternion representing the rotation from the optic's coordinate frame to the lab frame. 67 """ 68 return self._q 69 70 @q.setter 71 def q(self, NewOrientation): 72 """ 73 Set the orientation of the optical element. 74 This function normalizes the input quaternion before storing it. 75 If the input is not a quaternion, raise a TypeError. 76 """ 77 if isinstance(NewOrientation, np.quaternion): 78 self._q = NewOrientation.normalized() 79 else: 80 raise TypeError("Orientation must be a quaternion.") 81 82 @property 83 def position(self): 84 """ 85 Return the position of the basepoint of the optical element. Often it is the same as the optical centre. 86 This position is the one around which all rotations are performed. 87 """ 88 return self.r0 + self.r 89 90 @property 91 def orientation(self): 92 """ 93 Return the orientation of the optical element. 94 The utility of this method is unclear. 95 """ 96 return self.q 97 98 @property 99 def basis(self): 100 return self.r0, self.r, self.q 101 102 def __init__(self): 103 self._r = mgeo.Vector([0.0, 0.0, 0.0]) 104 self._q = np.quaternion(1, 0, 0, 0) 105 self.r0 = mgeo.Point([0.0, 0.0, 0.0]) 106 107 @abstractmethod 108 def centre(self): 109 """ 110 (0,0) point of the detector. 111 """ 112 pass 113 114 @abstractmethod 115 def refpoint(self): 116 """ 117 Reference point from which to measure the distance of the detector. 118 Usually the center of the previous optical element. 119 """ 120 pass 121 122 @property 123 def distance(self): 124 """ 125 Return distance of the Detector from its reference point Detector.refpoint. 126 """ 127 return (self.refpoint - self.centre).norm 128 129 @distance.setter 130 def distance(self, NewDistance: float): 131 """ 132 Shift the Detector to the distance NewDistance from its reference point Detector.refpoint. 133 134 Parameters 135 ---------- 136 NewDistance : number 137 The distance (absolute) to which to shift the Detector. 138 """ 139 vector = (self.centre - self.refpoint).normalized 140 self.centre = self.refpoint + vector * NewDistance 141 142 @abstractmethod 143 def get_2D_points(self, RayList): 144 pass 145 146 @abstractmethod 147 def get_3D_points(self, RayList): 148 pass 149 150 @abstractmethod 151 def __copy__(self): 152 pass 153 154 155# %% Infinite plane detector class 156class InfiniteDetector(Detector): 157 """ 158 Simple infinite plane. 159 Beyond being just a solid object, the 160 161 Attributes 162 ---------- 163 centre : np.ndarray 164 3D Point in the Detector plane. 165 166 refpoint : np.ndarray 167 3D reference point from which to measure the Detector distance. 168 """ 169 type = "Infinite Detector" 170 def __init__(self, index=-1): 171 super().__init__() 172 self._centre = mgeo.Origin 173 self._refpoint = mgeo.Origin 174 self.index = index 175 176 def __copy__(self): 177 """ 178 Returns a new Detector object with the same properties. 179 """ 180 result = InfiniteDetector() 181 result._centre = copy(self._centre) 182 result._refpoint = copy(self._refpoint) 183 result._r = copy(self._r) 184 result._q = copy(self._q) 185 result.r0 = copy(self.r0) 186 result.index = self.index 187 return result 188 189 @property 190 def normal(self): 191 return mgeo.Vector([0,0,1]).from_basis(*self.basis) 192 193 @property 194 def centre(self): 195 return self._centre 196 197 # a setter function 198 @centre.setter 199 def centre(self, Centre): 200 if isinstance(Centre, np.ndarray) and Centre.shape == (3,): 201 self._centre = mgeo.Point(Centre) 202 self.r = self._centre # Simply because r0 is 0,0,0 anyways 203 else: 204 raise TypeError("Detector Centre must be a 3D-vector, given as numpy.ndarray of shape (3,).") 205 206 @property 207 def refpoint(self): 208 return self._refpoint 209 210 # a setter function 211 @refpoint.setter 212 def refpoint(self, RefPoint): 213 if isinstance(RefPoint, np.ndarray) and RefPoint.shape == (3,): 214 self._refpoint = mgeo.Point(RefPoint) 215 else: 216 raise TypeError("Detector RefPoint must a 3D-Point") 217 218 # %% Detector placement methods 219 220 @property 221 def distance(self): 222 """ 223 Return distance of the Detector from its reference point Detector.refpoint. 224 """ 225 return (self.refpoint - self.centre).norm 226 227 @distance.setter 228 def distance(self, NewDistance: float): 229 """ 230 Shift the Detector to the distance NewDistance from its reference point Detector.refpoint. 231 232 Parameters 233 ---------- 234 NewDistance : number 235 The distance (absolute) to which to shift the Detector. 236 """ 237 vector = (self.centre - self.refpoint).normalized() 238 self.centre = self.refpoint + vector * NewDistance 239 240 def set_distance(self,x): 241 self.distance = x 242 243 def autoplace(self, RayList, DistanceDetector: float): 244 """ 245 Automatically place and orient the detector such that it is normal to the central ray 246 of the supplied RayList, at the distance DistanceDetector the origin point of that central ray. 247 248 Parameters 249 ---------- 250 RayList : list[Ray] 251 A list of objects of the ModuleOpticalRay.Ray-class. 252 253 DistanceDetector : float 254 The distance at which to place the Detector. 255 """ 256 #CentralRay = mp.FindCentralRay(RayList) 257 CentralRay = RayList[0] 258 if CentralRay is None: 259 logger.warning(f"Could not find central ray! The list of rays has a length of {len(RayList)}") 260 CentralPoint = mgeo.Origin 261 for k in RayList: 262 CentralPoint = CentralPoint + (k.point - mgeo.Origin) 263 CentralPoint = CentralPoint / len(RayList) 264 else: 265 logger.debug(f"Found central ray, using it to position detector: \n{CentralRay}") 266 CentralPoint = CentralRay.point 267 268 self.centre = CentralPoint + CentralRay.vector * DistanceDetector 269 self.refpoint = CentralPoint 270 self.q = mgeo.QRotationVector2Vector(mgeo.Vector([0,0,1]), -CentralRay.vector) 271 272 # %% Detector optimisation methods 273 def test_callback_distances(self, RayList, distances, callback, 274 provide_points=False, 275 detector_reference_frame=False): 276 LocalRayList = [k.to_basis(*self.basis) for k in RayList] 277 N = len(distances) 278 # Calculate the impact points 279 Points = mgeo.IntersectionRayListZPlane(LocalRayList, distances) 280 values = [] 281 for k in range(N): 282 Points[k] = Points[k]._add_dimension() + mgeo.Point([0,0,distances[k]]) 283 values += [callback(RayList, Points[k], self.basis)] 284 return values 285 286 def optimise_distance(self, RayList, Range, callback, 287 detector_reference_frame=False, 288 provide_points=False, 289 maxiter=5, 290 tol=1e-6, 291 splitting=50, 292 Nrays=1000, 293 callback_iteration=None 294 ): 295 """ 296 Optimise the position of the detector within the provided range, trying to 297 maximise some value calculated by the callback function. 298 299 The callback function receives a list of rays that already hit the detector. 300 If requested, the rays can be provided in the reference frame of the detector. 301 302 The callback function should return a single value that is to be minimised. 303 304 The function splits the range into `splitting` points and uses the `IntersectionRayListZPlane` 305 function to calculate the impact points for all of these. 306 If `provide_points` is set to True, the function will pass on the result to the callback function. 307 Otherwise it will generate rays for the callback function, including calculating the paths. 308 309 It will then find the best value and redo the iteration around that value. 310 Repeat until either the maximum number of iterations is reached or the tolerance is met. 311 312 Parameters 313 ---------- 314 RayList : list[Ray] 315 A list of objects of the ModuleOpticalRay.Ray-class. 316 317 Range : tuple 318 The range of distances to consider for the detector placement. 319 320 callback : function 321 The function to be minimised. It should take a list of rays and return a single value. 322 323 detector_reference_frame : bool 324 If True, the rays are provided in the reference frame of the detector. 325 326 provide_points : bool 327 If True, the callback function will receive the impact points directly. 328 329 maxiter : int 330 The maximum number of iterations to perform. 331 332 tol : float 333 The tolerance to reach before stopping the iteration. 334 335 splitting : int 336 The number of points to split the range into. 337 338 Returns 339 ---------- 340 The optimal distance of the detector. 341 """ 342 Range = [i-self.distance for i in Range] 343 if Nrays<len(RayList): 344 RayList = np.random.choice(RayList, Nrays, replace=False) 345 previous_best = None 346 values = [0]*splitting 347 for i in range(maxiter): 348 if callback_iteration is not None: 349 callback_iteration(i, Range) 350 # Split the range into `splitting` points 351 distances = np.linspace(*Range, splitting) 352 values = self.test_callback_distances(RayList, distances, callback, 353 provide_points=provide_points, 354 detector_reference_frame=detector_reference_frame) 355 best = np.argmin(values) 356 # Update the range 357 if best == 0: 358 Range = (distances[0], distances[2]) 359 elif best == splitting - 1: 360 Range = (distances[-3], distances[-1]) 361 else: 362 Range = (distances[best - 1], distances[best + 1]) 363 # Check if the tolerance is met 364 if i>0: 365 if np.abs(values[best] - previous_best) < tol: 366 break 367 previous_best = values[best] 368 self.distance -= distances[best] 369 370 def _spot_size(self, RayList, Points, basis): 371 """ 372 Returns focal spot size. 373 374 Parameters 375 ---------- 376 Points : mgeo.PointArray 377 The points where the rays hit the detector. 378 Returns 379 ---------- 380 float 381 """ 382 center = mgeo.Point(np.mean(Points[:,:2], axis=0)) 383 distances = (Points[:,:2] - center).norm 384 return np.std(distances) 385 386 def _delay_std(self, RayList, Points, basis): 387 """ 388 Returns the standard deviation of the delays of the rays 389 390 Parameters 391 ---------- 392 RayList : list[Ray] 393 A list of objects of the ModuleOpticalRay.Ray-class. 394 395 Returns 396 ---------- 397 float 398 """ 399 #paths = np.array([np.sum(k.path) for k in RayList]) 400 paths = np.sum(np.array([k.path for k in RayList], dtype=float), axis=1) 401 StartingPoints = mgeo.PointArray([k.point for k in RayList]) 402 XYZ = Points.from_basis(*basis) 403 LastDistances = (XYZ - StartingPoints).norm 404 return float(np.std(paths+LastDistances)) 405 406 def _over_intensity(self, RayList, Points, Basis): 407 """ 408 Calculates spot_size * delay_std for the given RayList. 409 """ 410 spot_size = self._spot_size(RayList, Points, Basis) 411 delay_std = self._delay_std(RayList, Points, Basis) 412 return spot_size**2 * delay_std 413 414 # %% Detector response methods 415 def get_3D_points(self, RayList) -> list[np.ndarray]: 416 """ 417 Returns the list of 3D-points in lab-space where the rays in the 418 list RayList hit the detector plane. 419 420 Parameters 421 ---------- 422 RayList : list[Ray] 423 A list of objects of the ModuleOpticalRay.Ray-class. 424 425 Returns 426 ---------- 427 ListPointDetector3D : list[np.ndarray of shape (3,)] 428 """ 429 return self.get_2D_points(RayList)[0]._add_dimension().from_basis(*self.basis) 430 431 def get_2D_points(self, RayList) -> list[np.ndarray]: 432 """ 433 Returns the list of 2D-points in the detector plane, with the origin at Detector.centre. 434 435 Parameters 436 ---------- 437 RayList : list[Ray] 438 A list of objects of the ModuleOpticalRay.Ray-class of length N 439 440 Returns 441 ---------- 442 XY : np.ndarray of shape (N,2) 443 """ 444 return mgeo.IntersectionRayListZPlane([r.to_basis(*self.basis) for r in RayList]) 445 446 def get_centre_2D_points(self, RayList) -> list[np.ndarray]: 447 """ 448 Returns the center of the 2D array of points. 449 450 Parameters 451 ---------- 452 RayList* : list[Ray] 453 A list of objects of the ModuleOpticalRay.Ray-class. 454 455 Returns 456 ---------- 457 ListPointDetector2DCentre : list[np.ndarray of shape (2,)] 458 """ 459 return np.mean(self.get_2D_points(RayList), axis=0) 460 461 def get_Delays(self, RayList) -> list[float]: 462 """ 463 Returns the list of delays of the rays of RayList as they hit the detector plane, 464 calculated as their total path length divided by the speed of light. 465 The delays are relative to the mean “travel time”, i.e. the mean path length of RayList 466 divided by the speed of light. 467 The list index corresponds to that of the RayList. 468 469 Parameters 470 ---------- 471 RayList : list[Ray] 472 A list of objects of the ModuleOpticalRay.Ray-class. 473 474 Returns 475 ---------- 476 DelayList : list[float] 477 """ 478 localRayList = RayList.to_basis(*self.basis) 479 paths = np.sum(RayList.path, axis=1) 480 StartingPoints = mgeo.PointArray([k.point for k in localRayList]) 481 XYZ = mgeo.IntersectionRayListZPlane(localRayList)[0]._add_dimension() 482 LastDistances = (XYZ - StartingPoints).norm 483 TotalPaths = paths+LastDistances 484 MeanPath = np.mean(TotalPaths) 485 return list((TotalPaths-MeanPath) / LightSpeed * 1e15) # in fs
logger =
<Logger ARTcore.ModuleDetector (WARNING)>
LightSpeed =
299792458000
class
Detector(abc.ABC):
37class Detector(ABC): 38 """ 39 Abstract class describing a detector and the properties/methods it has to have. 40 It implements the required components for the detector to behave as a 3D object. 41 Non-abstract subclasses of `Detector` should implement the remaining 42 functiuons such as the detection of rays 43 """ 44 type = "Generic Detector (Abstract)" 45 @property 46 def r(self): 47 """ 48 Return the offset of the optical element from its reference frame to the lab reference frame. 49 Not that this is **NOT** the position of the optical element's center point. Rather it is the 50 offset of the reference frame of the optical element from the lab reference frame. 51 """ 52 return self._r 53 54 @r.setter 55 def r(self, NewPosition): 56 """ 57 Set the offset of the optical element from its reference frame to the lab reference frame. 58 """ 59 if isinstance(NewPosition, mgeo.Point) and len(NewPosition) == 3: 60 self._r = NewPosition 61 else: 62 raise TypeError("Position must be a 3D mgeo.Point.") 63 @property 64 def q(self): 65 """ 66 Return the orientation of the optical element. 67 The orientation is stored as a unit quaternion representing the rotation from the optic's coordinate frame to the lab frame. 68 """ 69 return self._q 70 71 @q.setter 72 def q(self, NewOrientation): 73 """ 74 Set the orientation of the optical element. 75 This function normalizes the input quaternion before storing it. 76 If the input is not a quaternion, raise a TypeError. 77 """ 78 if isinstance(NewOrientation, np.quaternion): 79 self._q = NewOrientation.normalized() 80 else: 81 raise TypeError("Orientation must be a quaternion.") 82 83 @property 84 def position(self): 85 """ 86 Return the position of the basepoint of the optical element. Often it is the same as the optical centre. 87 This position is the one around which all rotations are performed. 88 """ 89 return self.r0 + self.r 90 91 @property 92 def orientation(self): 93 """ 94 Return the orientation of the optical element. 95 The utility of this method is unclear. 96 """ 97 return self.q 98 99 @property 100 def basis(self): 101 return self.r0, self.r, self.q 102 103 def __init__(self): 104 self._r = mgeo.Vector([0.0, 0.0, 0.0]) 105 self._q = np.quaternion(1, 0, 0, 0) 106 self.r0 = mgeo.Point([0.0, 0.0, 0.0]) 107 108 @abstractmethod 109 def centre(self): 110 """ 111 (0,0) point of the detector. 112 """ 113 pass 114 115 @abstractmethod 116 def refpoint(self): 117 """ 118 Reference point from which to measure the distance of the detector. 119 Usually the center of the previous optical element. 120 """ 121 pass 122 123 @property 124 def distance(self): 125 """ 126 Return distance of the Detector from its reference point Detector.refpoint. 127 """ 128 return (self.refpoint - self.centre).norm 129 130 @distance.setter 131 def distance(self, NewDistance: float): 132 """ 133 Shift the Detector to the distance NewDistance from its reference point Detector.refpoint. 134 135 Parameters 136 ---------- 137 NewDistance : number 138 The distance (absolute) to which to shift the Detector. 139 """ 140 vector = (self.centre - self.refpoint).normalized 141 self.centre = self.refpoint + vector * NewDistance 142 143 @abstractmethod 144 def get_2D_points(self, RayList): 145 pass 146 147 @abstractmethod 148 def get_3D_points(self, RayList): 149 pass 150 151 @abstractmethod 152 def __copy__(self): 153 pass
Abstract class describing a detector and the properties/methods it has to have.
It implements the required components for the detector to behave as a 3D object.
Non-abstract subclasses of Detector should implement the remaining
functiuons such as the detection of rays
r
45 @property 46 def r(self): 47 """ 48 Return the offset of the optical element from its reference frame to the lab reference frame. 49 Not that this is **NOT** the position of the optical element's center point. Rather it is the 50 offset of the reference frame of the optical element from the lab reference frame. 51 """ 52 return self._r
Return the offset of the optical element from its reference frame to the lab reference frame. Not that this is NOT the position of the optical element's center point. Rather it is the offset of the reference frame of the optical element from the lab reference frame.
q
63 @property 64 def q(self): 65 """ 66 Return the orientation of the optical element. 67 The orientation is stored as a unit quaternion representing the rotation from the optic's coordinate frame to the lab frame. 68 """ 69 return self._q
Return the orientation of the optical element. The orientation is stored as a unit quaternion representing the rotation from the optic's coordinate frame to the lab frame.
position
83 @property 84 def position(self): 85 """ 86 Return the position of the basepoint of the optical element. Often it is the same as the optical centre. 87 This position is the one around which all rotations are performed. 88 """ 89 return self.r0 + self.r
Return the position of the basepoint of the optical element. Often it is the same as the optical centre. This position is the one around which all rotations are performed.
orientation
91 @property 92 def orientation(self): 93 """ 94 Return the orientation of the optical element. 95 The utility of this method is unclear. 96 """ 97 return self.q
Return the orientation of the optical element. The utility of this method is unclear.
@abstractmethod
def
refpoint(self):
115 @abstractmethod 116 def refpoint(self): 117 """ 118 Reference point from which to measure the distance of the detector. 119 Usually the center of the previous optical element. 120 """ 121 pass
Reference point from which to measure the distance of the detector. Usually the center of the previous optical element.
distance
123 @property 124 def distance(self): 125 """ 126 Return distance of the Detector from its reference point Detector.refpoint. 127 """ 128 return (self.refpoint - self.centre).norm
Return distance of the Detector from its reference point Detector.refpoint.
157class InfiniteDetector(Detector): 158 """ 159 Simple infinite plane. 160 Beyond being just a solid object, the 161 162 Attributes 163 ---------- 164 centre : np.ndarray 165 3D Point in the Detector plane. 166 167 refpoint : np.ndarray 168 3D reference point from which to measure the Detector distance. 169 """ 170 type = "Infinite Detector" 171 def __init__(self, index=-1): 172 super().__init__() 173 self._centre = mgeo.Origin 174 self._refpoint = mgeo.Origin 175 self.index = index 176 177 def __copy__(self): 178 """ 179 Returns a new Detector object with the same properties. 180 """ 181 result = InfiniteDetector() 182 result._centre = copy(self._centre) 183 result._refpoint = copy(self._refpoint) 184 result._r = copy(self._r) 185 result._q = copy(self._q) 186 result.r0 = copy(self.r0) 187 result.index = self.index 188 return result 189 190 @property 191 def normal(self): 192 return mgeo.Vector([0,0,1]).from_basis(*self.basis) 193 194 @property 195 def centre(self): 196 return self._centre 197 198 # a setter function 199 @centre.setter 200 def centre(self, Centre): 201 if isinstance(Centre, np.ndarray) and Centre.shape == (3,): 202 self._centre = mgeo.Point(Centre) 203 self.r = self._centre # Simply because r0 is 0,0,0 anyways 204 else: 205 raise TypeError("Detector Centre must be a 3D-vector, given as numpy.ndarray of shape (3,).") 206 207 @property 208 def refpoint(self): 209 return self._refpoint 210 211 # a setter function 212 @refpoint.setter 213 def refpoint(self, RefPoint): 214 if isinstance(RefPoint, np.ndarray) and RefPoint.shape == (3,): 215 self._refpoint = mgeo.Point(RefPoint) 216 else: 217 raise TypeError("Detector RefPoint must a 3D-Point") 218 219 # %% Detector placement methods 220 221 @property 222 def distance(self): 223 """ 224 Return distance of the Detector from its reference point Detector.refpoint. 225 """ 226 return (self.refpoint - self.centre).norm 227 228 @distance.setter 229 def distance(self, NewDistance: float): 230 """ 231 Shift the Detector to the distance NewDistance from its reference point Detector.refpoint. 232 233 Parameters 234 ---------- 235 NewDistance : number 236 The distance (absolute) to which to shift the Detector. 237 """ 238 vector = (self.centre - self.refpoint).normalized() 239 self.centre = self.refpoint + vector * NewDistance 240 241 def set_distance(self,x): 242 self.distance = x 243 244 def autoplace(self, RayList, DistanceDetector: float): 245 """ 246 Automatically place and orient the detector such that it is normal to the central ray 247 of the supplied RayList, at the distance DistanceDetector the origin point of that central ray. 248 249 Parameters 250 ---------- 251 RayList : list[Ray] 252 A list of objects of the ModuleOpticalRay.Ray-class. 253 254 DistanceDetector : float 255 The distance at which to place the Detector. 256 """ 257 #CentralRay = mp.FindCentralRay(RayList) 258 CentralRay = RayList[0] 259 if CentralRay is None: 260 logger.warning(f"Could not find central ray! The list of rays has a length of {len(RayList)}") 261 CentralPoint = mgeo.Origin 262 for k in RayList: 263 CentralPoint = CentralPoint + (k.point - mgeo.Origin) 264 CentralPoint = CentralPoint / len(RayList) 265 else: 266 logger.debug(f"Found central ray, using it to position detector: \n{CentralRay}") 267 CentralPoint = CentralRay.point 268 269 self.centre = CentralPoint + CentralRay.vector * DistanceDetector 270 self.refpoint = CentralPoint 271 self.q = mgeo.QRotationVector2Vector(mgeo.Vector([0,0,1]), -CentralRay.vector) 272 273 # %% Detector optimisation methods 274 def test_callback_distances(self, RayList, distances, callback, 275 provide_points=False, 276 detector_reference_frame=False): 277 LocalRayList = [k.to_basis(*self.basis) for k in RayList] 278 N = len(distances) 279 # Calculate the impact points 280 Points = mgeo.IntersectionRayListZPlane(LocalRayList, distances) 281 values = [] 282 for k in range(N): 283 Points[k] = Points[k]._add_dimension() + mgeo.Point([0,0,distances[k]]) 284 values += [callback(RayList, Points[k], self.basis)] 285 return values 286 287 def optimise_distance(self, RayList, Range, callback, 288 detector_reference_frame=False, 289 provide_points=False, 290 maxiter=5, 291 tol=1e-6, 292 splitting=50, 293 Nrays=1000, 294 callback_iteration=None 295 ): 296 """ 297 Optimise the position of the detector within the provided range, trying to 298 maximise some value calculated by the callback function. 299 300 The callback function receives a list of rays that already hit the detector. 301 If requested, the rays can be provided in the reference frame of the detector. 302 303 The callback function should return a single value that is to be minimised. 304 305 The function splits the range into `splitting` points and uses the `IntersectionRayListZPlane` 306 function to calculate the impact points for all of these. 307 If `provide_points` is set to True, the function will pass on the result to the callback function. 308 Otherwise it will generate rays for the callback function, including calculating the paths. 309 310 It will then find the best value and redo the iteration around that value. 311 Repeat until either the maximum number of iterations is reached or the tolerance is met. 312 313 Parameters 314 ---------- 315 RayList : list[Ray] 316 A list of objects of the ModuleOpticalRay.Ray-class. 317 318 Range : tuple 319 The range of distances to consider for the detector placement. 320 321 callback : function 322 The function to be minimised. It should take a list of rays and return a single value. 323 324 detector_reference_frame : bool 325 If True, the rays are provided in the reference frame of the detector. 326 327 provide_points : bool 328 If True, the callback function will receive the impact points directly. 329 330 maxiter : int 331 The maximum number of iterations to perform. 332 333 tol : float 334 The tolerance to reach before stopping the iteration. 335 336 splitting : int 337 The number of points to split the range into. 338 339 Returns 340 ---------- 341 The optimal distance of the detector. 342 """ 343 Range = [i-self.distance for i in Range] 344 if Nrays<len(RayList): 345 RayList = np.random.choice(RayList, Nrays, replace=False) 346 previous_best = None 347 values = [0]*splitting 348 for i in range(maxiter): 349 if callback_iteration is not None: 350 callback_iteration(i, Range) 351 # Split the range into `splitting` points 352 distances = np.linspace(*Range, splitting) 353 values = self.test_callback_distances(RayList, distances, callback, 354 provide_points=provide_points, 355 detector_reference_frame=detector_reference_frame) 356 best = np.argmin(values) 357 # Update the range 358 if best == 0: 359 Range = (distances[0], distances[2]) 360 elif best == splitting - 1: 361 Range = (distances[-3], distances[-1]) 362 else: 363 Range = (distances[best - 1], distances[best + 1]) 364 # Check if the tolerance is met 365 if i>0: 366 if np.abs(values[best] - previous_best) < tol: 367 break 368 previous_best = values[best] 369 self.distance -= distances[best] 370 371 def _spot_size(self, RayList, Points, basis): 372 """ 373 Returns focal spot size. 374 375 Parameters 376 ---------- 377 Points : mgeo.PointArray 378 The points where the rays hit the detector. 379 Returns 380 ---------- 381 float 382 """ 383 center = mgeo.Point(np.mean(Points[:,:2], axis=0)) 384 distances = (Points[:,:2] - center).norm 385 return np.std(distances) 386 387 def _delay_std(self, RayList, Points, basis): 388 """ 389 Returns the standard deviation of the delays of the rays 390 391 Parameters 392 ---------- 393 RayList : list[Ray] 394 A list of objects of the ModuleOpticalRay.Ray-class. 395 396 Returns 397 ---------- 398 float 399 """ 400 #paths = np.array([np.sum(k.path) for k in RayList]) 401 paths = np.sum(np.array([k.path for k in RayList], dtype=float), axis=1) 402 StartingPoints = mgeo.PointArray([k.point for k in RayList]) 403 XYZ = Points.from_basis(*basis) 404 LastDistances = (XYZ - StartingPoints).norm 405 return float(np.std(paths+LastDistances)) 406 407 def _over_intensity(self, RayList, Points, Basis): 408 """ 409 Calculates spot_size * delay_std for the given RayList. 410 """ 411 spot_size = self._spot_size(RayList, Points, Basis) 412 delay_std = self._delay_std(RayList, Points, Basis) 413 return spot_size**2 * delay_std 414 415 # %% Detector response methods 416 def get_3D_points(self, RayList) -> list[np.ndarray]: 417 """ 418 Returns the list of 3D-points in lab-space where the rays in the 419 list RayList hit the detector plane. 420 421 Parameters 422 ---------- 423 RayList : list[Ray] 424 A list of objects of the ModuleOpticalRay.Ray-class. 425 426 Returns 427 ---------- 428 ListPointDetector3D : list[np.ndarray of shape (3,)] 429 """ 430 return self.get_2D_points(RayList)[0]._add_dimension().from_basis(*self.basis) 431 432 def get_2D_points(self, RayList) -> list[np.ndarray]: 433 """ 434 Returns the list of 2D-points in the detector plane, with the origin at Detector.centre. 435 436 Parameters 437 ---------- 438 RayList : list[Ray] 439 A list of objects of the ModuleOpticalRay.Ray-class of length N 440 441 Returns 442 ---------- 443 XY : np.ndarray of shape (N,2) 444 """ 445 return mgeo.IntersectionRayListZPlane([r.to_basis(*self.basis) for r in RayList]) 446 447 def get_centre_2D_points(self, RayList) -> list[np.ndarray]: 448 """ 449 Returns the center of the 2D array of points. 450 451 Parameters 452 ---------- 453 RayList* : list[Ray] 454 A list of objects of the ModuleOpticalRay.Ray-class. 455 456 Returns 457 ---------- 458 ListPointDetector2DCentre : list[np.ndarray of shape (2,)] 459 """ 460 return np.mean(self.get_2D_points(RayList), axis=0) 461 462 def get_Delays(self, RayList) -> list[float]: 463 """ 464 Returns the list of delays of the rays of RayList as they hit the detector plane, 465 calculated as their total path length divided by the speed of light. 466 The delays are relative to the mean “travel time”, i.e. the mean path length of RayList 467 divided by the speed of light. 468 The list index corresponds to that of the RayList. 469 470 Parameters 471 ---------- 472 RayList : list[Ray] 473 A list of objects of the ModuleOpticalRay.Ray-class. 474 475 Returns 476 ---------- 477 DelayList : list[float] 478 """ 479 localRayList = RayList.to_basis(*self.basis) 480 paths = np.sum(RayList.path, axis=1) 481 StartingPoints = mgeo.PointArray([k.point for k in localRayList]) 482 XYZ = mgeo.IntersectionRayListZPlane(localRayList)[0]._add_dimension() 483 LastDistances = (XYZ - StartingPoints).norm 484 TotalPaths = paths+LastDistances 485 MeanPath = np.mean(TotalPaths) 486 return list((TotalPaths-MeanPath) / LightSpeed * 1e15) # in fs
Simple infinite plane. Beyond being just a solid object, the
Attributes
centre : np.ndarray
3D Point in the Detector plane.
refpoint : np.ndarray
3D reference point from which to measure the Detector distance.
refpoint
Reference point from which to measure the distance of the detector. Usually the center of the previous optical element.
distance
221 @property 222 def distance(self): 223 """ 224 Return distance of the Detector from its reference point Detector.refpoint. 225 """ 226 return (self.refpoint - self.centre).norm
Return distance of the Detector from its reference point Detector.refpoint.
def
autoplace(self, RayList, DistanceDetector: float):
244 def autoplace(self, RayList, DistanceDetector: float): 245 """ 246 Automatically place and orient the detector such that it is normal to the central ray 247 of the supplied RayList, at the distance DistanceDetector the origin point of that central ray. 248 249 Parameters 250 ---------- 251 RayList : list[Ray] 252 A list of objects of the ModuleOpticalRay.Ray-class. 253 254 DistanceDetector : float 255 The distance at which to place the Detector. 256 """ 257 #CentralRay = mp.FindCentralRay(RayList) 258 CentralRay = RayList[0] 259 if CentralRay is None: 260 logger.warning(f"Could not find central ray! The list of rays has a length of {len(RayList)}") 261 CentralPoint = mgeo.Origin 262 for k in RayList: 263 CentralPoint = CentralPoint + (k.point - mgeo.Origin) 264 CentralPoint = CentralPoint / len(RayList) 265 else: 266 logger.debug(f"Found central ray, using it to position detector: \n{CentralRay}") 267 CentralPoint = CentralRay.point 268 269 self.centre = CentralPoint + CentralRay.vector * DistanceDetector 270 self.refpoint = CentralPoint 271 self.q = mgeo.QRotationVector2Vector(mgeo.Vector([0,0,1]), -CentralRay.vector)
Automatically place and orient the detector such that it is normal to the central ray of the supplied RayList, at the distance DistanceDetector the origin point of that central ray.
Parameters
RayList : list[Ray]
A list of objects of the ModuleOpticalRay.Ray-class.
DistanceDetector : float
The distance at which to place the Detector.
def
test_callback_distances( self, RayList, distances, callback, provide_points=False, detector_reference_frame=False):
274 def test_callback_distances(self, RayList, distances, callback, 275 provide_points=False, 276 detector_reference_frame=False): 277 LocalRayList = [k.to_basis(*self.basis) for k in RayList] 278 N = len(distances) 279 # Calculate the impact points 280 Points = mgeo.IntersectionRayListZPlane(LocalRayList, distances) 281 values = [] 282 for k in range(N): 283 Points[k] = Points[k]._add_dimension() + mgeo.Point([0,0,distances[k]]) 284 values += [callback(RayList, Points[k], self.basis)] 285 return values
def
optimise_distance( self, RayList, Range, callback, detector_reference_frame=False, provide_points=False, maxiter=5, tol=1e-06, splitting=50, Nrays=1000, callback_iteration=None):
287 def optimise_distance(self, RayList, Range, callback, 288 detector_reference_frame=False, 289 provide_points=False, 290 maxiter=5, 291 tol=1e-6, 292 splitting=50, 293 Nrays=1000, 294 callback_iteration=None 295 ): 296 """ 297 Optimise the position of the detector within the provided range, trying to 298 maximise some value calculated by the callback function. 299 300 The callback function receives a list of rays that already hit the detector. 301 If requested, the rays can be provided in the reference frame of the detector. 302 303 The callback function should return a single value that is to be minimised. 304 305 The function splits the range into `splitting` points and uses the `IntersectionRayListZPlane` 306 function to calculate the impact points for all of these. 307 If `provide_points` is set to True, the function will pass on the result to the callback function. 308 Otherwise it will generate rays for the callback function, including calculating the paths. 309 310 It will then find the best value and redo the iteration around that value. 311 Repeat until either the maximum number of iterations is reached or the tolerance is met. 312 313 Parameters 314 ---------- 315 RayList : list[Ray] 316 A list of objects of the ModuleOpticalRay.Ray-class. 317 318 Range : tuple 319 The range of distances to consider for the detector placement. 320 321 callback : function 322 The function to be minimised. It should take a list of rays and return a single value. 323 324 detector_reference_frame : bool 325 If True, the rays are provided in the reference frame of the detector. 326 327 provide_points : bool 328 If True, the callback function will receive the impact points directly. 329 330 maxiter : int 331 The maximum number of iterations to perform. 332 333 tol : float 334 The tolerance to reach before stopping the iteration. 335 336 splitting : int 337 The number of points to split the range into. 338 339 Returns 340 ---------- 341 The optimal distance of the detector. 342 """ 343 Range = [i-self.distance for i in Range] 344 if Nrays<len(RayList): 345 RayList = np.random.choice(RayList, Nrays, replace=False) 346 previous_best = None 347 values = [0]*splitting 348 for i in range(maxiter): 349 if callback_iteration is not None: 350 callback_iteration(i, Range) 351 # Split the range into `splitting` points 352 distances = np.linspace(*Range, splitting) 353 values = self.test_callback_distances(RayList, distances, callback, 354 provide_points=provide_points, 355 detector_reference_frame=detector_reference_frame) 356 best = np.argmin(values) 357 # Update the range 358 if best == 0: 359 Range = (distances[0], distances[2]) 360 elif best == splitting - 1: 361 Range = (distances[-3], distances[-1]) 362 else: 363 Range = (distances[best - 1], distances[best + 1]) 364 # Check if the tolerance is met 365 if i>0: 366 if np.abs(values[best] - previous_best) < tol: 367 break 368 previous_best = values[best] 369 self.distance -= distances[best]
Optimise the position of the detector within the provided range, trying to maximise some value calculated by the callback function.
The callback function receives a list of rays that already hit the detector. If requested, the rays can be provided in the reference frame of the detector.
The callback function should return a single value that is to be minimised.
The function splits the range into splitting points and uses the IntersectionRayListZPlane
function to calculate the impact points for all of these.
If provide_points is set to True, the function will pass on the result to the callback function.
Otherwise it will generate rays for the callback function, including calculating the paths.
It will then find the best value and redo the iteration around that value. Repeat until either the maximum number of iterations is reached or the tolerance is met.
Parameters
RayList : list[Ray]
A list of objects of the ModuleOpticalRay.Ray-class.
Range : tuple
The range of distances to consider for the detector placement.
callback : function
The function to be minimised. It should take a list of rays and return a single value.
detector_reference_frame : bool
If True, the rays are provided in the reference frame of the detector.
provide_points : bool
If True, the callback function will receive the impact points directly.
maxiter : int
The maximum number of iterations to perform.
tol : float
The tolerance to reach before stopping the iteration.
splitting : int
The number of points to split the range into.
Returns
The optimal distance of the detector.
def
get_3D_points(self, RayList) -> list[numpy.ndarray]:
416 def get_3D_points(self, RayList) -> list[np.ndarray]: 417 """ 418 Returns the list of 3D-points in lab-space where the rays in the 419 list RayList hit the detector plane. 420 421 Parameters 422 ---------- 423 RayList : list[Ray] 424 A list of objects of the ModuleOpticalRay.Ray-class. 425 426 Returns 427 ---------- 428 ListPointDetector3D : list[np.ndarray of shape (3,)] 429 """ 430 return self.get_2D_points(RayList)[0]._add_dimension().from_basis(*self.basis)
Returns the list of 3D-points in lab-space where the rays in the list RayList hit the detector plane.
Parameters
RayList : list[Ray]
A list of objects of the ModuleOpticalRay.Ray-class.
Returns
ListPointDetector3D : list[np.ndarray of shape (3,)]
def
get_2D_points(self, RayList) -> list[numpy.ndarray]:
432 def get_2D_points(self, RayList) -> list[np.ndarray]: 433 """ 434 Returns the list of 2D-points in the detector plane, with the origin at Detector.centre. 435 436 Parameters 437 ---------- 438 RayList : list[Ray] 439 A list of objects of the ModuleOpticalRay.Ray-class of length N 440 441 Returns 442 ---------- 443 XY : np.ndarray of shape (N,2) 444 """ 445 return mgeo.IntersectionRayListZPlane([r.to_basis(*self.basis) for r in RayList])
Returns the list of 2D-points in the detector plane, with the origin at Detector.centre.
Parameters
RayList : list[Ray]
A list of objects of the ModuleOpticalRay.Ray-class of length N
Returns
XY : np.ndarray of shape (N,2)
def
get_centre_2D_points(self, RayList) -> list[numpy.ndarray]:
447 def get_centre_2D_points(self, RayList) -> list[np.ndarray]: 448 """ 449 Returns the center of the 2D array of points. 450 451 Parameters 452 ---------- 453 RayList* : list[Ray] 454 A list of objects of the ModuleOpticalRay.Ray-class. 455 456 Returns 457 ---------- 458 ListPointDetector2DCentre : list[np.ndarray of shape (2,)] 459 """ 460 return np.mean(self.get_2D_points(RayList), axis=0)
Returns the center of the 2D array of points.
Parameters
RayList* : list[Ray]
A list of objects of the ModuleOpticalRay.Ray-class.
Returns
ListPointDetector2DCentre : list[np.ndarray of shape (2,)]
def
get_Delays(self, RayList) -> list[float]:
462 def get_Delays(self, RayList) -> list[float]: 463 """ 464 Returns the list of delays of the rays of RayList as they hit the detector plane, 465 calculated as their total path length divided by the speed of light. 466 The delays are relative to the mean “travel time”, i.e. the mean path length of RayList 467 divided by the speed of light. 468 The list index corresponds to that of the RayList. 469 470 Parameters 471 ---------- 472 RayList : list[Ray] 473 A list of objects of the ModuleOpticalRay.Ray-class. 474 475 Returns 476 ---------- 477 DelayList : list[float] 478 """ 479 localRayList = RayList.to_basis(*self.basis) 480 paths = np.sum(RayList.path, axis=1) 481 StartingPoints = mgeo.PointArray([k.point for k in localRayList]) 482 XYZ = mgeo.IntersectionRayListZPlane(localRayList)[0]._add_dimension() 483 LastDistances = (XYZ - StartingPoints).norm 484 TotalPaths = paths+LastDistances 485 MeanPath = np.mean(TotalPaths) 486 return list((TotalPaths-MeanPath) / LightSpeed * 1e15) # in fs
Returns the list of delays of the rays of RayList as they hit the detector plane, calculated as their total path length divided by the speed of light. The delays are relative to the mean “travel time”, i.e. the mean path length of RayList divided by the speed of light. The list index corresponds to that of the RayList.
Parameters
RayList : list[Ray]
A list of objects of the ModuleOpticalRay.Ray-class.
Returns
DelayList : list[float]