Skip to content

Enable dose file with transformed geometry #1165

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 7 commits into
base: develop
Choose a base branch
from
Open

Enable dose file with transformed geometry #1165

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
4850d0d
Increase maximum bremsstrahlung splitting to 5000
Jun 7, 2023
013892f
Fix #1131: clarify comment about ms conditional
ftessier May 23, 2024
aad9321
First changes.
Jul 28, 2023
6b7246b
Fix SOME transform errors.
Aug 8, 2023
e6bbeb5
Delete debug output statements.
Aug 9, 2023
aa07bd3
Mods to allow viewing of dose in transformed EGS_XYZGeometry.
Aug 9, 2024
f8df8c5
Formatting changes.
Aug 9, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
92 changes: 59 additions & 33 deletions HEN_HOUSE/egs++/ausgab_objects/egs_dose_scoring/egs_dose_scoring.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,6 +95,7 @@
#include "egs_dose_scoring.h"
#include "egs_input.h"
#include "egs_functions.h"
#include "egs_transformations.h"

EGS_DoseScoring::EGS_DoseScoring(const string &Name,
EGS_ObjectFactory *f) :
Expand Down Expand Up @@ -156,7 +157,7 @@ void EGS_DoseScoring::setApplication(EGS_Application *App) {
if (d_region.size() && d_region.size() < nreg) { // specific dose regions provided
// Get maximum dose scoring region
for (vector<int>::iterator it = d_region.begin(); it < d_region.end(); it++) {
if (*it > max_dreg) {
if ( *it > max_dreg) {
max_dreg = *it;
}
}
Expand Down Expand Up @@ -218,29 +219,42 @@ void EGS_DoseScoring::setApplication(EGS_Application *App) {
//determine the region no. in the EGS_XYZGeometry and corresponding
//global reg. no.
EGS_Vector tp;
int nx=dose_geom->getNRegDir(0);
int ny=dose_geom->getNRegDir(1);
int nz=dose_geom->getNRegDir(2);
EGS_BaseGeometry *d_geom = dose_geom;
vector < EGS_AffineTransform * > t_form;
//now see if this is a transformed EGS_XYZGeometry
//if so, find the base EGS_XYZGeometry and use that to set up the scoring array
while (d_geom->getType().find_last_of("T") == d_geom->getType().length()-1) {
t_form.push_back(d_geom->getTransform()); //have to get transform first
d_geom = d_geom->getBaseGeom();
}
int nx=d_geom->getNRegDir(0);
int ny=d_geom->getNRegDir(1);
int nz=d_geom->getNRegDir(2);
EGS_Float minx,maxx,miny,maxy,minz,maxz;
for (int k=0; k<nz; k++) {
for (int j=0; j<ny; j++) {
for (int i=0; i<nx; i++) {
minx=dose_geom->getBound(0,i);
maxx=dose_geom->getBound(0,i+1);
miny=dose_geom->getBound(1,j);
maxy=dose_geom->getBound(1,j+1);
minz=dose_geom->getBound(2,k);
maxz=dose_geom->getBound(2,k+1);
minx=d_geom->getBound(0,i);
maxx=d_geom->getBound(0,i+1);
miny=d_geom->getBound(1,j);
maxy=d_geom->getBound(1,j+1);
minz=d_geom->getBound(2,k);
maxz=d_geom->getBound(2,k+1);
tp.x=(minx+maxx)/2.;
tp.y=(miny+maxy)/2.;
tp.z=(minz+maxz)/2.;
//now transform tp if this is a transformed EGS_XYZGeometry
//do innermost transformation first
for (int m=t_form.size()-1; m>=0; m--) {
t_form[m]->transform(tp);
}
int g_reg = app->isWhere(tp);
df_reg[g_reg]=i+j*nx+k*nx*ny;
df_reg[g_reg]=i+j *nx+k *nx *ny;
}
}
}
//create an egs_scoring_array of the appropriate size
doseF = new EGS_ScoringArray(nx*ny*nz);
doseF = new EGS_ScoringArray(nx *ny *nz);
}

description = "\n*******************************************\n";
Expand Down Expand Up @@ -313,7 +327,7 @@ void EGS_DoseScoring::reportResults() {
EGS_Float F = app->getFluence();
egsInformation("=> last case = %lld fluence = %g\n", m_lastCase, F);
/* Normalize to actual source fluence */
normE = m_lastCase/F*norm_u;
normE = m_lastCase/F *norm_u;
normD = 1.602e-10*normE;
int irmax_digits = getDigits(max_dreg);
if (irmax_digits < 2) {
Expand Down Expand Up @@ -354,7 +368,7 @@ void EGS_DoseScoring::reportResults() {
dr=1;
}
egsInformation("%*d %-*s %11.8f %13.6f %12.4e +/- %-7.3f%% %10.4e +/- %-7.3f%%\n",
irmax_digits,ireg,max_medl,app->getMediumName(imed),rho,vol[ireg],r*normE,dr*100.,r*normD/mass,dr*100.);
irmax_digits,ireg,max_medl,app->getMediumName(imed),rho,vol[ireg],r *normE,dr*100.,r *normD/mass,dr*100.);
}
}
egsInformation("%s\n",line.c_str());
Expand Down Expand Up @@ -395,7 +409,7 @@ void EGS_DoseScoring::reportResults() {
dr = dr/r;
egsInformation(
"%-*s %10.4e +/- %-7.3f%% %10.4e +/- %-7.3f%%\n",
max_medl,app->getMediumName(im),r*normE,dr*100.,r*normD/massM[im],dr*100.);
max_medl,app->getMediumName(im),r *normE,dr*100.,r *normD/massM[im],dr*100.);
}
}
egsInformation("%s\n",line.c_str());
Expand Down Expand Up @@ -425,39 +439,45 @@ void EGS_DoseScoring::outputDoseFile(const EGS_Float &normD) {
exit(1);
}
//output data
int nx=dose_geom->getNRegDir(0);
int ny=dose_geom->getNRegDir(1);
int nz=dose_geom->getNRegDir(2);
EGS_BaseGeometry *d_geom = dose_geom;
//now see if this is a transformed EGS_XYZGeometry
//if so, find the base EGS_XYZGeometry and use that for data output
while (d_geom->getType().find_last_of("T") == d_geom->getType().length()-1) {
d_geom = d_geom->getBaseGeom();
}
int nx=d_geom->getNRegDir(0);
int ny=d_geom->getNRegDir(1);
int nz=d_geom->getNRegDir(2);
//output no. of voxels in x,y,z
df_out << nx << " " << ny << " " << nz << "\n";
//use single precision real for output
float bound, dose, doseun;
//output voxel boundaries
for (int i=0; i<=nx; i++) {
bound=dose_geom->getBound(0,i);
bound=d_geom->getBound(0,i);
df_out << bound << " ";
}
df_out << "\n";
for (int j=0; j<=ny; j++) {
bound=dose_geom->getBound(1,j);
bound=d_geom->getBound(1,j);
df_out << bound << " ";
}
df_out << "\n";
for (int k=0; k<=nz; k++) {
bound=dose_geom->getBound(2,k);
bound=d_geom->getBound(2,k);
df_out << bound << " ";
}
df_out << "\n";
//divide dose by mass and output
for (int i=0; i<nx*ny*nz; i++) {
for (int i=0; i<nx *ny *nz; i++) {
doseF->currentResult(i,r,dr);
EGS_Float mass = dose_geom->getVolume(i)*getRealRho(i); //local reg.
dose=r*normD/mass;
EGS_Float mass = d_geom->getVolume(i)*getRealRho(i); //local reg.
dose=r *normD/mass;
df_out << dose << " ";
}
df_out << "\n";
//output uncertainties
for (int i=0; i<nx*ny*nz; i++) {
for (int i=0; i<nx *ny *nz; i++) {
doseF->currentResult(i,r,dr);
if (r > 0) {
dr = dr/r;
Expand Down Expand Up @@ -529,24 +549,30 @@ int EGS_DoseScoring::addTheStates(istream &data) {
if (!tmp.setState(data)) {
return 4402;
}
(*dose) += tmp;
( *dose) += tmp;
}
if (doseM) {
EGS_ScoringArray tmpM(nmedia);
if (!tmpM.setState(data)) {
return 4403;
}
(*doseM) += tmpM;
( *doseM) += tmpM;
}
if (doseF) {
int nx=dose_geom->getNRegDir(0);
int ny=dose_geom->getNRegDir(1);
int nz=dose_geom->getNRegDir(2);
EGS_ScoringArray tmpF(nx*ny*nz);
EGS_BaseGeometry *d_geom = dose_geom;
//now see if this is a transformed EGS_XYZGeometry
//if so, find the base EGS_XYZGeometry and use its specs to output data
while (d_geom->getType().find_last_of("T") == d_geom->getType().length()-1) {
d_geom = d_geom->getBaseGeom();
}
int nx=d_geom->getNRegDir(0);
int ny=d_geom->getNRegDir(1);
int nz=d_geom->getNRegDir(2);
EGS_ScoringArray tmpF(nx *ny *nz);
if (!tmpF.setState(data)) {
return 4404;
}
(*doseF) += tmpF;
( *doseF) += tmpF;
}
return 0;
}
Expand Down Expand Up @@ -672,7 +698,7 @@ extern "C" {
if (!dgeom) {
egsFatal("EGS_DoseScoring: Output dose file: %s does not name an existing geometry\n",gname.c_str());
}
else if (dgeom->getType()!="EGS_XYZGeometry") {
else if (dgeom->getType().find("EGS_XYZGeometry")==std::string::npos) {
egsFatal("EGS_DoseScoring: Output dose file: %s is not an EGS_XYZGeometry.\n",gname.c_str());
}
else {
Expand Down
23 changes: 20 additions & 3 deletions HEN_HOUSE/egs++/egs_base_geometry.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,7 @@ using std::vector;

class EGS_Application; // forward declaration
class EGS_Input;
class EGS_AffineTransform; // forward declaration
struct EGS_GeometryIntersections;

#ifdef BPROPERTY64
Expand Down Expand Up @@ -249,6 +250,22 @@ class EGS_EXPORT EGS_BaseGeometry {
return 0.0;
}

/*! \brief Returns base geometry of a transformed geometry

Currently only implemented in EGS_TransformedGeometry
*/
virtual EGS_BaseGeometry *getBaseGeom() {
return NULL;
}

/*! \brief Returns the transform for a transformed geometry

Currently only implemented in EGS_TransformedGeometry
*/
virtual EGS_AffineTransform *getTransform() {
return NULL;
}

/*! Returns number of planar slabs/cylinders/etc in direction idir

Currently only implemented in EGS_XYZGeometry, where idir=0--> X-boundaries,
Expand Down Expand Up @@ -442,7 +459,7 @@ class EGS_EXPORT EGS_BaseGeometry {
virtual EGS_Float getBScaling(int ireg) const {
if (has_Ref_rho && has_B_scaling) {
if (bfactor && ireg >= 0 && ireg < nreg) {
return getMediumRho(medium(ireg))/rhoRef*bfactor[ireg];
return getMediumRho(medium(ireg))/rhoRef *bfactor[ireg];
}
else {
return 1.0;
Expand Down Expand Up @@ -612,9 +629,9 @@ class EGS_EXPORT EGS_BaseGeometry {
*/
virtual bool hasBooleanProperty(int ireg, EGS_BPType prop) const {
if (!bp_array) {
return (prop & bproperty);
return (prop &bproperty);
}
return ireg >= 0 && ireg < nreg ? prop & bp_array[ireg] : false;
return ireg >= 0 && ireg < nreg ? prop &bp_array[ireg] : false;
};

/*! \brief Set the boolean properties of the entire geometry to \a prop
Expand Down
12 changes: 10 additions & 2 deletions HEN_HOUSE/egs++/geometry/egs_gtransformed/egs_gtransformed.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -195,7 +195,7 @@ class EGS_GTRANSFORMED_EXPORT EGS_TransformedGeometry :

EGS_Float howfarToOutside(int ireg, const EGS_Vector &x,
const EGS_Vector &u) {
return ireg >= 0 ? g->howfarToOutside(ireg,x*T,u*T.getRotation()) : 0;
return ireg >= 0 ? g->howfarToOutside(ireg,x *T,u *T.getRotation()) : 0;
};
int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u,
EGS_Float &t, int *newmed=0, EGS_Vector *normal=0) {
Expand All @@ -205,7 +205,7 @@ class EGS_GTRANSFORMED_EXPORT EGS_TransformedGeometry :
int inew = g->howfar(ireg,xt,ut,t,newmed,normal);
//int inew = g->howfar(ireg,x*T,u*T.getRotation(),t,newmed,normal);
if (inew != ireg && normal) {
*normal = T.getRotation()*(*normal);
*normal = T.getRotation()*( *normal);
}
return inew;
};
Expand Down Expand Up @@ -263,6 +263,14 @@ class EGS_GTRANSFORMED_EXPORT EGS_TransformedGeometry :

virtual void getLabelRegions(const string &str, vector<int> &regs);

EGS_BaseGeometry *getBaseGeom() {
return g;
}

EGS_AffineTransform *getTransform() {
return &T;
}

protected:

/*! \brief Don't define media in the transformed geometry definition.
Expand Down
Loading