From 240b10fa122846ae6862e3101f209a20085c2064 Mon Sep 17 00:00:00 2001 From: Reid Townson Date: Wed, 21 Sep 2016 14:22:03 -0400 Subject: [PATCH] Add egs++ volume calculation Add a Monte Carlo volume calculation to egs++. It can be turned on using :start volume calculation: in any input file. The recursive algorithm builds an octree with higher resolution in geometries with more regions. A user specified box is used as the boundary for sampling. --- HEN_HOUSE/doc/src/pirs898-egs++/geometry.doxy | 45 ++++ HEN_HOUSE/egs++/egs_application.cpp | 20 +- HEN_HOUSE/egs++/egs_base_geometry.cpp | 188 ++++++++++++++++ HEN_HOUSE/egs++/egs_base_geometry.h | 200 ++++++++++++++++++ 4 files changed, 447 insertions(+), 6 deletions(-) diff --git a/HEN_HOUSE/doc/src/pirs898-egs++/geometry.doxy b/HEN_HOUSE/doc/src/pirs898-egs++/geometry.doxy index 977b27f8f..8b55a4b9e 100644 --- a/HEN_HOUSE/doc/src/pirs898-egs++/geometry.doxy +++ b/HEN_HOUSE/doc/src/pirs898-egs++/geometry.doxy @@ -317,6 +317,51 @@ developing new geometries and has helped to find numerous bugs in the initial implementations of the various geometry classes provided with \c egspp. +\section volume_calculation Calculating volumes + +The volumes of regions or media in your geometry can be calculated via +Monte Carlo sampling using the built in volume calculation. To perform the +volume calculation, you define a box within which test points will be generated. +Make sure that the box encloses your entire volume(s) of interest. + +The volumes of each region and medium are automatically printed as output. +The user may also specify regions to calculate the volume in using labels. +After setting labels in geometries of interest, list the labels under +the \c labels input. For each label, the sum of the region volumes is output. + +To turn on volume calculation, include a new control block in your input +file: +\verbatim +:start volume calculation: + passes = the number of recursion passes to create the octree + samples = the number samples to perform per node per pass + box min = 0 0 0 + box max = 100 100 100 + labels = labels you have added to your geometry using `set label` +:stop volume calculation: +\endverbatim + +For example: +\verbatim +:start volume calculation: + passes = 20 + samples = 1e7 + box min = 0 0 0 + box max = 100 100 100 + labels = detector +:stop volume calculation: +\endverbatim + +For the volume calculation, the sampling box is divided into progressively +smaller voxels as the geometry is analyzed. +Voxels are created in an \em octree, to provide higher sampling resolution in +areas with a greater density of geometry regions. + +As the volume calculation progresses, the number of samplings in each region is +adjusted according to the uncertainty on its volume. On each pass, regions +with a higher uncertainty are sampled more often, while regions with lower +uncertainty are sampled less often. + \section geometry_view The geometry viewer: egs_view The geometry package comes with a viewer that can be used diff --git a/HEN_HOUSE/egs++/egs_application.cpp b/HEN_HOUSE/egs++/egs_application.cpp index 92e3a9a3a..6f96caed3 100644 --- a/HEN_HOUSE/egs++/egs_application.cpp +++ b/HEN_HOUSE/egs++/egs_application.cpp @@ -690,6 +690,11 @@ int EGS_Application::initGeometry() { } geometry->ref(); geometry->setApplication(this); + + if(rndm) { + geometry->printVolumes(input,rndm); + } + return 0; } @@ -749,9 +754,17 @@ int EGS_Application::initRNG() { int EGS_Application::initSimulation() { //if( !input ) { egsWarning("%s no input\n",__egs_app_msg2); return -1; } - egsInformation("In EGS_Application::initSimulation()\n"); + + egsInformation("\n================================================================================\n"); + egsInformation("EGS_Application::initSimulation()\n"); + egsInformation("Initializing simulation...\n\n"); int err; bool ok = true; + err = initRNG(); + if (err) { + egsWarning("\n\n%s RNG initialization failed\n",__egs_app_msg2); + ok = false; + } err = initGeometry(); if (err) { egsWarning("\n\n%s geometry initialization failed\n",__egs_app_msg2); @@ -762,11 +775,6 @@ int EGS_Application::initSimulation() { egsWarning("\n\n%s source initialization failed\n",__egs_app_msg2); ok = false; } - err = initRNG(); - if (err) { - egsWarning("\n\n%s RNG initialization failed\n",__egs_app_msg2); - ok = false; - } err = initRunControl(); if (err) { egsWarning("\n\n%s run control initialization failed\n",__egs_app_msg2); diff --git a/HEN_HOUSE/egs++/egs_base_geometry.cpp b/HEN_HOUSE/egs++/egs_base_geometry.cpp index e622679d1..ee9ef2ca2 100644 --- a/HEN_HOUSE/egs++/egs_base_geometry.cpp +++ b/HEN_HOUSE/egs++/egs_base_geometry.cpp @@ -46,6 +46,8 @@ #include "egs_library.h" #include "egs_input.h" #include "egs_application.h" +#include "egs_rndm.h" +#include "egs_timer.h" #include #include @@ -1170,3 +1172,189 @@ int EGS_BaseGeometry::setLabels(const string &inp) { return 1; } + +void EGS_BaseGeometry::printVolumes(EGS_Input *input, EGS_RandomGenerator *rndm) { + + EGS_Timer timer; + + EGS_Input *inp; + inp = input->takeInputItem("volume calculation"); + if (!inp) { + return; + } + + egsInformation("\n"); + egsInformation("================================================================================\n"); + egsInformation("EGS_BaseGeometry::printVolumes()\n"); + egsInformation("Performing a Monte Carlo volume calculation...\n"); + egsInformation("================================================================================\n"); + + double npass; + if (inp->getInput("passes", npass)) { + egsFatal("ERROR: passes = undefined.\nAborting."); + } + double nsample; + if (inp->getInput("samples", nsample)) { + egsFatal("ERROR: passes = undefined.\nAborting."); + } + std::vector boxmin, boxmax; + if (inp->getInput("box min", boxmin)) { + egsFatal("ERROR: box min = undefined.\nAborting."); + } + if (inp->getInput("box max", boxmax)) { + egsFatal("ERROR: box max = undefined.\nAborting."); + } + + EGS_Vector min(boxmin[0], boxmin[1], boxmin[2]); + EGS_Vector max(boxmax[0], boxmax[1], boxmax[2]); + + std::vector volumeLabels; + inp->getInput("labels", volumeLabels); + + // get number of regions + int nreg = regions(); + + egsInformation("Number of passes = %.0f\n", npass); + egsInformation("Number of samples = %.0f\n", nsample); + egsInformation("Box min = %f %f %f\n", boxmin[0], boxmin[1], boxmin[2]); + egsInformation("Box max = %f %f %f\n", boxmax[0], boxmax[1], boxmax[2]); + egsInformation("Number of regions = %d\n", nreg); + + // create octree root node + Volume_Info *vInfo = new Volume_Info(min, max, this, rndm); + Volume_Node *root = new Volume_Node(vInfo, NULL, -1); + + // sample points inside the node (this is where all the work gets done) + int nsampleLeaf = nsample; + egsInformation("\n"); + egsInformation("================================================================================\n"); + egsInformation("%15s %15s %15s %15s\n", "PASS", "octree leaves", "samples/leaf", "total points"); + egsInformation("================================================================================\n"); + vector sigma, volume, nhit, ntry; + double avgSigma = 0, ntot = 0; + for (int i=0; icountLeaves(), nsampleLeaf, root->countLeaves()*nsampleLeaf); + + if(i > 0) { + double varSum = 0; + for (int j=0; j 0) { + sigma[j] = sqrt(sigma[j])/volume[j]*100.; + } else { + sigma[j] = 0; + } + varSum += sigma[j]; + } + avgSigma = varSum / nreg; + } + + root->sampleVolumes(vInfo, nsampleLeaf, sigma, avgSigma); + + ntot = 0; + for (int j=0; jcalculateVolumes(vInfo, volume, sigma, nreg, ntot, nhit, ntry); + + nsampleLeaf = nsample/root->countLeaves() + 1; + if (nsampleLeaf < 10) nsampleLeaf = 10; + } + + // report regions + egsInformation("\n"); + egsInformation("================================================================================\n"); + egsInformation("%15s %15s %15s %15s %%\n", "REGION", "volume", "sigma", "sigma"); + egsInformation("================================================================================\n"); + double total_volume = 0; + double total_sigma = 0; + for (int i=0; iisRealRegion(i)) continue; + double vol = volume[i]; + double sig = sigma[i]; + double err = sqrt(sig); + double relerr = 0; + if (vol > 0) relerr = err/vol*100; + egsInformation("%15d %15.6f %15.6f %15.6f %%\n", i, vol, err, relerr); + total_volume += vol; + total_sigma += sig; + } + double total_err = sqrt(total_sigma); + double relerr = 0; + if (total_volume > 0) relerr = total_err/total_volume*100; + egsInformation("\n"); + egsInformation("%15s %15.6f %15.6f %15.6f %%\n", "TOTAL", total_volume, total_err, relerr); + + // report media + egsInformation("\n"); + egsInformation("================================================================================\n"); + egsInformation("%15s %15s %15s %15s %%\n", "MEDIUM", "volume", "sigma", "sigma"); + egsInformation("================================================================================\n"); + int nmed = this->nMedia(); + vector vol_med, sig_med; + for (int i=0; iisRealRegion(i)) continue; + int med = this->medium(i); + if (med >=0 && med < nmed) { + vol_med[med] += volume[i]; + sig_med[med] += sigma[i]; + } + } + for (int i=0; i 0) relerr = err/vol*100; + egsInformation("%15s %15.6f %15.6f %15.6f %%\n", this->getMediumName(i), vol, err, relerr); + } + egsInformation("\n"); + + // report labels + if (volumeLabels.size() > 0) { + egsInformation("================================================================================\n"); + egsInformation("%15s %15s %15s %15s %%\n", "LABEL", "volume", "sigma", "sigma"); + egsInformation("================================================================================\n"); + for (int k=0; k label_regions; + this->getLabelRegions(volumeLabels[k], label_regions); + double label_volume = 0; + double label_sigma = 0; + for (int i=0; i 0) relerr = err/label_volume*100; + egsInformation("%15s %15.6f %15.6f %15.6f %%\n", volumeLabels[k].c_str(), label_volume, err, relerr); + } + egsInformation("\n"); + } + + EGS_Float cpu = timer.time(); + + // report total number of samples + egsInformation("Total number of points sampled = %.0f\n", ntot); +// egsInformation("Memory used = %.2f MB\n", root->countMem()); + egsInformation("Calculation time = %g seconds\n\n", cpu); + + egsInformation("EGS_BaseGeometry::printVolumes(): Volume calculation complete.\n\n"); + + root->deleteChildren(); + delete root, vInfo; +} diff --git a/HEN_HOUSE/egs++/egs_base_geometry.h b/HEN_HOUSE/egs++/egs_base_geometry.h index 0305b1295..c6915dbec 100644 --- a/HEN_HOUSE/egs++/egs_base_geometry.h +++ b/HEN_HOUSE/egs++/egs_base_geometry.h @@ -43,6 +43,7 @@ #define EGS_BASE_GEOMETRY_ #include "egs_vector.h" +#include "egs_rndm.h" #include #include @@ -729,6 +730,9 @@ class EGS_EXPORT EGS_BaseGeometry { /*! \brief Set the labels from an input string */ int setLabels(const string &inp); + + /*! \brief Calculate the volumes of regions in the geometry */ + void printVolumes(EGS_Input *input, EGS_RandomGenerator *rndm); protected: @@ -1335,5 +1339,201 @@ struct EGS_GeometryIntersections { and the complete implementation: */ +class Volume_Info { +public: + EGS_Vector min, max; + EGS_BaseGeometry *g; + EGS_RandomGenerator *rng; + + // constructor + Volume_Info(EGS_Vector &bbmin, EGS_Vector &bbmax, EGS_BaseGeometry *geom, EGS_RandomGenerator *generator) { + min = bbmin; + max = bbmax; + g = geom; + rng = generator; + } +}; + +class Volume_Node { +public: + double ntry; + vector nhit; + int ix, iy, iz; ///< octree indices, representing the binary location code of the node in x, y, z; + int region; ///< region number of the node + unsigned short level; ///< depth of the node (root node is level 0) + Volume_Node *parent; ///< pointer to the parent node (only root node can have parent set to NULL) + Volume_Node **child; ///< pointer to children nodes, NULL is there are no children + + // constructor + Volume_Node(Volume_Info *v, Volume_Node *p, int index) { + + // init child + child = NULL; + + // set parent + parent = NULL; + if (p != NULL) parent = p; + + // set level + if (parent != NULL) { + level = parent->level + 1; + ix = (parent->ix << 1) | (index>>0 & 0x1); + iy = (parent->iy << 1) | (index>>1 & 0x1); + iz = (parent->iz << 1) | (index>>2 & 0x1); + } + else { + level = 0; + ix = iy = iz = 0; + } + + // set region number of the node + setRegion(v); + + // allocate array of region hits + int nreg = v->g->regions(); + for (int i=0; imax - v->min) * (1.0/(double)(1<min; + r0.x += (ix+0.5)*d.x; + r0.y += (iy+0.5)*d.y; + r0.z += (iz+0.5)*d.z; + region = v->g->isWhere(r0); + } + + // create children nodes + void createChildren(Volume_Info *v) { // create children to the this node + if (!child) { // ensure there are no children already + child = new Volume_Node* [8]; // allocate memory for 8 new children nodes + for (int i=0; i<8; i++) { + child[i] = new Volume_Node(v, this, i); + } + } + } + + // delete children (recursively) + void deleteChildren() { // delete children of this node + if (child) { // if this node has children, delete them + for (int i=0; i<8; i++) { + child[i]->deleteChildren(); // recursive calls to delete all children branches + delete child[i]; + } + delete [] child; // free the memory allocated for the 8 children + child = NULL; // set children pointer to NULL to indicate that there are no children anymore + } + } + + // count the number of octree leaves + int countLeaves() { + int count = 1; + if (child) { + count -= 1; + for (int i=0; i<8; i++) { + count += child[i]->countLeaves(); + } + } + return count; + } + + // count memory + double countMem() { + int count = 1; + if (child) { + for (int i=0; i<8; i++) { + count += child[i]->countLeaves(); + } + } + return (double)count*sizeof(Volume_Node)/1024.0/1024.0; + } + + // sample region volumes + void sampleVolumes(Volume_Info *v, int nsample, vector sigma, double avgSigma) { + + if (child) { + for (int i=0; i<8; i++) { + child[i]->sampleVolumes(v, nsample, sigma, avgSigma); + } + } + else { + EGS_Vector d = (v->max - v->min) * (1.0/(double)(1<min.x + ix*d.x; + double miny = v->min.y + iy*d.y; + double minz = v->min.z + iz*d.z; + int nreg = v->g->regions(); + int nsampleAdjusted; + if(avgSigma > 0) { + // Adjust the number sampled by this leaf using the + // deviation of the uncertainty from the average + // This increases the number of samples in regions with + // higher uncertainty + // The efficiency improvement is a factor of about 3.8 + nsampleAdjusted = int(double(nsample) * sigma[region]*sigma[region] / (avgSigma*avgSigma)); + //nsampleAdjusted = int(double(nsample) * sigma[region] / avgSigma); + //nsampleAdjusted = nsample; + } else { + nsampleAdjusted = nsample; + } + + for (long long i=0; irng->getUniform(); + double ry = v->rng->getUniform(); + double rz = v->rng->getUniform(); + r.x = minx + (d.x*rx); + r.y = miny + (d.y*ry); + r.z = minz + (d.z*rz); + int ireg = v->g->isWhere(r); + if (level < 8 && ireg != region) { + createChildren(v); + break; + } + else if (ireg >= 0 && ireg < nreg) { + nhit[ireg]++; + } + ntry++; + } + } + } + + // calculate region volumes + void calculateVolumes(Volume_Info *v, vector &volume, vector &sigma, int nreg, double &ntot, vector ®ionHits, vector ®ionTrys) { + + if (child) { + for (int i=0; i<8; i++) { + child[i]->calculateVolumes(v, volume, sigma, nreg, ntot, regionHits, regionTrys); + } + } + else { + EGS_Vector d = (v->max - v->min) * (1.0/(double)(1<= 0 && nhit[region] == ntry) { + volume[region] += vol; + regionHits[region] += nhit[region]; + regionTrys[region] += ntry; + } + else if (ntry > 0) { + for (int i=0; i