EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_rz.cpp
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ rz geometry
5 # Copyright (C) 2016 Randle E. P. Taylor, Rowan M. Thomson,
6 # Marc J. P. Chamberland, Dave W. O. Rogers
7 #
8 # This file is part of EGSnrc.
9 #
10 # EGSnrc is free software: you can redistribute it and/or modify it under
11 # the terms of the GNU Affero General Public License as published by the
12 # Free Software Foundation, either version 3 of the License, or (at your
13 # option) any later version.
14 #
15 # EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
16 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 # FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for
18 # more details.
19 #
20 # You should have received a copy of the GNU Affero General Public License
21 # along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
22 #
23 ###############################################################################
24 #
25 # Author: Randle Taylor, 2016
26 #
27 # Contributors: Marc Chamberland
28 # Rowan Thomson
29 # Dave Rogers
30 # Martin Martinov
31 #
32 ###############################################################################
33 #
34 # egs_rz was developed for the Carleton Laboratory for Radiotherapy Physics.
35 #
36 ###############################################################################
37 */
38 
39 
45 #include <map>
46 #include "egs_input.h"
47 #include "egs_rz.h"
48 #include "../egs_cylinders/egs_cylinders.h"
49 #include "../egs_planes/egs_planes.h"
50 
51 
52 string EGS_RZGeometry::RZType = "EGS_RZ";
53 
54 EGS_RZGeometry::EGS_RZGeometry(vector<EGS_BaseGeometry *> geoms,
55  vector<EGS_Float> rads, vector<EGS_Float> zbs, const string &name) :
56  EGS_NDGeometry(geoms, name), radii(rads), zbounds(zbs) {
57 
58  /* we always set inner most radii to 0 (simplifies volume calcs) */
59  if (radii[0] != 0) {
60  radii.insert(radii.begin(), 0.);
61  }
62 
63  vector<EGS_Float> vol;
64  for (size_t r=0; r < radii.size()-1; r++) {
65 
66  EGS_Float rmin = radii[r];
67  EGS_Float rmax = radii[r+1];
68 
69  EGS_Float area = M_PI*(rmax*rmax - rmin*rmin);
70 
71  for (size_t plane = 0; plane < zbounds.size()-1; plane++) {
72  EGS_Float zmin = zbounds[plane];
73  EGS_Float zmax = zbounds[plane+1];
74  reg_vol.push_back((zmax-zmin)*area);
75  }
76  }
77 
78 };
79 
80 EGS_Float EGS_RZGeometry::getBound(int idir, int ind) {
81  if (idir == ZDIR && ind >= 0 && ind < (int)zbounds.size()) {
82  return zbounds[ind];
83  }
84  else if (idir == RDIR && ind >= 0 && ind < (int)radii.size()) {
85  return radii[ind];
86  }
87  return 0.;
88 }
89 
90 
92  if (dir == ZDIR) {
93  return zbounds.size() - 1;
94  }
95  else if (dir == RDIR) {
96  return radii.size() - 1;
97  }
98  return 0;
99 }
100 
101 EGS_Float EGS_RZGeometry::getVolume(int ireg) {
102  if (ireg < 0 || ireg >= nreg) {
103  return -1;
104  }
105  return reg_vol[ireg];
106 }
107 
108 
109 /* Initialization helpers for EGS_RZ */
110 namespace egs_rz {
111 
112 vector<EGS_Float> getRadiiByShells(EGS_Input *input) {
113 
114  vector<EGS_Float> radii;
115 
116  vector<int> nshells;
117  int err = input->getInput("number of shells", nshells);
118  if (err) {
119  return radii;
120  }
121 
122  vector<EGS_Float> thick;
123  err = input->getInput("shell thickness", thick);
124  if (err) {
125  return radii;
126  }
127 
128  EGS_Float cur_r = 0;
129  for (size_t shell_group=0; shell_group < min(thick.size(), nshells.size()); shell_group++) {
130 
131  for (int shell = 0; shell < nshells[shell_group]; shell++) {
132  cur_r += thick[shell_group];
133  radii.push_back(cur_r);
134 
135  }
136  }
137 
138  return radii;
139 
140 }
141 
142 vector<EGS_Float> getRadii(EGS_Input *input) {
143 
144  vector<EGS_Float> radii;
145  int err = input->getInput("radii", radii);
146  if (err) {
147  /* no explicit radii input so assume user will be specifying using shells */
148  radii = getRadiiByShells(input);
149  }
150 
151  return radii;
152 
153 }
154 
155 
156 vector<EGS_Float> EGS_RZ_LOCAL getZPlanesBySlabs(EGS_Input *input) {
157 
158  vector<EGS_Float> zplanes;
159 
160  EGS_Float zo = 0;
161 
162  int err = input->getInput("first plane", zo);
163  if (err) {
164  egsWarning("RZ: missing 'first plane' input. Assuming zo=0");
165  zo = 0;
166  }
167 
168  vector<int> nslabs;
169  err = input->getInput("number of slabs", nslabs);
170  if (err) {
171  return zplanes;
172  }
173 
174  vector<EGS_Float> thick;
175  err = input->getInput("slab thickness", thick);
176  if (err) {
177  return zplanes;
178  }
179 
180  EGS_Float cur_z = zo;
181  zplanes.push_back(zo);
182  for (size_t slab_group=0; slab_group < min(thick.size(), nslabs.size()); slab_group++) {
183 
184  for (int slab = 0; slab < nslabs[slab_group]; slab++) {
185  cur_z += thick[slab_group];
186  zplanes.push_back(cur_z);
187  }
188  }
189 
190  return zplanes;
191 
192 }
193 
194 vector<EGS_Float> EGS_RZ_LOCAL getZPlanes(EGS_Input *input) {
195 
196  vector<EGS_Float> zplanes;
197  int err = input->getInput("z-planes", zplanes);
198  if (err) {
199  zplanes = getZPlanesBySlabs(input);
200  }
201 
202  return zplanes;
203 
204 }
205 
206 bool allIncreasing(vector<EGS_Float> vec) {
207 
208  if (vec.size() == 0) {
209  return true;
210  }
211 
212  EGS_Float last = vec[0];
213  for (size_t i=1; i < vec.size(); i++) {
214  if (vec[i] <= last) {
215  return false;
216  }
217  last = vec[i];
218  }
219 
220  return true;
221 
222 }
223 
224 
225 };
226 
227 extern "C" {
228 
229  EGS_RZ_EXPORT EGS_BaseGeometry *createGeometry(EGS_Input *input) {
230 
231  if (!input) {
232  egsWarning("createGeometry(egs_rz): null input?\n");
233  return 0;
234  }
235 
236  vector<EGS_Float> radii = egs_rz::getRadii(input);
237  if (radii.size() == 0) {
238  egsWarning("createGeometry(rz): wrong/missing radii inputs\n");
239  return 0;
240  }
241  if (!egs_rz::allIncreasing(radii)) {
242  egsWarning("createGeometry(rz): radii must be monotonically increasing\n");
243  return 0;
244  }
245 
246  vector<EGS_Float> zplanes = egs_rz::getZPlanes(input);
247  if (zplanes.size() == 0) {
248  egsWarning("createGeometry(rz): wrong/missing z plane inputs\n");
249  return 0;
250  }
251 
253  EGS_PlanesZ *planes = new EGS_PlanesZ(zplanes, EGS_BaseGeometry::getUniqueName(), EGS_ZProjector("z-planes"));
254 
255  vector<EGS_BaseGeometry *> rz_geoms;
256  rz_geoms.push_back(planes);
257  rz_geoms.push_back(cyl);
258 
259  EGS_BaseGeometry *rz = new EGS_RZGeometry(rz_geoms, radii, zplanes);
260 
261  if (!rz) {
262  egsWarning("createGeometry(rz): failed to create nd geometry\n");
263  return 0;
264  }
265 
266  rz->setName(input);
267  rz->setMedia(input);
268  rz->setLabels(input);
269 
270  return rz;
271 
272  }
273 
274 }
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
void setMedia(EGS_Input *inp)
Set the media in the geometry from the input pointed to by inp.
int nreg
Number of local regions in this geometry.
void setName(EGS_Input *inp)
Set the name of the geometry from the input inp.
static string getUniqueName()
Get a unique geometry name.
int setLabels(EGS_Input *input)
Set the labels from an input block.
A set of concentric cylinders.
A class for storing information in a tree-like structure of key-value pairs. This class is used throu...
Definition: egs_input.h:182
int getInput(const string &key, vector< string > &values) const
Assign values to an array of strings from an input identified by key.
Definition: egs_input.cpp:338
A class modeling a N-dimensional geometry.
A set of parallel planes.
Definition: egs_planes.h:167
a subclass of EGS_NDGeometry for conveniently defining an RZ geometry
Definition: egs_rz.h:116
EGS_Float getBound(int idir, int ind)
get RZ boundary for a given direction and directional index
Definition: egs_rz.cpp:80
EGS_Float getVolume(int ireg)
get mass of a given region
Definition: egs_rz.cpp:101
EGS_RZGeometry(vector< EGS_BaseGeometry * > geoms, vector< EGS_Float > rads, vector< EGS_Float > zbs, const string &name="")
RZ geometry constructor.
Definition: egs_rz.cpp:54
int getNRegDir(int dir)
get number of regions in a given direction
Definition: egs_rz.cpp:91
A class representing 3D vectors.
Definition: egs_vector.h:57
A projector into the z-plane.
EGS_GLIB_EXPORT EGS_BaseGeometry * createGeometry(EGS_Input *input)
Definition: egs_glib.cpp:57
EGS_Input class header file.
An egs_nd_geometry wrapper to simplify RZ geometry creation.
EGS_InfoFunction EGS_EXPORT egsWarning
Always use this function for reporting warnings.