-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathPMRMC.m
193 lines (167 loc) · 6.21 KB
/
PMRMC.m
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
function record = PMRMC(numOfPositron, magneticField, ...
annihilateThreshold, scatterThreshold, vInitial, isTrack)
% PMRMC (PET/MR Monte Carlo)
% Xiangxi Meng (Beijing Cancer Hospital) and Lujia Jin (Peking Univeristy)
% A Monte Carlo simulation of positron propagation in the magnetic field.
%
% This code presents a minimum environment for Monte Carlo simulation of
% positron propagation, emphasizing on the effect of the magnetic field.
% It is meant to qualitatively demonstrate the physical effect in medical
% PET/MR scanners.
%
% Input parameters:
% numOfPositron - number of positrons simulated;
% magneticField - the unitless magnetic field;
% annihilatedTreshold - annihilation threshold, [0,1], the higher the more
% likely to annihilate;
% scatterThreshold - scattering threshold, [0,1], the higher the more like-
% ly to scatter;
% vInitial - initial velocity distribution parameter;
% isTrack - track mode or point cloud mode; 0 - point cloud output, 1 -
% positron trajectory tracking.
%
% Output parameter:
% record - the record of the annihilated point or positron trajectory.
% record can either be a cell (containing matrice) or a matrix, depending
% on isTrack.
% The contents of record are in the format specified.
% The first three columns of the matrix are the corresponding coordinates.
% The next three columns of the matrix are the corresponding velocity.
% These are in the Cartesian coordinates.
% The last column records the number of scattering events.
%
% mengxiangxi_{a]_pku.edu.cn
% Created: 20200801
%
% Load and initiate the parameters
global NUM_OF_POSITRON MAGNETIC_FIELD ANNIHILATE_THRESHOLD ...
SCATTER_THRESHOLD V_INITIAL STEP_TIME ELECTRICITY_OF_POSITRON ...
MASS_OF_POSITRON;
NUM_OF_POSITRON = numOfPositron;
MAGNETIC_FIELD = magneticField;
ANNIHILATE_THRESHOLD = annihilateThreshold;
SCATTER_THRESHOLD = scatterThreshold;
V_INITIAL = vInitial;
STEP_TIME = 0.001;
ELECTRICITY_OF_POSITRON = 1;
MASS_OF_POSITRON = 1;
% Create the output parameter
if isTrack
record = cell(1);
else
record = zeros(NUM_OF_POSITRON, 7);
end
% Launch the positrons
for ii = 1:NUM_OF_POSITRON
p = Positron(V_INITIAL); % Specified in the class Positron
[p, track] = trackPositron(p,isTrack);
% Record the results
if isTrack
record{ii} = track; % record is a cell
else
record(ii,:) = [p.coord, p.dir, p.scatter]; % record is a matrix
end
end
end
function [p_out, track] = trackPositron(p_in, recordTrack)
% trackPositron
% Input Argument:
% p_in - the input positron, (class Positron);
% recordTrack - track trajectory or annihilation.
% Output Argument:
% p_out - the positron after propogation.
p = p_in;
terminate = 0; % termination flag
track = [p_in.coord, p_in.dir, p_in.scatter]; % Initiate using the input
while terminate == 0
p = propagate(p); % Positron propagate
p = scattering(p); % Positron scattering
if recordTrack == 1
track = [track; p.coord, p.dir, p.scatter];
end
terminate = isTerminate(p);
end
p_out = p;
end
function p_out = propagate(p_in)
% propagate
% Take care of the positron propagate.
% Input Argument:
% p_in - the input positron, (class Positron).
% Output Argument:
% p_out - the positron after propogation.
% The influence of Lorentz force on positron motion is considered in the
% propagation of positron.
% the default direction of magnetic is parallel to the Z-axis
global MAGNETIC_FIELD ELECTRICITY_OF_POSITRON MASS_OF_POSITRON STEP_TIME;
p = p_in;
if MAGNETIC_FIELD ~= 0
[dir_theta, dir_phi, dir_r] = p.getDirSph();
% the angular velocity of circular motion
w = MAGNETIC_FIELD * ELECTRICITY_OF_POSITRON / MASS_OF_POSITRON;
% the radius of circular motion
R = MASS_OF_POSITRON * dir_r * cos(dir_phi) / ...
ELECTRICITY_OF_POSITRON / MAGNETIC_FIELD;
v_x = dir_r * cos(dir_phi) * cos(dir_theta - w * STEP_TIME);
v_y = dir_r * cos(dir_phi) * sin(dir_theta - w * STEP_TIME);
v_z = p.dir(3);
x = R * (sin(w * STEP_TIME - dir_theta) + sin(dir_theta)) + ...
p.coord(1);
y = R * (cos(w * STEP_TIME - dir_theta) - cos(dir_theta)) + ...
p.coord(2);
z = p.dir(3) * STEP_TIME + p.coord(3);
p.coord = [x, y, z];
p.dir = [v_x, v_y, v_z];
else % No magnetic field, just propagate
for i = 1:3
p.coord(i) = p.dir(i) * STEP_TIME + p.coord(i);
end
end
p_out = p;
end
function p_out = scattering(p_in)
% scatter
% Scatter the positron
% Input Argument:
% p_in - the input positron, (class Positron).
% Output Argument:
% p_out - the positron after scattering.
global SCATTER_THRESHOLD;
p = p_in;
[dir_theta, dir_phi, dir_r] = p.getDirSph();
if rand() < SCATTER_THRESHOLD
% In this simplified model, the difference between spatial
% uniformity and time uniformity is ignored.
scat_ang_dev = asin(2 * rand() - 1);
scat_ang_azim = 2 * pi * rand();
% The Angle between the post-scatter velocity and the positive Y axis
% in the new coordinate system.
v_x_temp = dir_r * cos(scat_ang_dev);
v_y_temp = dir_r * sin(scat_ang_dev) * cos(scat_ang_azim);
v_z_temp = dir_r * sin(scat_ang_dev) * sin(scat_ang_azim);
v_x = (- v_z_temp * sin(dir_phi) + v_x_temp * cos(dir_phi)) * ...
cos(dir_theta) - v_y_temp * sin(dir_theta);
v_y = (- v_z_temp * sin(dir_phi) + v_x_temp * cos(dir_phi)) * ...
sin(dir_theta) - v_y_temp * cos(dir_theta);
v_z = v_z_temp * cos(dir_phi) + v_x_temp * sin(dir_phi);
p.dir = [v_x, v_y, v_z];
p.scatter = p.scatter + 1;
end
p_out = p;
end
function terminate_flag = isTerminate(p)
% isAnnihilate
% Determine if the positron is annihilated
% Output paramter:
% terminate_flag - 0 - not teminate, 1 - terminate
global ANNIHILATE_THRESHOLD;
[~, ~, v_p] = p.getDirSph();
annihilate_probab = atan(ANNIHILATE_THRESHOLD./v_p)./pi*2;
if rand() < annihilate_probab
% In this simplified model, the difference between spatial uniformity
% and time uniformity is ignored.
terminate_flag = 1;
else
terminate_flag = 0;
end
end