Skip to content

Commit

Permalink
Improve volume function design
Browse files Browse the repository at this point in the history 
Split printVolumes() into two functions so that
apps can use results of the volume calculation
instead of just printing the output. Also just
tidy the code a bit.
  • Loading branch information
@rtownson @ftessier
rtownson authored and ftessier committed Apr 14, 2021
1 parent 9e3ab84 commit aac8f51
Show file tree
Hide file tree
Showing 2 changed files with 162 additions and 90 deletions.
176 changes: 115 additions & 61 deletions HEN_HOUSE/egs++/egs_base_geometry.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<EGS_Float> boxmin, boxmax;
Expand All @@ -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<std::string> 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<double> sigma, volume, nhit, ntry;
double avgSigma = 0, ntot = 0;
for (int i=0; i<nreg; i++) {
sigma.push_back(0);
volume.push_back(0);
nhit.push_back(0);
ntry.push_back(0);
}
for (long long i=0; i<npass; i++) {
egsInformation("%15d %15d %15d %15d\n", i, root->countLeaves(), nsampleLeaf, root->countLeaves()*nsampleLeaf);

if(i > 0) {
double varSum = 0;
for (int j=0; j<nreg; j++) {
if(volume[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; j<nreg; j++) {
sigma[j] = 0;
volume[j] = 0;
nhit[j] = 0;
ntry[j] = 0;
}

root->calculateVolumes(vInfo, volume, sigma, nreg, ntot, nhit, ntry);

nsampleLeaf = nsample/root->countLeaves() + 1;
if (nsampleLeaf < 10) nsampleLeaf = 10;
vector<double> sigma, volume;
vector<EGS_I64> 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
Expand All @@ -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; i<nreg; i++) {
if (!this->isRealRegion(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);

Expand All @@ -1306,7 +1275,9 @@ void EGS_BaseGeometry::printVolumes(EGS_Input *input, EGS_RandomGenerator *rndm)
sig_med.push_back(0);
}
for (int i=0; i<nreg; i++) {
if (!this->isRealRegion(i)) continue;
if (!this->isRealRegion(i)) {
continue;
}
int med = this->medium(i);
if (med >=0 && med < nmed) {
vol_med[med] += volume[i];
Expand All @@ -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");
Expand All @@ -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");
Expand All @@ -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<EGS_Float> &boxmin,
vector<EGS_Float> &boxmax,
unsigned int &nreg, vector<double> &sigma,
vector<double> &volume,
vector<EGS_I64> &nhit, vector<EGS_I64> &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; i<nreg; i++) {
sigma.push_back(0);
volume.push_back(0);
nhit.push_back(0);
ntry.push_back(0);
}

// Start executing passes
for (EGS_I64 i=0; i<npass; i++) {
if (verbose) {
egsInformation("%15d %15d %15d %15d\n", i, root->countLeaves(), nsampleLeaf, root->countLeaves()*nsampleLeaf);
}

if (i > 0) {
double varSum = 0;
for (unsigned int j=0; j<nreg; j++) {
if (volume[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; j<nreg; j++) {
sigma[j] = 0;
volume[j] = 0;
nhit[j] = 0;
ntry[j] = 0;
}

root->calculateVolumes(vInfo, volume, sigma, nreg, ntot, nhit, ntry);

nsampleLeaf = nsample/root->countLeaves() + 1;
if (nsampleLeaf < 10) {
nsampleLeaf = 10;
}
}

root->deleteChildren();
delete root, vInfo;
Expand Down
Loading

0 comments on commit aac8f51

Please sign in to comment.