Merge branch develop into master for 2023 release #965

Remove superfluous command in entry point source_check_emin copied from entry point source_get_emax during implementation, which was setting the value for emax. Not really clear to me why this was happening but it turns out that having emax as an argument in one entry point would force to also have in the other entry point to avoid SEGFAULT errors. However, line has been removed as it shouldn't have been there in the first place!

Before this change, there was a small risk of buffer overflow in the tutor7pp fluence scoring arrays. The value r2 was checked to be less than 400, and then the scoring bin number (from 0 to 199) is calculated using the expression (int)(sqrt(r2)*10.). The problem with this check is that certain double precision values less than 400 can nevertheless have a square root equal to exactly 20.0 due to floating-point rounding. For example, taking the square root of the double 399.99999999999994 results in 20.0 under certain conditions, not a value slightly less than 20.0 as would be expected. The default g++ rounding mode is round-to-nearest, so even with the square root's true value being under 20.0, the resulting double can be rounded upwards to 20.0 if the true value is closer to 20.0 than the previous representable float. Then, scoring in bin 200 results in buffer overflow in the scoring array. This change ensures the resulting scoring array bin integer is in bounds (< 200), eliminating the potential for buffer overflow.

Check that the current region index (ir) is not -1 for region dose scoring. Before this change, there was an out-of-bounds read in d_reg_index when the region was -1 (outside the geometry).

Vacuum has a special medium number -1, while all other media get offsets from 0 to nmedia. Before this change, there was a buffer overflow if media dose scoring was used with a geometry that had vacuum as one of the media. The vector massM is indexed by medium number, but for any vacuum regions, it would read and write memory at the memory location one before the buffer bounds: massM[-1]. Skipping the vacuum region masses shouldn't matter, as vacuum dose is always zero. massM is only used later on with non-vacuum media, so skipping vacuum regions is the right choice instead of an alternative like appending the vacuum mass to massM.

Bug: results depend on the order in which the material composition is defined. Reported by Iymad Mansour who observed that egs_brachy simulations of an irradiation with an isotropic 35 keV photon point source, embedded in a steel applicator, produced 5% differences in the surface dose of a water phantom when the order of the steel components was changed. Reproduced with egs_kerma for a 35 keV photon point source surrounded by a 1 mm steel encapsulation embedded in a spherical water phantom. Kerma at the spherical water shell next to the encapsulation differs by 4%. Cause: sampling element and atomic shell for photo-electric interaction with MCDF-XCOM shell-wise using unsorted array zelem rather than pe_zsorted. Solution: using sorted array pe_zsorted fixes the issue.

Fix a mistake in the previous fix for non-alpha numeric characters in the job name in egs-parallel jobs resulting in a crash. The previous fix removed non-alphanum characters from the start of the original job name, but actually we needed to remove from the start of the trimmed 14 character job name.

Change a few occurrences of endl to \n in egs++ scoring routines, to increase performance by about 50%. The reason for this is that endl performs a buffer flush each time, which is not necessary. The performance hit was very noticeable for large arrays, like a 300x300x300 3ddose file.

Add a record of the git hash, branch and compile time to the c++ makefiles. This allows the c++ codes to access those parameters as defines. Advanced applications now output those values during initEGSnrcBackEnd().

Add the compilation time of an application to the printout of git hash and branch, and add these messages to the mortran codes. Also add a check for applications in the makefile to compare with the HEN_HOUSE versions and print a warning if there are differences. With the git information now printed to machine.mortran, it is not necessary to pass it as c++ compile flags. However, I will leave these flags for future use.

Also replace an old reference to nmake with mingw32-make.

Change the git branch and commit reporting to be passed as defines at compile time, instead of being generated during the configuration of EGSnrc. This fixes the issue where users could switch branches but not run the configuration script, resulting in the wrong git branch and commit being output.

As discussed in issue #745 ensure that one EGS_TrackScoring object is define to prevent bugs at write time when multiple objects try to write to the same file.

Fix memory leaks in EGS_AdvancedApplication: interpolation arrays are never freed, and transportp is allocated by EGS_Input::getInputItem but never freed. Also Fix memory leak in EGS_ObjectFactory::createObjects: EGS_Input* input is allocated by EGS_Input::takeInputItem but never freed. Memory leaks reported by AddressSanitizer.

