EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_base_geometry.h
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ base geometry headers
5 # Copyright (C) 2015 National Research Council Canada
6 #
7 # This file is part of EGSnrc.
8 #
9 # EGSnrc is free software: you can redistribute it and/or modify it under
10 # the terms of the GNU Affero General Public License as published by the
11 # Free Software Foundation, either version 3 of the License, or (at your
12 # option) any later version.
13 #
14 # EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
15 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
16 # FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for
17 # more details.
18 #
19 # You should have received a copy of the GNU Affero General Public License
20 # along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
21 #
22 ###############################################################################
23 #
24 # Author: Iwan Kawrakow, 2005
25 #
26 # Contributors: Frederic Tessier
27 # Blake Walters
28 # Marc Chamberland
29 # Reid Townson
30 # Ernesto Mainegra-Hing
31 # Hugo Bouchard
32 #
33 ###############################################################################
34 */
35 
36 
42 #ifndef EGS_BASE_GEOMETRY_
43 #define EGS_BASE_GEOMETRY_
44 
45 #include "egs_vector.h"
46 
47 #include <string>
48 #include <vector>
49 #include <iostream>
50 
51 using std::string;
52 using std::vector;
53 
54 class EGS_Application; // forward declaration
55 class EGS_Input;
57 
58 #ifdef BPROPERTY64
59  typedef EGS_I64 EGS_BPType;
60 #elif defined BPROPERTY32
61  typedef unsigned int EGS_BPType;
62 #elif defined BPROPERTY16
63  typedef unsigned short EGS_BPType;
64 #else
65  typedef unsigned char EGS_BPType;
66 #endif
67 
68 class label {
69 public:
70  string name;
71  vector<int> regions;
72 };
73 
94 
95 public:
96 
102  EGS_BaseGeometry(const string &Name);
103 
109  virtual ~EGS_BaseGeometry();
110 
124  inline bool isConvex() const {
125  return is_convex;
126  };
127 
133  virtual int inside(const EGS_Vector &x) = 0;
134 
141  virtual bool isInside(const EGS_Vector &x) = 0;
142 
150  virtual int isWhere(const EGS_Vector &x) = 0;
151 
160  static int findRegion(EGS_Float xp, int np, const EGS_Float *p) {
161  int ml = 0, mu = np;
162  while (mu - ml > 1) {
163  int mav = (ml+mu)/2;
164  if (xp <= p[mav]) {
165  mu = mav;
166  }
167  else {
168  ml = mav;
169  }
170  }
171  return mu - 1;
172  };
173 
204  virtual int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u,
205  EGS_Float &t, int *newmed=0, EGS_Vector *normal=0) = 0;
206 
219  virtual EGS_Float howfarToOutside(int ireg, const EGS_Vector &x,
220  const EGS_Vector &u);
221 
232  virtual EGS_Float hownear(int ireg, const EGS_Vector &x) = 0;
233 
239  virtual EGS_Float getVolume(int ireg) {
240  return 1.0;
241  }
242 
248  virtual EGS_Float getBound(int idir, int ind) {
249  return 0.0;
250  }
251 
257  virtual int getNRegDir(int idir) {
258  return 0;
259  }
260 
267  int regions() const {
268  return nreg;
269  };
270 
278  virtual bool isRealRegion(int ireg) const {
279  return (ireg >= 0 && ireg < nreg);
280  };
281 
288  virtual int medium(int ireg) const {
289  return region_media ? region_media[ireg] : med;
290  };
291 
297  virtual int getMaxStep() const {
298  return nreg+1;
299  };
300 
313  virtual int computeIntersections(int ireg, int n, const EGS_Vector &x,
314  const EGS_Vector &u, EGS_GeometryIntersections *isections);
315 
325  void setMedium(const string &Name);
326 
334  void setMedium(int start, int end, const string &Name, int delta=1);
335 
340  void setMedium(int imed) {
341  med = imed;
342  };
343 
351  void setMedium(int istart, int iend, int imed, int delta=1);
352 
363  void setMedia(EGS_Input *inp);
364 
372  static int nMedia();
373 
379  static const char *getMediumName(int ind);
380 
389  static int addMedium(const string &medname);
390 
396  static int getMediumIndex(const string &medname);
397 
401  inline bool hasRhoScaling() const {
402  return has_rho_scaling;
403  };
404 
408  virtual EGS_Float getRelativeRho(int ireg) const {
409  return rhor && ireg >= 0 && ireg < nreg ? rhor[ireg] : 1;
410  };
411 
417  virtual void setRelativeRho(int start, int end, EGS_Float rho);
418 
428  virtual void setRelativeRho(EGS_Input *);
429 
430  EGS_Float getMediumRho(int ind) const;
431 
432  virtual void setApplication(EGS_Application *app);
433 
436  inline bool hasBScaling() const {
437  return (has_B_scaling || has_Ref_rho);
438  };
439 
442  virtual EGS_Float getBScaling(int ireg) const {
443  if (has_Ref_rho && has_B_scaling) {
444  if (bfactor && ireg >= 0 && ireg < nreg) {
445  return getMediumRho(medium(ireg))/rhoRef*bfactor[ireg];
446  }
447  else {
448  return 1.0;
449  }
450  }
451  else if (has_Ref_rho && !has_B_scaling) {
452  if (ireg >= 0 && ireg < nreg) {
453  return getMediumRho(medium(ireg))/rhoRef;
454  }
455  else {
456  return 1.0;
457  }
458  }
459  else if (!has_Ref_rho && has_B_scaling) {
460  if (bfactor && ireg >= 0 && ireg < nreg) {
461  return bfactor[ireg];
462 
463  }
464  else {
465  return 1.0;
466  }
467  }
468  else {
469  return 1.0;
470  }
471  }
472 
478  virtual void setBScaling(int start, int end, EGS_Float bf);
479 
489  virtual void setBScaling(EGS_Input *);
490 
496  const string &getName() const {
497  return name;
498  };
499 
506  virtual const string &getType() const = 0;
507 
527 
544  static EGS_BaseGeometry *createSingleGeometry(EGS_Input *inp);
545 
553  static void clearGeometries();
554 
560  void setDebug(bool deb) {
561  debug = deb;
562  };
563 
570  static EGS_BaseGeometry *getGeometry(const string &Name);
571 
572  static EGS_BaseGeometry **getGeometries();
573 
574  static int getNGeometries();
575 
582  static string getUniqueName();
583 
592  void setName(EGS_Input *inp);
593 
602  void setBoundaryTolerance(EGS_Input *inp);
603 
606  void setBoundaryTolerance(EGS_Float tol) {
607  boundaryTolerance = tol;
608  halfBoundaryTolerance = tol/2.;
609  }
610 
613  virtual bool hasBooleanProperty(int ireg, EGS_BPType prop) const {
614  if (!bp_array) {
615  return (prop & bproperty);
616  }
617  return ireg >= 0 && ireg < nreg ? prop & bp_array[ireg] : false;
618  };
619 
625  virtual void setBooleanProperty(EGS_BPType prop);
626 
634  virtual void addBooleanProperty(int bit);
635 
642  virtual void setBooleanProperty(EGS_BPType prop, int start, int end,
643  int step=1);
644 
651  virtual void addBooleanProperty(int bit, int start, int end, int step=1);
652 
663  virtual void printInfo() const;
664 
670  static void describeGeometries();
671 
678  inline int ref() {
679  return ++nref;
680  };
681 
682 
690  inline int deref() {
691  return --nref;
692  };
693 
697  static void setActiveGeometryList(int list);
698 
699  static int getLastError() {
700  return error_flag;
701  };
702 
703  static void resetErrorFlag() {
704  error_flag = 0;
705  };
706 
708  EGS_Float getBoundaryTolerance() {
709  return boundaryTolerance;
710  };
711 
713  virtual void getNumberRegions(const string &str, vector<int> &regs);
714 
716  virtual void getLabelRegions(const string &str, vector<int> &regs);
717 
719  virtual const string &getLabelName(const int i) {
720  return labels[i].name;
721  }
722 
724  virtual int getLabelCount() {
725  return labels.size();
726  }
727 
729  int setLabels(EGS_Input *input);
730 
732  int setLabels(const string &inp);
733 
734 protected:
735 
741  int nreg;
742 
748  string name;
749 
756  short *region_media;
757 
764  int med;
765 
766 
771 
775  EGS_Float *rhor;
776 
780  bool has_B_scaling, has_Ref_rho;
781 
785  EGS_Float *bfactor;
786 
790  EGS_Float rhoRef;
791 
796  int nref;
797 
806  virtual void setMedia(EGS_Input *inp, int nmed, const int *med_ind);
807 
812  bool debug;
813 
822  bool is_convex;
823 
828  EGS_BPType bproperty;
829 
835  EGS_BPType *bp_array;
836 
838  EGS_Float boundaryTolerance, halfBoundaryTolerance;
839 
841  static int error_flag;
842 
858  vector<label> labels;
859 
862 
863 private:
864 
874  static int active_glist;
875 
876 };
877 
879  EGS_Float t;
880  EGS_Float rhof;
881  EGS_I32 ireg;
882  EGS_I32 imed;
883 };
884 
1130 /* \example geometry/example1/geometry_example1.cpp
1131 
1132  Suppose that you frequently use the same complex geometry and you
1133  are tired of always having to include the long definition of this geometry
1134  into your input file. In this case it may be worthwhile to implement your
1135  own geometry DSO that constructs the geometry of interest with very little
1136  input (or even no input at all). This example illustrates how to do this
1137  programatically.
1138 
1139  As this is just an example, the geometry we will be constructing is not very
1140  complex: it consists of a cylinder inscribed in a sphere and the sphere
1141  inscribed in a box as shown in the figure ???. We want to be able to define
1142  the media of the cylinder, the sphere and the box in the input file
1143  (so that we can easily use different PEGS data sets, for example). All other
1144  dimensions of our example geometry are fixed.
1145 
1146  To create our new geometry, we make a new subdirectory in the egspp
1147  geometries subdirectory called \c example1. To be able to use the
1148  predefined rules for building geometry libraries we create a dummy
1149  header file \c %geometry_example1.h that does nothing else except to
1150  define the \c EGS_EXAMPLE1_DLL
1151  macro:
1152  \include geometry/example1/geometry_example1.h
1153  and put the implementation in \c %geometry_example1.cpp.
1154  \dontinclude geometry/example1/geometry_example1.cpp
1155  In this file we include the just created header file and the EGS_Input
1156  class header file to get access to its definition,
1157  needed in the \c %createGeometry() function (see below):
1158  \until egs_input
1159  We also need the declarations of the various geometry classes that we
1160  will use to construct our geometry:
1161  \until egs_envelope_geometry.h
1162  We then define the dimensions of the various geometrical structures
1163  \until box_size
1164 
1165  The remaining task is to implement a C-style \c %createGeometry() function
1166  exported by the DSO
1167  \until EGS_EXAMPLE1_EXPORT
1168  We check that the input object is valid and return \c null if it is not
1169  \until }
1170  We then look for a definition of the cylinder medium and print a warning if
1171  such a definition is missing:
1172  \skipline cyl_medium
1173  \until }
1174  We repeat basically the same code to obtain the sphere and the box
1175  media
1176  \until }
1177  \until }
1178  and return \c null if any of the media definition was missing:
1179  \until if
1180  Now we are ready to construct our geometry. We first create a set of
1181  parallel z-planes
1182  \until EGS_ZProjector
1183  As the planes will only be used by our geometry object internally, we don't
1184  care how the plane set is named and therefore give it a name using
1185  the static function EGS_BaseGeometry::getUniqueName().
1186  In a similar way we create a set of z-cylinders consisting of a single
1187  cylinder
1188  \until ZProjector
1189  and make a closed cylinder modeled as a 2D-geometry from the planes and the
1190  cylinder surface:
1191  \until getUniqueName
1192  We then set the medium of the cylinder
1193  \until setMedium
1194  The next step is to create the sphere as a set of spheres consisting of
1195  a single sphere
1196  \until getUniqueName
1197  and to fill it with the sphere medium defined in the input
1198  \until setMedium
1199  We then use an envelope geometry to inscribe the cylinder into the just
1200  created sphere:
1201  \until getUniqueName
1202  The next step is to create the outer box and fill it with the box medium
1203  defined in the input
1204  \until setMedium
1205  (as our geometry will always be centered about the origin, we pass a
1206  \c null affine transformation to the constructor of the box).
1207  We now inscribe the \c the_sphere geometry object (which is the sphere
1208  with the cylinder inscribed in it) into the just created box using
1209  again an envelope geometry,
1210  \until the_box
1211  set the name of \c the_box from the input using the inherited
1212  \link EGS_BaseGeometry::setName(EGS_Input*) setName() \endlink method,
1213  \until setName
1214  and return a pointer to the just created geometry
1215  \until }
1216  \until }
1217  That's all.
1218 
1219  In the Makefile we include the egspp make config files,
1220  \dontinclude geometry/example1/Makefile
1221  \skipline EGS_CONFIG
1222  \until my_machine
1223  set the C-preprocessor defines to the standard set of defines
1224  for our egspp configuration,
1225  \until DEFS
1226  define the name of our DSO,
1227  \until library
1228  and and the set of files needed to build the DSO
1229  \until lib_files
1230  Because we are using planes, cylinders, spheres, a ND-geometry, a box
1231  geometry and an envelope geometry to construct our geometry, we must
1232  link against the geometry libraries containing these geometries.
1233  We accomplish this by seting the \c link2_libs make variable to the
1234  list of needed libraries
1235  \until egs_genvelope
1236  We then include the standard rules for building a DSO
1237  \until include
1238  and set the dependencies of our object files
1239  \until make_depend
1240  We must also add the dependencies on the various geometry header files
1241  that we are including
1242  \until egs_envelope_geometry.h
1243  We can now build our geometry DSO by typing \c make and use
1244  it by including input like this in the input file
1245  \verbatim
1246  :start geometry:
1247  library = geometry_example1
1248  name = some_name
1249  cylinder medium = H2O521ICRU
1250  sphere medium = AL521ICRU
1251  box medium = AIR521ICRU
1252  :stop geometry:
1253  \endverbatim
1254  It is worth noting that we could have constructed this geometry as the
1255  \link EGS_UnionGeometry union \endlink of the cylinder pointed to
1256  by \c the_cyl, the sphere pointed to by \c sphere and the
1257  box pointed to by \c box. However, such an implementation would be
1258  much slower at run time because for each invocation of a geometry
1259  method the union will have to interogate the geometry methods of all
1260  3 objects. In contrast, the envelope geometry will be checking only
1261  against the box for points outside of the box, the box and the sphere for
1262  points inside the box but outside the sphere, the sphere and the cylinder
1263  for points inside the sphere but outside the cylinder and the cylinder
1264  only for points inside the cylinder.
1265 
1266  The complete source code of this example geometry is below
1267  \include geometry/example1/Makefile
1268  \include geometry/example1/geometry_example1.h
1269 
1270 */
1271 
1272 /* \example geometry/example2/geometry_example2.cpp
1273 
1274  This example illustrates how to create the same geometry as
1275  described in
1276  <a href="geometry_2example1_2geometry__example1_8cpp-example.html">
1277  this example</a> but using EGS_BaseGeometry::createSingleGeometry
1278  to create the geometry from definitions stored in an EGS_Input objects.
1279  The header file is similar as in the previous
1280  <a href="geometry_2example1_2geometry__example1_8cpp-example.html">
1281  example</a>, but unlike in the previous example we don't need the
1282  header files of the various geometries being used.
1283  Our implementation consists of a single C-style function
1284  \c %createGeometry()
1285  \dontinclude geometry/example2/geometry_example2.cpp
1286  \skipline createGeometry
1287  As with the previous
1288  <a href="geometry_2example1_2geometry__example1_8cpp-example.html">
1289  example</a>, we check for valid input and extract the
1290  media of the cylinder, sphere and the box from the input
1291  \until !ok
1292  We will create the various geometry objects needed in our composite
1293  geometry by passing EGS_Input objects containing geometry definitions
1294  to the EGS_BaseGeometry::createSingleGeometry function.
1295  We start with the set of planes needed to close the cylinder.
1296  We first create an EGS_Input object named \c geometry that we will
1297  be passing to the geometry creation function:
1298  \until plane_input
1299  We then create properties for the library name,
1300  \until plane_library
1301  the plane-set type,
1302  \until plane_type
1303  the name of the set of planes,
1304  \until plane_name
1305  and the plane positions
1306  \until plane_positions
1307  and add them to the geometry definition
1308  \until plane_positions
1309  We then create the set of planes
1310  \until createSingleGeometry
1311  The definition and creation of the set of cylinders is very similar
1312  \until createSingleGeometry
1313  We now make a closed cylinder as a 2D geometry
1314  \until createSingleGeometry
1315  and set its medium
1316  \until setMedium
1317  In a similar fashion we create the sphere, use an envelope geometry
1318  to inscribe the cylinder into the sphere, create a box, and use again
1319  an envelope to inscribe the sphere into the box. The code is not shown
1320  here as it is lengthy and boring but can be found at the end of this page.
1321  Finally we set the name of the newly created geometry \c the_box and
1322  return it:
1323  \skipline the_box->setName
1324  \until }
1325  \until }
1326 
1327  The Makefile is now simpler as we don't need to link against the
1328  various geometry DSOs as they are loaded dynamically at run time
1329  by the EGS_BaseGeometry::createSingleGeometry function and
1330  our implementation also does not depend on the header files of the
1331  various geometry classes being used:
1332  \include geometry/example2/Makefile
1333 
1334  For completeness here is the header file:
1335  \include geometry/example2/geometry_example2.h
1336  and the complete implementation:
1337 */
1338 
1339 
1340 #endif
#define EGS_EXPORT
Export symbols from the egspp library.
Definition: egs_libconfig.h:90
virtual bool hasBooleanProperty(int ireg, EGS_BPType prop) const
Is the boolean property prop set for region ireg ?
int regions() const
Returns the number of local regions in this geometry.
int deref()
Decrease the reference count to this geometry.
virtual int getNRegDir(int idir)
virtual int medium(int ireg) const
Returns the medium index in region ireg.
EGS_Float t
distance to next region boundary
void setMedium(int imed)
Set all regions to a medium with index imed.
EGS_Float * bfactor
Array with B field scaling factors.
EGS_Float rhoRef
Reference density for B field scaling.
EGS_Float rhof
relative mass density in that region
EGS_BPType bproperty
A bit mask of boolean properties for the entire geometry.
string name
Name of this geometry.
EGS_BPType * bp_array
An array of boolean properties on a region by region basis.
EGS_Vector methods for the manipulation of 3D vectors in cartesian co-ordinates.
virtual int getMaxStep() const
Returns the maximum number of steps through the geometry.
A class representing 3D vectors.
Definition: egs_vector.h:56
virtual EGS_Float getVolume(int ireg)
Calculates the volume of region ireg.
bool hasRhoScaling() const
Does this geometry object have a mass density scaling feature?
int nref
Number of references to this geometry.
static int findRegion(EGS_Float xp, int np, const EGS_Float *p)
Find the bin to which xp belongs, given np bin edges p.
EGS_I32 imed
medium index
int med
Medium index.
EGS_GLIB_EXPORT EGS_BaseGeometry * createGeometry(EGS_Input *input)
Definition: egs_glib.cpp:57
vector< label > labels
Labels.
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
virtual EGS_Float getBound(int idir, int ind)
Returns region boundaries in direction determined by idir.
EGS_Float getBoundaryTolerance()
Get the value of the boundary tolerance.
bool isConvex() const
Is the geometry convex?
virtual EGS_Float getBScaling(int ireg) const
Get the B field scaling factor in region ireg.
virtual const string & getLabelName(const int i)
Get the name of the i-th explicit label in the geometry.
bool is_convex
Is this geometry convex?
void setBoundaryTolerance(EGS_Float tol)
Set the value of the boundary tolerance from argument.
bool hasBScaling() const
Does this geometry object have a B field scaling feature?
virtual EGS_Float getRelativeRho(int ireg) const
Get the relative mass density in region ireg.
EGS_Float boundaryTolerance
Boundary tolerance for geometries that need it.
virtual int getLabelCount()
Get the number of explicit labels in the geometry.
short * region_media
Array of media indeces.
bool has_rho_scaling
Does this geometry have relative mass density scvaling?
int ref()
Increase the reference count to this geometry.
virtual bool isRealRegion(int ireg) const
Returnes true if ireg is a real region, false otherwise.
static int error_flag
Set to non-zero status if a geometry problem is encountered.
bool debug
Debugging flag.
EGS_I32 ireg
region index
A class for storing information in a tree-like structure of key-value pairs. This class is used throu...
Definition: egs_input.h:182
const string & getName() const
Get the name of this geometry.
EGS_Application * app
The application this object belongs to.
void setDebug(bool deb)
Turn debugging on.
int nreg
Number of local regions in this geometry.
EGS_Float * rhor
Array with relative mass densities.
Base class for advanced EGSnrc C++ applications.
bool has_B_scaling
Does this geometry has B field scaling factor?