Skip to content

Commit

Permalink
Update radionuclide source to rely on base source
Browse files Browse the repository at this point in the history 
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.
  • Loading branch information
@MartinMartinov @ftessier
MartinMartinov authored and ftessier committed Jul 22, 2022
1 parent 8b2fc94 commit c986fd0
Show file tree
Hide file tree
Showing 2 changed files with 78 additions and 366 deletions.
291 changes: 36 additions & 255 deletions HEN_HOUSE/egs++/sources/egs_radionuclide_source/egs_radionuclide_source.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@
#
# Author: Reid Townson, 2016
#
# Contributors:
# Contributors: Martin Martinov
#
###############################################################################
*/
Expand All @@ -41,11 +41,7 @@

EGS_RadionuclideSource::EGS_RadionuclideSource(EGS_Input *input,
EGS_ObjectFactory *f) : EGS_BaseSource(input,f),
sourceType(""), shape(0),
geom(0), regions(0), nrs(0), min_theta(0), max_theta(M_PI),
min_phi(0), max_phi(2*M_PI), gc(IncludeAll),
source_shape(0), target_shape(0), ctry(0), dist(1),
q_allowed(0), decays(0), activity(1) {
baseSource(0), q_allowed(0), decays(0), activity(1), sCount(0) {

int err;
vector<int> tmp_q;
Expand Down Expand Up @@ -73,6 +69,7 @@ EGS_RadionuclideSource::EGS_RadionuclideSource(EGS_Input *input,

// Create the decay spectra
count = 0;
sCount = 0;
Emax = 0;
unsigned int i = 0;
EGS_Float spectrumWeightTotal = 0;
Expand Down Expand Up @@ -104,8 +101,7 @@ EGS_RadionuclideSource::EGS_RadionuclideSource(EGS_Input *input,
++i;
}
if (decays.size() < 1) {
egsWarning("\nEGS_RadionuclideSource: Error: No spectrum of type EGS_RadionuclideSpectrum was defined.\n\n");
return;
egsFatal("\nEGS_RadionuclideSource: Error: No spectrum of type EGS_RadionuclideSpectrum was defined.\n");
}

// Normalize the spectrum weights
Expand Down Expand Up @@ -152,154 +148,25 @@ EGS_RadionuclideSource::EGS_RadionuclideSource(EGS_Input *input,
// Get the active application
app = EGS_Application::activeApplication();

// ==========================
// Source type = isotropic
// ==========================

err = input->getInput("source type",sourceType);
if (err || input->compare(sourceType,"isotropic")) {
sourceType = "isotropic";

egsInformation("EGS_RadionuclideSource: Source type: %s\n",sourceType.c_str());

// Create the shape for source emissions
vector<EGS_Float> pos;
EGS_Input *ishape = input->takeInputItem("shape");
if (ishape) {
shape = EGS_BaseShape::createShape(ishape);
delete ishape;
}
if (!shape) {
string sname;
err = input->getInput("shape name",sname);
if (err)
egsWarning("EGS_RadionuclideSource: missing/wrong inline shape "
"definition and missing wrong 'shape name' input\n");
else {
shape = EGS_BaseShape::getShape(sname);
if (!shape) egsWarning("EGS_RadionuclideSource: a shape named %s"
" does not exist\n");
}
}
string geom_name;
err = input->getInput("geometry",geom_name);
if (!err) {
geom = EGS_BaseGeometry::getGeometry(geom_name);
if (!geom) egsWarning("EGS_RadionuclideSource: no geometry named %s\n",
geom_name.c_str());
else {
vector<string> reg_options;
reg_options.push_back("IncludeAll");
reg_options.push_back("ExcludeAll");
reg_options.push_back("IncludeSelected");
reg_options.push_back("ExcludeSelected");
gc = (GeometryConfinement) input->getInput("region "
"selection",reg_options,0);
if (gc == IncludeSelected || gc == ExcludeSelected) {
vector<int> regs;
err = input->getInput("selected regions",regs);
if (err || regs.size() < 1) {
egsWarning("EGS_RadionuclideSource: region selection %d "
"used but no 'selected regions' input "
"found\n",gc);
gc = gc == IncludeSelected ? IncludeAll : ExcludeAll;
egsWarning(" using %d\n",gc);
}
nrs = regs.size();
regions = new int [nrs];
for (int j=0; j<nrs; j++) {
regions[j] = regs[j];
}
}
}
}
EGS_Float tmp_theta;
err = input->getInput("min theta", tmp_theta);
if (!err) {
min_theta = tmp_theta/180.0*M_PI;
}

err = input->getInput("max theta", tmp_theta);
if (!err) {
max_theta = tmp_theta/180.0*M_PI;
}

err = input->getInput("min phi", tmp_theta);
if (!err) {
min_phi = tmp_theta/180.0*M_PI;
}

err = input->getInput("max phi", tmp_theta);
if (!err) {
max_phi = tmp_theta/180.0*M_PI;
}

buf_1 = cos(min_theta);
buf_2 = cos(max_theta);
// Check for deprecated inputs
string dummy;
err = input->getInput("source type",dummy);
int err2 = input->getInput("geometry",dummy);
if(!err || !err2) {
egsWarning("\nEGS_RadionuclideSource: Warning: Inputs for defining the radionuclide source as an isotropic or collimated source (e.g. 'source type') are deprecated. Please see the documentation - define the type of source separately and then refer to it using the new 'base source' input.\n");
}

// ==========================
// Source type = collimated
// ==========================

else if (input->compare(sourceType,"collimated")) {
sourceType = "collimated";

egsInformation("EGS_RadionuclideSource: Source type: %s\n",sourceType.c_str());

EGS_Input *ishape = input->takeInputItem("source shape");
if (ishape) {
source_shape = EGS_BaseShape::createShape(ishape);
delete ishape;
}
if (!source_shape) {
string sname;
int err = input->getInput("source shape name",sname);
if (err)
egsWarning("EGS_RadionuclideSource: missing/wrong inline source "
"shape definition and missing/wrong 'source shape name' input\n");
else {
source_shape = EGS_BaseShape::getShape(sname);
if (!source_shape)
egsWarning("EGS_RadionuclideSource: a shape named %s"
" does not exist\n",sname.c_str());
}
}
ishape = input->takeInputItem("target shape");
if (ishape) {
target_shape = EGS_BaseShape::createShape(ishape);
delete ishape;
}
if (!target_shape) {
string sname;
int err = input->getInput("target shape name",sname);
if (err)
egsWarning("EGS_RadionuclideSource: missing/wrong inline target"
"shape definition and missing/wrong 'target shape name' input\n");
else {
target_shape = EGS_BaseShape::getShape(sname);
if (!target_shape)
egsWarning("EGS_RadionuclideSource: a shape named %s"
" does not exist\n",sname.c_str());
}
}
if (target_shape) {
if (!target_shape->supportsDirectionMethod())
egsWarning("EGS_RadionuclideSource: the target shape %s, which is"
" of type %s, does not support the getPointSourceDirection()"
" method\n",target_shape->getObjectName().c_str(),
target_shape->getObjectType().c_str());
}
EGS_Float auxd;
int errd = input->getInput("distance",auxd);
if (!errd) {
dist = auxd;
}

// Import base source
err = input->getInput("base source",sName);
if (err) {
egsFatal("\nEGS_RadionuclideSource: Error: Base source must be defined\n"
"using 'base source = some_name'\n");
}
else {
egsFatal("EGS_RadionuclideSource: a source type named %s"
" does not exist\n",sourceType.c_str());
baseSource = EGS_BaseSource::getSource(sName);
if (!baseSource) {
egsFatal("\nEGS_RadionuclideSource: Error: no source named %s"
" is defined\n",sName.c_str());
}

// Initialize emission type to signify nothing has happened yet
Expand Down Expand Up @@ -360,8 +227,10 @@ EGS_I64 EGS_RadionuclideSource::getNextParticle(EGS_RandomGenerator *rndm, int
}

q = decays[i]->getCharge();

getPositionDirection(rndm,x,u,wt);
int qTemp(q), latchTemp(latch);
EGS_Float ETemp(E);
baseSource->getNextParticle(rndm, qTemp, latchTemp, ETemp, wt, x, u);
sCount++;
latch = 0;

if (disintegrationOccurred) {
Expand Down Expand Up @@ -400,103 +269,14 @@ EGS_I64 EGS_RadionuclideSource::getNextParticle(EGS_RandomGenerator *rndm, int
return ++count;
}

void EGS_RadionuclideSource::getPositionDirection(EGS_RandomGenerator *rndm,
EGS_Vector &x, EGS_Vector &u, EGS_Float &wt) {

if (sourceType == "isotropic") {
bool ok = true;
do {
x = shape->getRandomPoint(rndm);
if (geom) {
if (gc == IncludeAll) {
ok = geom->isInside(x);
}
else if (gc == ExcludeAll) {
ok = !geom->isInside(x);
}
else if (gc == IncludeSelected) {
ok = false;
int ireg = geom->isWhere(x);
for (int j=0; j<nrs; ++j) {
if (ireg == regions[j]) {
ok = true;
break;
}
}
}
else {
ok = true;
int ireg = geom->isWhere(x);
for (int j=0; j<nrs; ++j) {
if (ireg == regions[j]) {
ok = false;
break;
}
}
}
}
}
while (!ok);
u.z = buf_1 - rndm->getUniform()*(buf_1 - buf_2);
EGS_Float sinz = 1-u.z*u.z;
if (sinz > epsilon) {
sinz = sqrt(sinz);
EGS_Float cphi, sphi;
EGS_Float phi = min_phi +(max_phi - min_phi)*rndm->getUniform();
cphi = cos(phi);
sphi = sin(phi);
u.x = sinz*cphi;
u.y = sinz*sphi;
}
else {
u.x = 0;
u.y = 0;
}
wt = 1;
}
else if (sourceType == "collimated") {
x = source_shape->getRandomPoint(rndm);
int ntry = 0;
do {
target_shape->getPointSourceDirection(x,rndm,u,wt);
ntry++;
if (ntry > 10000) {
egsFatal("EGS_RadionuclideSource::getPositionDirection:\n"
" my target shape %s, which is of type %s, failed to\n"
" return a positive weight after 10000 attempts\n",
target_shape->getObjectName().c_str(),
target_shape->getObjectType().c_str());
}
}
while (wt <= 0);

if (disintegrationOccurred) {
ctry += ntry-1;
}
}
}

void EGS_RadionuclideSource::setUp() {
otype = "EGS_RadionuclideSource";
if (!isValid()) {
description = "Invalid radionuclide source";
}
else {
description = "Radionuclide source of type ";

if (sourceType == "isotropic") {
description += sourceType;
description += " and a shape of type ";
description += shape->getObjectType();

}
else if (sourceType == "collimated") {
description += sourceType;
description += " from a shape of type ";
description += source_shape->getObjectType();
description += " onto a shape of type ";
description += target_shape->getObjectType();
}
description = "Radionuclide production in base source ";
description += sName;

description += " with:";
if (std::find(q_allowed.begin(), q_allowed.end(), -1) !=
Expand All @@ -513,11 +293,8 @@ void EGS_RadionuclideSource::setUp() {
description += " alphas";
}

if (sourceType == "isotropic") {
if (geom) {
geom->ref();
}
}
description += "\nBase source description:\n";
description += baseSource->getSourceDescription();
}
}

Expand All @@ -528,7 +305,8 @@ bool EGS_RadionuclideSource::storeState(ostream &data_out) const {
}
}
egsStoreI64(data_out,count);
egsStoreI64(data_out,ctry);
egsStoreI64(data_out,sCount);
baseSource->storeState(data_out);

return true;
}
Expand All @@ -542,9 +320,10 @@ bool EGS_RadionuclideSource::addState(istream &data) {
}
EGS_I64 tmp_val;
egsGetI64(data,tmp_val);
count = tmp_val;
count += tmp_val;
egsGetI64(data,tmp_val);
ctry = tmp_val;
sCount += tmp_val;
baseSource->addState(data);

return true;
}
Expand All @@ -555,7 +334,8 @@ void EGS_RadionuclideSource::resetCounter() {
}
ishower = 0;
count = 0;
ctry = 0;
sCount = 0;
baseSource->resetCounter();
}

bool EGS_RadionuclideSource::setState(istream &data) {
Expand All @@ -565,7 +345,8 @@ bool EGS_RadionuclideSource::setState(istream &data) {
}
}
egsGetI64(data,count);
egsGetI64(data,ctry);
egsGetI64(data,sCount);
baseSource->setState(data);

return true;
}
Expand Down
Loading

0 comments on commit c986fd0

Please sign in to comment.