After compiling egs_view with asan support, there was a buffer overflow error traced back to the memcpy call in ImageWindow::paintEvent. ERROR: AddressSanitizer: dynamic-stack-buffer-overflow ... READ of size 400 at 0x7ffe71889e40 thread T0 #0 0x7f543595ecdf (/usr/lib/x86_64-linux-gnu/libasan.so.5+0x99cdf) #1 0x55dc3acdbb15 in memcpy #2 0x55dc3acdbb15 in ImageWindow::paintEvent(QPaintEvent*) #3 0x7f5434f86047 in QWidget::event(QEvent*) The issue is that memcpy will always copy sizeof(lastRegions) bytes into the array regions. But before this change, regions could be shorter than lastRegions, leading to a buffer overflow. After this change, maxreg is always set to N_REG_MAX, the length of lastRegions.

@mjukdisk

In the solution to one of the questions in the getting-started manual, the total dose was calculated as the sum of dose to regions, instead of appropriately weighted by the relative volume of each region. This was corrected, as well as another typo in the instructions. Reported by @mjukdisk in discussion #780.

Fix a GUI error that occurred in DYNJAWS jaw definitions when clicking the button to calculate the jaw settings from the field size and SSD. The button triggered a function that was duplicated with the same name in SYNCJAWS, leading to the error. As the fix, the syncjaws version of the function has been renamed. This bug has existed since SYNCJAWS was added to the GUI.

Set the minimum distance between component modules in BEAMnrc to 1e-5 cm (macro $MIN_GAP). Previously this minimum gap was set to 0.01 cm. Given that the code is normally used in double-precision nowadays, there is no point in leaving such a large value as the default. This change follows for looking into issue #773.

Add the source-axis distance to the z position of the particle before it is written to the phase space file. This change only affects IAEA phsp sources, and only when scoring particles in BEAMnrc coordinates.

Fix a minor bug in FLURZnrc where it expected the wrong number of inputs when reading in the SLOTE parameter from the egsinp file. This did not result in any actual issues, but an error message was printed to the screen about the user providing the wrong number of inputs.

Implement a graceful exit from the addphsp application if the total number of particles is greater than 2^32-1, the maximum count that can be stored in the header. Only enforce this for a phase space in EGSnrc format.

Restore some code from DYNVMLC that seemed necessary to avoid changing regions when the particle had not actually traveled as far as the region boundary.

Note that commented out debug write statements are still present in the source code.

Also comment out debug write statements.

Add $SYNCHDMLC_ORIENT options to flip leaf numbering and cross-sections with respect to the axis perpendicular to the opening axis, so that the MLC coordinate system matches that in the manufacturer (Varian) specs for 120MLC and HD120. This reflection reverses the order of leaf numbering, hence leaf opening sequences are automatically "flipped" so that the opening coordinate applied to leaf i is instead applied to leaf N-i+1, where N is the total number of leaves. See the updated BEAMnrc manual for more details. Also, add new ORIENT=2,3 options to SYNCHDMLC gui

Reverse leaf opening sequence definition for option ORIENT_$SYNCHDMLC=2,3

Reorder cross section y boundaries when input restrictions are relaxed.

Add an egs++ geometry called egs_mesh that implements particle transport methods for imported CAD mesh files in the Gmsh format. Note that the mesh filename extension must be ".msh".

This is a combination of 8 original commits to start implementing a .msh file parser, and associated unit tests: - parse msh file version header - parse msh2 $Nodes section - parse msh2 and msh4 physical groups - parse msh2 $Elements section - parse msh4 $Entities section

The original intent for the msh file parser was to support both the legacy (but stable) v2.2 and today's default (v4.1). Thinking again, it's better if we only need to maintain a single parser. v4.1 is the default for Gmsh today. There are no planned changes to the sections we care about ($Entities, $PhysicalNames, $Elements and $Nodes). Any future 4.x revisions will hopefully not break the parser and we could "support" them by accepting the new versions in parse_msh_version.

This is a combination of 11 original commits related to the implementing the msh4 parser, and associated unit tests: - parse msh4 $Nodes entity bloc - parse entire msh4 $Nodes section - parse single msh4 element bloc - parse entire msh4 $Elements section - add default values for mesh data structs - Use namespace constant SIZET_MAX for default value - add test for duplicate msh4 entities - fix container size warnings - add minimum required input sections for msh v4.1

