Repository navigation
Expand file tree
/
Copy pathPyRate_BeamParameters.py
More file actions
228 lines (182 loc) · 8.03 KB
/
Copy pathPyRate_BeamParameters.py
File metadata and controls
228 lines (182 loc) · 8.03 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
import numpy as np
import math
def calculate_spot_extent(spot_sigma, resolution):
# Work out half xy extent of spot in mm
# Work out pixel value of half spot extent using resolution
spot_fold_sigma_size = 3
spot_sigma_x = spot_sigma[0]
spot_sigma_y = spot_sigma[1]
spot_extent_y = spot_sigma_y * spot_fold_sigma_size
spot_extent_x = spot_sigma_x * spot_fold_sigma_size
spot_pixels_y = math.ceil(spot_extent_y / resolution)
spot_pixels_x = math.ceil(spot_extent_x / resolution)
spot_pixels_max = max(spot_pixels_x, spot_pixels_y)
return spot_pixels_max
def calculate_transmission(energy, dataset):
######################
# MEASUREMENTS HERE! #
######################
if dataset == 'Pre-2023':
# Pre-2023
t1 = 1.46139e-10
t2 = 3.50567
t3 = 2.38552e-18
t4 = 0.155346
elif dataset == 'January-2023':
# Post-2023
t1 = 3.98E-14
t2 = 4.95E+00
t3 = 1.36E-18
t4 = 1.55E-01
elif dataset == 'January-2024':
# 2024
t1 = 3.613e-12
t2 = 4.086
t3 = 1.235e-13
t4 = 1.094e-01
elif dataset == 'January-2025':
# 2025
t1 = 5.64e-13
t2 = 4.46
t3 = 1.196e-17
t4 = 0.147
transmission = (t1*energy**t2) + (t3*math.exp(t4*energy))
return transmission
def calculate_spot_sigma(energy, dataset):
######################
# MEASUREMENTS HERE! #
######################
if dataset == 'Pre-2023':
# pre 2023
spot_sigma_y = 134.60 * energy ** (-0.5515) # mm
spot_sigma_x = 211.56 * energy ** (-0.6939) # mm
elif dataset == 'January-2023':
# January 2023
spot_sigma_y = 190.80 * energy ** (-0.6362) # mm
spot_sigma_x = 344.61 * energy ** (-0.7999) # mm
elif dataset == 'January-2024':
#January 2024
spot_sigma_y = 188.17674274740347 * energy ** (-0.6312043426412497) # mm
spot_sigma_x = 239.49411653405235 * energy ** (-0.7229534938242351) # mm
elif dataset == 'January-2025':
#January 2025
spot_sigma_y = 202.22 * energy ** (-0.6478) # mm
spot_sigma_x = 253.12 * energy ** (-0.7377) # mm
return spot_sigma_x, spot_sigma_y
def gaussian_value(X, Y, mean, sigma_x, sigma_y):
# https://stackoverflow.com/questions/28342968/how-to-plot-a-2d-gaussian-with-different-sigma
# define the covariance matrix
variance_x = sigma_x * sigma_x
variance_y = sigma_y * sigma_y
covariance = np.array([[variance_x, 0.],
[0., variance_y]])
# inverse of covariance matrix
cov_inv = np.linalg.inv(covariance)
# determinant of covariance matrix
cov_det = np.linalg.det(covariance)
# calculate the coefficient
coefficient = 1.0 / ((2 * np.pi) ** 2 * cov_det) ** 0.5
Z = coefficient * np.e ** (-0.5 * (cov_inv[0, 0] * (X - mean[0]) ** 2 + (cov_inv[0, 1] + cov_inv[1, 0]) * (X - mean[0]) * (Y - mean[1]) + cov_inv[1, 1] * (Y - mean[1]) ** 2))
return Z
def simulate_spot(dose_map, spot_position, spot_pixels_max, spot_sigma, resolution, signal_multiplier):
# 11) Work out the i & j range of dose_map to iterate over by using spot position
# and the pixel extent of the spot
min_i = spot_position[0] - spot_pixels_max # x
max_i = spot_position[0] + spot_pixels_max # x
min_j = spot_position[1] - spot_pixels_max # y
max_j = spot_position[1] + spot_pixels_max # y
# drop values less than 0.1% in kernel
max_gaus_value = gaussian_value(0, 0, (0, 0), spot_sigma[0], spot_sigma[1])
cutoff = 0.001 * max_gaus_value
spot_position_mm = [spot_position[0] * resolution, spot_position[1] * resolution]
if min_i < 0:
min_i = 0
if min_j < 0:
min_j = 0
if max_i > dose_map.shape[1]:
max_i = dose_map.shape[1]
if max_j > dose_map.shape[0]:
max_j = dose_map.shape[0]
for j in range(min_j, max_j):
for i in range(min_i, max_i):
# Generate gaussian values
value = gaussian_value(i * resolution, # X
j * resolution, # Y
spot_position_mm, # [mx, my]
spot_sigma[0], # sx
spot_sigma[1]) # sy
if value < cutoff:
dose_map[j][i] = 0
else:
dose_map[j][i] += (value * signal_multiplier)
return dose_map
def create_spot_kernel(spot_details, resolution, signal_multiplier):
spot_pixels_max = spot_details[1]
spot_sigma = spot_details[0]
map = np.zeros((spot_pixels_max * 2, spot_pixels_max * 2)) # (y, x)
spot_position = (int(map.shape[0] / 2), int(map.shape[1] / 2))
map = simulate_spot(map, spot_position, spot_pixels_max, spot_sigma, resolution, signal_multiplier)
return map
def reference_peak_dose_rate_calculations(dataset):
######################
# MEASUREMENTS HERE! #
######################
# Spot size when these measurements were done: 4.651480583164612, 6.477687077054052
# Jack measured dose rate to be 112.7 (works out to 112.655 peak)
if dataset == 'Pre-2023':
measured_spot_sigma_horizontal = 4.651
measured_spot_sigma_vertical = 6.478
measured_udiamond_dose_rate = 112.7
measured_transmission = 0.06893865376297725
measured_cyclotron_current = 800
else:
# Nick's Measurement 08/24
measured_spot_sigma_horizontal = 4.08
measured_spot_sigma_vertical = 5.41
measured_udiamond_dose_rate = 183.7 # +/- 0.4
measured_transmission = 0.07485
measured_cyclotron_current = 621.452830333339 # Calculated as weighted average of the cyclotron logs
######################
# Calculate the relevant parameters for a 245 MeV beam:
# Microdiamond has a 1.1 mm radius (3.8 mm^2 area). Sum the
# contributions from spots across this area and work out the
# average. Measured micro-diamond dose rate should be equated
# to the average and then the peak dose rate can be calculated
# from this.
average_signal = 0.
pixels = 0
min = 999999999
max = 0
for x in np.linspace(-1.2, 1.2, 24):
for y in np.linspace(-1.2, 1.2, 24):
if np.sqrt(x**2 + y**2) <= 1.1:
value = gaussian_value(x,y,(0,0), measured_spot_sigma_horizontal, measured_spot_sigma_vertical)
if value < min:
min = value
elif value > max:
max = value
average_signal += value
pixels += 1
average_signal /= pixels
dose_rate_per_signal = measured_udiamond_dose_rate / average_signal
peak_dose_rate_measured = gaussian_value(0, 0, (0, 0), measured_spot_sigma_horizontal, measured_spot_sigma_vertical) * dose_rate_per_signal
nozzle_current = measured_cyclotron_current * measured_transmission
return dose_rate_per_signal, nozzle_current, peak_dose_rate_measured
def calculate_spot_peak_dose_rate(map, energy, current, clinical_transmission, max_limit_nozzle_current, dataset):
reference_rate_per_signal, reference_nozzle_current, reference_peak_dose_rate = reference_peak_dose_rate_calculations(dataset)
# Work out the dose rate per signal and then apply this to the
# current spot to work out what the peak dose rate for this
# spot is.
centre = [map.shape[0] / 2, map.shape[1] / 2]
signal_at_centre = map[int(centre[0])][int(centre[1])]
peak_dose_rate_this_spot = signal_at_centre * reference_rate_per_signal
transmission = calculate_transmission(energy, dataset)
if clinical_transmission:
transmission = 0.10
nozzle_current_this_beam = current * transmission
if max_limit_nozzle_current:
nozzle_current_this_beam = 215
transmission = nozzle_current_this_beam / 800
nozzle_current_scaling = nozzle_current_this_beam / reference_nozzle_current
dose_rate_at_centre = peak_dose_rate_this_spot * nozzle_current_scaling
return dose_rate_at_centre, transmission, nozzle_current_this_beam, reference_peak_dose_rate