diff --git a/HEN_HOUSE/egs++/egs_base_geometry.cpp b/HEN_HOUSE/egs++/egs_base_geometry.cpp index ee9ef2ca2..e037c5ba0 100644 --- a/HEN_HOUSE/egs++/egs_base_geometry.cpp +++ b/HEN_HOUSE/egs++/egs_base_geometry.cpp @@ -1189,12 +1189,12 @@ void EGS_BaseGeometry::printVolumes(EGS_Input *input, EGS_RandomGenerator *rndm) egsInformation("Performing a Monte Carlo volume calculation...\n"); egsInformation("================================================================================\n"); - double npass; - if (inp->getInput("passes", npass)) { + double npass_double; + if (inp->getInput("passes", npass_double)) { egsFatal("ERROR: passes = undefined.\nAborting."); } - double nsample; - if (inp->getInput("samples", nsample)) { + double nsample_double; + if (inp->getInput("samples", nsample_double)) { egsFatal("ERROR: passes = undefined.\nAborting."); } std::vector boxmin, boxmax; @@ -1204,70 +1204,33 @@ void EGS_BaseGeometry::printVolumes(EGS_Input *input, EGS_RandomGenerator *rndm) 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); + EGS_I64 npass = (EGS_I64)npass_double; + EGS_I64 nsample = (EGS_I64)nsample_double; + // get number of regions - int nreg = regions(); + unsigned int nreg = regions(); - egsInformation("Number of passes = %.0f\n", npass); - egsInformation("Number of samples = %.0f\n", nsample); + egsInformation("Number of passes = %lld\n", npass); + egsInformation("Number of samples = %lld\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); + egsInformation("Number of regions = %u\n", nreg); // 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; + vector sigma, volume; + vector nhit, ntry; + EGS_I64 ntot; + if (getVolumes(npass, nsample, boxmin, boxmax, nreg, sigma, volume, + nhit, ntry, ntot, rndm, true)) { + egsWarning("EGS_BaseGeometry::printVolumes(): Error: Failed to calculate volumes.\n"); + return; } // report regions @@ -1278,19 +1241,25 @@ void EGS_BaseGeometry::printVolumes(EGS_Input *input, EGS_RandomGenerator *rndm) double total_volume = 0; double total_sigma = 0; for (int i=0; iisRealRegion(i)) continue; + if (!this->isRealRegion(i)) { + continue; + } double vol = volume[i]; double sig = sigma[i]; double err = sqrt(sig); double relerr = 0; - if (vol > 0) relerr = err/vol*100; + 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; + 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); @@ -1306,7 +1275,9 @@ void EGS_BaseGeometry::printVolumes(EGS_Input *input, EGS_RandomGenerator *rndm) sig_med.push_back(0); } for (int i=0; iisRealRegion(i)) continue; + if (!this->isRealRegion(i)) { + continue; + } int med = this->medium(i); if (med >=0 && med < nmed) { vol_med[med] += volume[i]; @@ -1318,7 +1289,9 @@ void EGS_BaseGeometry::printVolumes(EGS_Input *input, EGS_RandomGenerator *rndm) double sig = sig_med[i]; double err = sqrt(sig); double relerr = 0; - if (vol > 0) relerr = err/vol*100; + if (vol > 0) { + relerr = err/vol*100; + } egsInformation("%15s %15.6f %15.6f %15.6f %%\n", this->getMediumName(i), vol, err, relerr); } egsInformation("\n"); @@ -1340,7 +1313,9 @@ void EGS_BaseGeometry::printVolumes(EGS_Input *input, EGS_RandomGenerator *rndm) } double relerr = 0; double err = sqrt(label_sigma); - if (label_volume > 0) relerr = err/label_volume*100; + if (label_volume > 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"); @@ -1354,6 +1329,85 @@ void EGS_BaseGeometry::printVolumes(EGS_Input *input, EGS_RandomGenerator *rndm) egsInformation("Calculation time = %g seconds\n\n", cpu); egsInformation("EGS_BaseGeometry::printVolumes(): Volume calculation complete.\n\n"); +} + +int EGS_BaseGeometry::getVolumes(EGS_I64 &npass, EGS_I64 &nsample, + vector &boxmin, + vector &boxmax, + unsigned int &nreg, vector &sigma, + vector &volume, + vector &nhit, vector &ntry, + EGS_I64 &ntot, + EGS_RandomGenerator *rndm, bool verbose = true) { + + // Check box dimensions + if (boxmin.size() != 3 || boxmax.size() != 3) { + egsWarning("EGS_BaseGeometry::getVolumes(): Error: Box dimensions for volume calculation 'min' and 'max' must have sizes of 3.\n"); + return 1; + } + + // Check vectors haven't already been initialized + if (sigma.size() > 0 || volume.size() > 0 || + nhit.size() > 0 || ntry.size() > 0) { + egsWarning("EGS_BaseGeometry::getVolumes(): Error: Expected empty vectors but they already have size (sigma, volume, nhit, ntry).\n"); + return 1; + } + + EGS_Vector min(boxmin[0], boxmin[1], boxmin[2]); + EGS_Vector max(boxmax[0], boxmax[1], boxmax[2]); + + // Create octree root node + Volume_Info *vInfo = new Volume_Info(min, max, this, rndm); + Volume_Node *root = new Volume_Node(vInfo, NULL, -1); + + // Initialize + double avgSigma = 0; + ntot = 0; + int nsampleLeaf = nsample; + for (unsigned int i=0; icountLeaves(), nsampleLeaf, root->countLeaves()*nsampleLeaf); + } + + if (i > 0) { + double varSum = 0; + for (unsigned int j=0; j 1e-10) { + 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 (unsigned int j=0; jcalculateVolumes(vInfo, volume, sigma, nreg, ntot, nhit, ntry); + + nsampleLeaf = nsample/root->countLeaves() + 1; + if (nsampleLeaf < 10) { + nsampleLeaf = 10; + } + } root->deleteChildren(); delete root, vInfo; diff --git a/HEN_HOUSE/egs++/egs_base_geometry.h b/HEN_HOUSE/egs++/egs_base_geometry.h index 5fe23cb1d..6fe4f9377 100644 --- a/HEN_HOUSE/egs++/egs_base_geometry.h +++ b/HEN_HOUSE/egs++/egs_base_geometry.h @@ -731,10 +731,17 @@ 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 */ + + /*! \brief Calculate and print the volumes of regions in the geometry */ void printVolumes(EGS_Input *input, EGS_RandomGenerator *rndm); + /*! \brief Calculate the volumes of regions in the geometry */ + int getVolumes(EGS_I64 &npass, EGS_I64 &nsample, vector &min, + vector &max, unsigned int &nreg, + vector &sigma, vector &volume, + vector &nhit, vector &ntry, EGS_I64 &ntot, + EGS_RandomGenerator *rndm, bool verbose); + protected: /*! \brief Number of local regions in this geometry @@ -1347,7 +1354,8 @@ class Volume_Info { EGS_RandomGenerator *rng; // constructor - Volume_Info(EGS_Vector &bbmin, EGS_Vector &bbmax, EGS_BaseGeometry *geom, EGS_RandomGenerator *generator) { + Volume_Info(EGS_Vector &bbmin, EGS_Vector &bbmax, EGS_BaseGeometry *geom, + EGS_RandomGenerator *generator) { min = bbmin; max = bbmax; g = geom; @@ -1357,13 +1365,13 @@ class Volume_Info { 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 + EGS_I64 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 if there are no children // constructor Volume_Node(Volume_Info *v, Volume_Node *p, int index) { @@ -1373,7 +1381,9 @@ class Volume_Node { // set parent parent = NULL; - if (p != NULL) parent = p; + if (p != NULL) { + parent = p; + } // set level if (parent != NULL) { @@ -1397,13 +1407,13 @@ class Volume_Node { } ntry = 0; } - + ~Volume_Node() { deleteChildren(); } // setRegion - void setRegion (Volume_Info *v) { + void setRegion(Volume_Info *v) { EGS_Vector d = (v->max - v->min) * (1.0/(double)(1<min; r0.x += (ix+0.5)*d.x; @@ -1413,9 +1423,11 @@ class Volume_Node { } // 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 + void createChildren(Volume_Info *v) { // create children to the this node + // ensure there are no children already + if (!child) { + // allocate memory for 8 new children nodes + child = new Volume_Node* [8]; for (int i=0; i<8; i++) { child[i] = new Volume_Node(v, this, i); } @@ -1423,14 +1435,18 @@ class Volume_Node { } // delete children (recursively) - void deleteChildren() { // delete children of this node - if (child) { // if this node has children, delete them + void deleteChildren() { // delete children of this node + // if this node has children, delete them + if (child) { for (int i=0; i<8; i++) { - child[i]->deleteChildren(); // recursive calls to delete all children branches + // recursive calls to delete all children branches + child[i]->deleteChildren(); 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 + // free the memory allocated for the 8 children + delete [] child; + // set children pointer to NULL to indicate that there are no children anymore + child = NULL; } } @@ -1472,19 +1488,18 @@ class Volume_Node { double minz = v->min.z + iz*d.z; int nreg = v->g->regions(); int nsampleAdjusted; - if(avgSigma > 0) { + 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 { + } + else { nsampleAdjusted = nsample; } - + for (long long i=0; irng->getUniform(); @@ -1507,10 +1522,13 @@ class Volume_Node { } // calculate region volumes - void calculateVolumes(Volume_Info *v, vector &volume, vector &sigma, int nreg, double &ntot, vector ®ionHits, vector ®ionTrys) { + void calculateVolumes(Volume_Info *v, vector &volume, + vector &sigma, unsigned int nreg, + EGS_I64 &ntot, vector ®ionHits, + vector ®ionTrys) { if (child) { - for (int i=0; i<8; i++) { + for (unsigned int i=0; i<8; i++) { child[i]->calculateVolumes(v, volume, sigma, nreg, ntot, regionHits, regionTrys); } } @@ -1523,7 +1541,7 @@ class Volume_Node { regionTrys[region] += ntry; } else if (ntry > 0) { - for (int i=0; i