This is a combination of 14 original commits related to catching errors and adding exception throws in the msh4 parser: - catch node tag overflow - catch duplicate tags for all sections - bubble up section parsing errors - document that parsing functions ensure unique tags - check that each mesh volume has a known physical group tag - check that each tetrahedron has a valid volume (entity) tag - throw exceptions for errors in parse_msh_version - throw exceptions for errors in parse_msh4_groups - throw exceptions for errors in parse_msh4_node_bloc - throw exceptions for errors in parse_msh4_nodes - throw exceptions for errors in parse_msh4_element_bloc - throw exceptions for errors in parse_msh4_elements - throw exceptions for errors in parse_msh4_entities - throw exceptions for errors in parse_msh_file and parse_msh4_body

Remove all files not directly related to msh parsing to reduce the scope of this branch. Also, temporarily remove egs_mesh library from egs++ Makefile: it will be added back when egs_mesh.h and egs_mesh.cpp are added back as well.

Add msh_parser namespaces to split the API into public and private sections. Also, remove the redundant msh4 label in names of functions now in the msh41 namespace, and arrange data types at top of msh41 namespace.

Add pass-through constructors for all plain-old-data (POD) mesh data types. We previously relied on aggregate initialization for brevity but that is not actually supported in C++11 if structs have default members (it was added in C++14). The new constructors are boilerplate but keep us compliant with C++11.

Added static to all functions in msh_parser so it can be included in multiple translation units without causing linker errors. This should probably be refactored into a traditional header-implementation pair.

We're using the standard Möller-Trumbore ray-triangle intersection algorithm for each side of the tetrahedron. 1. Check if the velocity is parallel to the triangle, exit if so. 2. Calculate the barycentric coordinates u, v of the potential intersection point. 3. If the barycentric coordinates are inside the triangle, return the distance to the intersection point. Both clockwise and counter-clockwise point windings are allowed.

This is a combination of 4 original commits related to node ordering and storing the neighbour face boundaries: 1. Update neighbour algorithm to use int for elements: they must fit into an int to work with EGSnrc. 2. Order the nodes consistently ABCD. 3. Guarantee neighbour-face node ordering: To help implement howfar and hownear, guarantee that face 0 is made up of nodes b, c, d, face 1 is nodes a, c, d. This is an invariant we have to maintain for each mesh tetrahedron to avoid subtle transport errors. 4. Store which faces are boundaries: This will ensure only boundary faces are considered for intersection with an exterior point, assuming robustness in the intersection calculation.

Loop over the boundary faces and find the closest intersection along the given ray.

By rejecting negative distances we try to avoid negative ustep errors encountered with the initial intersection algorithm.

We'd like to write an output file based on the mesh when the simulation is done. This is most convenient to do in something similar to Mevegs_Application::outputResults but requires specific knowledge of the EGS_Mesh geometry (e.g., region volumes). We dynamic_cast the base geometry to a EGS_Mesh pointer to access its data and prepare to save a Gmsh file. Also, copy the input file to a new output file and output result data vectors. We're writing out scoring data to Gmsh for quick inspection after a simulation. Increase the maximum number of media to 200 for ICRP 145: the MRCP_AF and AM models have 52 media each, and ICRP 145 meshes have ~190 organs (media) and are both around 8M elements.

See issue #863

With this fix, the results from the two mesh libraries appear nearly identical. The from-scratch implementation seems 20-30% slower which may be due to not pre-calculating face normal vectors. This and other potential optimizations will be investigated. There has not been any profiling done on either implementation. Also, write out element doses in Gy and remove stray debugging print statements

Before (without inlining): 15.21% 15.21% mevegs libegs_mesh.so [.] EGS_Mesh::is_boundary After precalculation (without inlining): 4.64% 4.64% mevegs libegs_mesh.so [.] EGS_Mesh::is_boundary After switching to std::uint8_t to avoid bit_reference costs (wi): 3.89% 3.89% mevegs libegs_mesh.so [.] EGS_Mesh::is_boundary ...and moving the function to the header allowed it to be inlined.

The 1% overall gain for vector<uint8_t> was for a pathological case where the source was pointed away from the mesh. Space savings may become more important as the mesh grows larger.

Implement mesh reordering to try and reduce time spent in the isWhere loop. Simulation results are unchanged (as they should be). A first measurement is promising: a reduction from 75% to 5% of simulation time spent in isWhere. There is a noticeable difference when running simulations too: time spent in one batch went from ~13s to 4s. The EGS_Vector used as the reference point is hard-coded for the current test. Later commits will instead use the position of a particle generated by the source.

Reordering was commented out for a clean comparison. The isWhere cost was reduced from ~75% without reordering to ~3.5%.

This is a combination of 6 original commits related to small improvements in egs_mesh and the mevegs application: - don't print results for large meshes (>100) - improve egs_mesh initialization warnings - silence int <-> size_t comparison warning - check that the output file was actually opened - make output file name based on input and msh name - set mesh labels and boundary tolerance from input

The implementation is somewhat hacky but is intended for the most common case: wrapping an EGS_Mesh in a single simple, base geometry (e.g. some CAD model in an air box). This removes the requirement for users to mesh a bounding box at the Gmsh stage.

EGS_Mesh wasn't displaying properly in egs_view, because the normals weren't being set correctly.

- count overlapping elements in both bounding boxes - better octant numbering for reduced lookup cost - helper functions galore!

Rarely (~1.5% for a 600,000 element mesh), part of an element may be contained within an octant but not counted as a bounded element because none of the nodes or centroid are within the octant. The previous implementation wrongly considered these points as outside (-1), while the brute force method disagreed and found the correct element. This issue was overcome by not immediately halting the search, but by storing the size of the largest tetrahedral edge and comparing it to the bounding box of the octant searched so far. If the search fails in the most likely octant, its siblings are searched and then the parent's siblings, etc, expanding the search area by moving up the octree. Only when this "largest possible tetrahedron search" is completed, does the octree conclude the point is indeed not inside the mesh.

Refactor initialization code to carefully initialize static members.

See http://realtimecollisiondetection.net/books/rtcd/errata/ "Errata 4/27/05"

We also use a better tetrahedron-box intersection check to avoid requiring as many fallback searches for tetrahedra that are barely overlapping an octant. There are likely still some latent errors in the implementation because of the fragility of the floating-point routines.

Compared to both howfar and isWhere, hownear is much more expensive even in the best case. We have to find the minimum distance to the mesh in any direction, not just along one direction like howfar. So even when a candidate element is found, we have to continue searching until we conclude that no closer elements exist by testing against octree bounding boxes. If the minimum distance to an octant is larger than the current best distance, we can skip searching that octant. Using this, we can avoid checking much of the total surface mesh, especially if we find a good candidate element early on. Even so, many calls to hownear will harm the overall simulation efficiency. This is currently the case for envelope simulations, since hownear is called many times in the envelope before it reaches the mesh. We're currently doing a full check for each hownear call, but could potentially do a global bounding box check and return the distance to the outer bounding box if the point is outside the mesh. Hopefully this fudged distance won't impact correctness.

This change improves the overall simulation efficiency for envelope electron simulations by avoiding an exhaustive check to find the real minimum distance to the mesh. Instead, the distance to one of the surface tree's bounding boxes may returned. For example, this caused the efficiency of the MRCP_AM model in an air box with 1MeV photons, to increase from ~400 to ~5000 histories per second. This is allowed by the HOWNEAR spec (section 3.6 PIRS-701) which states: "In complex geometries, the mathematics of HOWNEAR can become difficult and sometimes almost impossible! If it is easier for the user to compute some lower bound to the nearest distance, this could be used..." Another alternative given in that document is setting the hownear result to 0.0, forcing ELECTR to assume it is very close to a boundary at all times, but that was found to be less efficient than this change. So, given an input point to hownear, there are two possibilities: 1. If the point is outside the global bounding box, return the distance to the closest point on the global bounding box and skip descending the tree. 2. If the point is inside the global bounding box, descend the tree until you reach a leaf node that contains the query point. A lower bound can then be obtained by returning the minimum of the distance to the bounding box sides and any boundary elements contained in the box.

Before this change, we were doing a full closest-triangle-point calculation for query points inside the tetrahedron, which is overkill because if a point is inside a tetrahedron, the shortest point will be along a face normal. This translates to an approximate 10% performance improvement for electron simulations. This commit changes the minimum-interior calculation to a dot product and 3D vector subtraction (i.e. distance to the triangle plane). We are precalculating the face normals to further improve simulation performance (5% faster compared to the baseline, calculating the normals on the fly), but this comes with a fairly large increase in memory. Each tetrahedron face normal has 3 double values, so the additional memory per element is (4 x 3 x 8 =) 96 bytes. For the 8-million element MRCP models, that is around 800 MB. We will offset this increase by no longer storing the node coordinate data inline, but instead use offsets into a node data vector. This precalculation step will be somewhat of a waste for photon simulations and if memory use becomes a problem it might be worthwhile making this into a optional flag to be configured in the input file.

This change avoids copying node coordinate data over and over in the octree constructor and instead looks it up using the EGS_Mesh. This reduces the memory spike observed during initialization, which was as high as 5GB but is now lowered to 4GB for the MRCP_AM model. The steady state is around 2GB but should be lower after node coordinate information is de-duplicated.

This change reduces the total memory used by only storing node coordinate information once instead of each time a node is used in a tetrahedron. For the 8-million tetrahedron MRCP_AM model, peak memory dropped from 4GB to 3.4GB and steady-state memory usage dropped from 2GB to 1.4GB, according to the massif heap profiler. This does not appear to have caused simulation runtime performance to worsen at all.

Besides being more straightforward at a high-level, this change is intended to reduce peak memory usage as much as possible in the constructor. Before this change, there were a number of large temporaries (e.g. the input vector of elements) that were still alive during memory-intensive operations like building the octrees and constructing the face-normal vector. The important changes here are: 1. The potentially large input vectors to the constructor are dropped at the end of initializeElements as soon as the EGS_Mesh equivalents are built. 2. The octree initialization step was moved before the face normal initialization step, because the octree step is fairly memory intensive and the face normal data is potentially large (~0.8GB for an 8 million element mesh). The face normals aren't used by the octree setup so leave calculating them to the end of the constructor after all the extra memory used to construct the octrees is freed up. For the MRCP_AM model (8.2 million elements) these changes brought the peak memory usage down from 3.4GB to 2.2GB, which is around 1.5 times the steady state usage.

- Use std::array on stack instead of heap vector in Face constructor. - Since we know both node and element numbers fit into an int, use int instead of std::size_t in neighbour-finding algorithms. This change lowered the memory used during the neighbour-finding step down from 2.2GB to ~1.6 GB and shifted the location of the global memory peak. The new global memory peak is during the face normal calculation step (~1.9GB).

This avoids any potential copies of these large vectors.

This is a combination of 4 original commits devoted to removing dead code in the egs_mesh library: 1. Remove dead code. 2. Remove unnecessary isWhere logic. With the introduction of the better box-tetrahedron intersection check, we know, barring any intersection check bugs, that if a point isn't found in an octant, we can immediately exit instead of searching until a minimum bounding box around the largest edge is reached. This simplifies both the isWhere implementation and lets us get rid of the octree's static members. 3. Remove obsolete kdtree classes. 4. Remove obselete reordering code: While reordering did improve performance over a brute search, performance got worse and worse farther away from the source as expected. The octree methods are faster and have better scaling.

Before this change, I tried to compile a unit test suite for EGS_Mesh, that included "egs_mesh.h". There was a confusing compiler error: error: invalid application of ‘sizeof’ to incomplete type ‘EGS_Mesh_Octree’ which stems from defining the destructor in the header file, where it can't properly destruct the forward declared EGS_Mesh_Octree class yet. For some reason, this was the first time this error came up even though we include "egs_mesh.h" in the mevegs user code for example. Also, use unique geometry name in EGS_Mesh constructor.

Fix a major bug in howfar_interior. The previous implementation did not ensure that the particle had to be "in front" of the tetrahedron face it was intersecting with. So if the particle was slightly outside of the tetrahedron as in the case of numerical undershoot, the wrong face would be intersected. Here is a diagram of the situation: /\ * -> / \ /____\ The particle is numerically outside the tetrahedron, but EGS considers it inside and we are required to calculate howfar as if it were inside (see Bielajew's "HOWFAR and HOWNEAR: Geometry Modelling for Monte Carlo Particle Transport"). So the near side must not be intersected. But the previous implementation did not account for this and so sometimes a particle travelling along one direction would be transported into a neighbour element in the opposite direction. We now ensure the particle must be on the "right" side of the tetrahedron face to intersect it. In all there are three cases now handled in howfar_interior: 1. The position is inside the region: calculate the intersection without further complications. 2. The position is outside the region but will intersect one of the i interior faces: calculate the intersection, ignoring any backwards facing faces. 3. The position is outside the region and won't intersect any of the interior faces: return 0.0 as the howfar distance and set the new region number to match the position (even if it's not an immediate neighbour of the region in question). Case 1 | Case 2 | Case 3 | | /\ | /\ | /\ / \ | / \ | / \ / \ | / \ | / \ / * -> X | * -> / X | <- * / \ /________\ | /________\ | /________\ | | Intersection | Intersection | No intersection Cases 1 and 2 are both handled by the same logic for now, but it might be more efficient to use different equations for both (e.g. ray-plane intersection algorithms for case 1).

Add a dedicated "parser-result" class (EGS_MeshSpec) has a few benefits: - It allows us to differentiate parsing errors from EGS_Mesh constructor errors. - It makes writing new parsers more straightforward by defining EGS_MeshSpec as the interface between the parser and EGS_Mesh. - It makes the EGS_Mesh constructor cleaner.

By avoiding buffering output data into result vectors, we avoid an increase in heap memory previously observed while writing the Gmsh output file.

Output results to a VTK file. Avoid writing a useless intermediate VTK file parallel job and instead only write it out for serial jobs (including calculate=combine runs). Also, make EGS_Mesh file output methods const correct.

Instead of crashing when a particle is stuck at an interface between regions (an infinite loop is detected), push particles along their momentum vector to ensure the simulation can proceed.

Before this change, simple two-element infinite loops were detected, but not those involving more elements. This change tracks small steps and if a potential stuck particle is detected, pushes the particle along the momentum vector to continue the simulation. This case was detected with an internal point source simulation where the point source was exactly on a node shared by many tetrahedra.

Move to thick planes to make howfar_interior more robust. The previous version used ray-triangle intersections, which proved in practice to be fragile (infinite loops, missed intersections at corners, etc.). It is probably fine to use a modified ray-triangle intersection test, but using thick planes is easier to understand and simpler to implement. I believe this is also what is used in Geant, PHITS, and MCNP. I'm also using a small step epsilon (EGS_Mesh::min_step_size) to ensure a minimum step is taken. This could hang if this epsilon is near the rounding error (i.e. for very large mesh cells). But if the epsilon is too large, then very thin tetrahedra might be skipped over. If it's not feasible to find a minimum step automatically, it should at least be made user-adjustable (i.e. using EGS_BaseGeometry::boundaryTolerance). Don't check the smallest absolute distance to the thick plane, which can consider points far away from a tetrahedron (since the point could be closer to another plane). Instead, check the largest negative signed distance isn't larger than the tolerance. That way, points very far away from the tetrahedron but close to one of its infinite face planes aren't considered part of the thick plane. First testing of this new version was very encouraging: good performance and no hangs for known problem cases.

Use addMedium to get the correct offset for media initialization, instead of starting from 0. If addMedium is called with a name already in the media list, the corresponding offset is returned. Otherwise, the media is added to the list and the new offset is returned. Before this change, EGS_Mesh used to always define element media offsets starting at 0. So if any media were defined before the mesh was created, the mesh media offsets would be wrong.

For large meshes over 500,000 elements, the EGS_Mesh setup time can take a few minutes. Add percent indicators for the setup steps that take the most time: - file parsing - neighbour finding - octree building The indicators let the user know everything is working during construction and that longer pauses are normal. Only log EGS_Mesh cstor progress for large meshes (smaller ones are built almost instantly so require no logging). The indicators are drawn on a single tty line using carriage returns. If the program is running as a batch job, only the final summary line is printed to avoid cluttering the log file. This feature has been tested on Windows and Linux. Example (interactive, all lines are redrawn on the same line): EGS_Mesh: reading 1030301 nodes (0%) # start EGS_Mesh: reading 1030301 nodes (1%) # in progress ... ... EGS_Mesh: reading 1030301 nodes (100%) EGS_Mesh: read 1030301 nodes in 1.541s # finish Example (batch, only the final summary is printed): EGS_Mesh: read 1030301 nodes in 1.541s # finish

This is a combination of 6 original commits related to fixing some small issues in the egs_mesh library: - fix EGS_Mesh setup error formatting - add performance note to EGS_Mesh_Octree routine - fix move constructor and assignment definitions (can't be in header) - add a proper EGS_MeshSpec constructor - fix parse_msh_file interface definition - make EGS_Mesh constructor explicit - pass element index to EGS_Mesh neighbour function - fix EGS_Mesh and mevegs makefile dependencies

This is a combination of 12 original commits related to cleaning up the egs_mesh library source code: - only throw std::runtime_error for uniform behaviour - improve missing mesh file error message - fix and clean up comments - remove unused class members - silence mesh_neighbours sign warning - fix library warnings for -Wall - fix inside definition to return region index, for correctness - fix bounding box print function - remove unused parameters from howfar_exterior - remove overly defensive range check in hownear - remove noisy warning in howfar recovery, now obsolete - set base geometry error flag instead of crashing for rare errors

Exclude implementation details such as EGS_Mesh_Octree, mesh_neighbours and msh_parser.

Add a scale key so mesh files that aren't in cm can be scaled appropriately.

This commit adds a TetGen node and ele file parser to EGS_Mesh. This lets users provide meshes in the TetGen format, which is used by the ICRP 145 mesh phantoms, among others. Previously a conversion step was required to convert the ICRP 145 meshes to msh4.1. Meshes from TetGen can now be directly simulated using EGS_Mesh. The TetGen format uses a .node file to store node information and .ele file for element information. Either the node or ele file can be provided for the EGS_Mesh `file` input, the other will be opened automatically. Apart from adding direct support for important meshes from the literature, the TetGen format is simpler than the Gmsh format, offering an alternative for users generating their own mesh files. The Geant4 and PHITS tetrahedral mesh implementations also have support for TetGen files, so the same mesh can be used with each code.

Add the sed 's' command which was missing when setting the system name. Also, change the sed separator from ':' to '|' to avoid ambiguity with sed labels.

Add the density correction files that were used in the ICRM high purity germanium Monte Carlo comparison. Also update the materials in the example input file to refer to the correct file names.

Set the c++ standard to c++14 if the g++ version is 5 or greater, otherwise fall back on c++11. Some systems still rely on version 4.8 which only supports up to c++11 features.

Remove the git branch information and provide the commit hash only. The branch name is just a label that can point to any commit, so it remains equivocal. This also solves the issue with long branch names where the fortran code would fail to compile because of exceeding the line length. Also, fix the comparisons between egs_home and HEN_HOUSE/user_code apps, so that it doesn't report a difference when the HEN_HOUSE version doesn't exist.

Change the egs_radionuclude_source to take an input "base source" instead of defining a collimated or isotropic source. This way, any of the egs++ sources or a custom source can now be used to define the distribution of the radionuclide. Also, add warning messages for deprecated radionuclide source inputs. This change was originally motivated by simulations of targeted radionuclide therapy.

Add tests in the configure scripts to avoid setting the gcc memory model flag -mdmodel=medium if a darwin (macos) system is detected. This flag causes compilation failures at the linking stage.

Change the binary format to the one suggested by Randle Taylor, and add information about the sparse matrix binary dose format in the DOSXYZnrc documentation.

Add an option to select the dose output format in the gui. Just emulating what was done for the howfarless option.

Change the egs_base_geometry getMass virtual function to a getVolume function that returns region volume instead of mass. Typically, initGeometry() is invoked before any media data is read in the application, hence getMass was unreliable. The getMass function was replaced with a getVolume function, which only depends on internal gemoetry definitions. The geometries in egs_autoenvelope, egs_nd_geometry, egs_rz, and egs_spheres are affected. The ausgab dose scoring object which outputs a 3ddose file was also changed to reflect the change from getMass to getVolume. This modification of egs_dose_scoring also now properly accounts for relativeRho in the mass calculation. Also, handle vacuum medium index appropriately in getMediumRho.

Fix compilation failure of BEAMnrc accelerators since the addition of git and compile time logging. This also adds back printing of the git hash, which was mistakenly removed in PR #898. Fix this for the compilation of the shared library as well.

The g++ compiler executable may have a longer name, for different installed versions (e.g., g++-11). Also, fix the wrong variable name printed in the option list.

Revert to the egs_brachy encoding for egsphant files, allowing for 62 media encoded with the set of alphanumeric characters (where 0 is reserved for vacuum). Briefly, #633 expanded the number of media in egsphant files to 95, using all printable ascii characters. However, the number of media was also independently expanded in egs_glib for the development of egs_brachy, but using a different encoding consisting of only the 62 alphanumeric ascii characters. This led to inconsistent egsphant files that are no longer interchangeable. This issue was originally reported and discussed in clrp-code/egs_brachy#22.

Fix the egs_chamber examples so that they now run without needing to compile any BEAMnrc accelerators or make changes to the input file. Also fix a typo in the example2 water box size.

- EGS_Interpolator: increase array size when creating an interpolator based on an array of value pairs to have an extra bin to prevent round-off errors. This is also done by PEGS4 when determining fit parameters for each interval. The extra bin uses the same parameters as the previous one, i.e., it extrapolates. - Add auxiliary methods to application classes to get electron step size, energy deposited and app source type. - Fluence scoring: at a plane, in regions for all particles.

Get proper stopping power interpolator depending on charge.

- Add "clean" target for dose and fluence scoring objects (instead of manually deleting files when building the ausgab object). Need to add to the rest of the ausgab objects! - One can now define DEBUG=-DDEBUG to turn on some debugging messages that I found useful. - Add MYDEF for user defines: use MYDEF=-DDEBUG to activate debugging messages and variables.

Reworked the classes for planar and volumetric fluence scoring to descend from a common base class which contains many commonalities, to avoid code duplication. Should have started that way!

Add auxiliary functions to application classes, required by the fluence scoring ausgab object.

- Add a method to get the charge of a base source, which returns -99 by default, matching the unknown type of the fluence scoring ausgab object. A base source does not know the charge of its particles unless it is later derived from a simple source. For instance a collection of sources, a beam source, or a phase space source will contain multiple type of particles. - Sources derived from the class EGS_SimpleSource return the source particle's charge. - Add a method to get source charge from the application.

Add an application-level method to obtain the maximum energy of the source. This proves useful for ausgab objects that need this information.

- Add option to define regions from where to score planar fluence to improve efficiency of charged particle calculations. - Move many common tasks to the base class EGS_FluenceScoring. - Fix group scoring region definitions. - Improve logic for defining scoring regions in planar scoring class. - Document classes with Doxygen comments.

- Add the ability to divide fluence in total and primary fluence. - Add the ability to calculate integral and differential fluence. - Improve output format. - Investigate the effect of scoring integral charged particle fluence using either EDEP (within differential scoring block) or directly TVSTEP: a) Identical for FLURZ-like calculation b) Small, negligible differences using stprO3 and stpwrO5 This allows extraction of the integral fluence scoring from the blocks for differential fluence. When one is only interested in the integral fluence, the calculation is up to 10% faster using a TVSTEP-based approach. Current implementation allows switching between EDEP-based and TVSTEP-based integral fluence calculation via the pre-processor define USETVSTEP. To use TVSTEP-based integral dose calculation use: make MYDEFS=-DUSETVSTEP Next iteration will only use the TVSTEP-based approach for integral fluence as both approaches result in negligible differences. - Add mechanism to track the distribution of number of energy bins covered when spreading EDEP. This is activated via the define variable DEBUG: to use it, build the fluence scoring ausgab object using: make MYDEFS=-DDEBUG One can also trace the difference between the step derived from spreading EDEP over the energy grid and TVSTEP via DEBUG.

- Add input blocks to define the energy grid and the scoring options based on the scoring type. - Improve output format.

Remove example for base class and add specific blocks for planar an volumetric scoring inputs.

- Remove initial approach to score charged particle integral fluence from each contribution to the differential fluence. Now is uses TVSTEP and computes integral fluence in one go. - Remove checks for pre-processor directive USETVSTEP as it not needed anymore. - Polish documentation.

This new ausgab object scores fluence in user-defined regions (volumetric) or at a circular or rectangular field (planar) located in an arbitrary position in space. Scoring can be done for either photons, electrons or positrons. Scoring for multiple particle types is accomplished by defining multiple fluence scoring ausgab objects.

Add a missing include of the array c++ class library. This avoids compiler errors that occurred on macOS.

Add ausgab calls before and after a photonuclear interaction, in order to handle cross section enhancement (xcse) correctly. Without these calls, there is no compensation for the enhanced cross section, hence the photonuclear cross section is effectively increased by the xcse factor. Also add the photonuclear event trapping in Fano calculations, for consistency. However, the Fano calculation fails when photonuclear attenuation is turned on anyway, because the photon is discarded without depositing its energy. Also update a couple loops that print out the ausgab flags (IAUSFL) values, to include the photonuclear attenuation option value.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Merge branch develop into master for 2023 release #965

Merge branch develop into master for 2023 release #965

Commits on Jun 27, 2022

Commits on Jun 28, 2022

Commits on Jun 30, 2022

Commits on Jul 4, 2022

Commits on Jul 8, 2022

Commits on Jul 12, 2022

Commits on Jul 15, 2022

Commits on Jul 16, 2022

Commits on Jul 18, 2022

Commits on Jul 22, 2022

Commits on Jul 25, 2022

Commits on Jul 27, 2022

Commits on Oct 18, 2022

Commits on Oct 31, 2022

Commits on Jan 10, 2023

Commits on Jan 16, 2023

Commits on Jan 23, 2023

Commits on Feb 15, 2023