99 EGS_DoseScoring::EGS_DoseScoring(
const string &Name,
102 norm_u(1.0), nreg(0), nmedia(0), max_dreg(-1), max_medl(0),
103 m_lastCase(-1),score_medium_dose(false), score_region_dose(false), output_dose_file(false) {
104 otype =
"EGS_DoseScoring";
107 EGS_DoseScoring::~EGS_DoseScoring() {
125 if (d_regionString.length() > 0 && d_region.size() < 1) {
126 getNumberRegions(d_regionString, d_region);
127 getLabelRegions(d_regionString, d_region);
131 nreg =
app->getnRegions();
133 nmedia =
app->getnMedia();
139 for (imed=0; imed < nmedia; imed++) {
140 sprintf(buf,
"%s%n",
app->getMediumName(imed),&count);
141 if (count > max_medl) {
145 if (d_region.size()>nreg)
146 egsWarning(
"\n*********************************************"
148 "\nRequesting %d dose scoring regions, but there"
149 "\nare only %d geometrical regions!"
150 "\nscoring will be done in all regions!"
151 "\n*********************************************",
152 d_region.size(),nreg);
155 if (score_region_dose) {
156 if (d_region.size() && d_region.size() < nreg) {
158 for (vector<int>::iterator it = d_region.begin(); it < d_region.end(); it++) {
159 if (*it > max_dreg) {
163 for (
int j=0; j<nreg; j++) {
164 d_reg_index.push_back(-1);
165 vol.push_back(vol_list[0]);
167 for (
int i=0; i<d_region.size(); i++) {
168 d_reg_index[d_region[i]]=i;
169 if (i < vol_list.size()) {
170 vol[d_region[i]] = vol_list[i];
173 if (score_region_dose) {
179 for (
int j=0; j<nreg; j++) {
180 d_reg_index.push_back(j);
181 if (j < vol_list.size()) {
182 vol.push_back(vol_list[j]);
185 vol.push_back(vol_list[0]);
188 if (score_region_dose) {
194 if (score_medium_dose) {
197 if (!score_region_dose) {
198 if (vol_list.size() == 1) {
199 vol.push_back(vol_list[0]);
202 for (
int j=0; j<nreg; j++) {
203 if (j < vol_list.size()) {
204 vol.push_back(vol_list[j]);
207 vol.push_back(vol_list[0]);
213 if (output_dose_file) {
215 for (
int i=0; i<nreg; i++) {
216 df_reg.push_back(-1);
224 EGS_Float minx,maxx,miny,maxy,minz,maxz;
225 for (
int k=0; k<nz; k++) {
226 for (
int j=0; j<ny; j++) {
227 for (
int i=0; i<nx; i++) {
237 int g_reg =
app->isWhere(tp);
238 df_reg[g_reg]=i+j*nx+k*nx*ny;
246 description =
"\n*******************************************\n";
250 description +=
"*******************************************\n";
257 description +=
" - Medium dose will be calculated\n";
259 description +=
"\n--------------------------------------\n";
260 sprintf(buf,
"%*s rho (g/cm^3)\n",max_medl,
"medium");
262 description +=
"--------------------------------------\n";
263 for (imed=0; imed < nmedia; imed++) {
266 sprintf(buf,
"%11.8f",
app->getMediumRho(imed));
270 description +=
"--------------------------------------\n";
272 description +=
" Non-unity user-requested normalization = ";
273 sprintf(buf,
"%g\n",norm_u);
276 description +=
"\n*******************************************\n\n";
278 if (d_region.size()) {
298 void EGS_DoseScoring::getNumberRegions(
const string &str, vector<int> ®s) {
302 void EGS_DoseScoring::getLabelRegions(
const string &str, vector<int> ®s) {
306 void EGS_DoseScoring::reportResults() {
307 egsInformation(
"\n======================================================\n");
309 egsInformation(
"======================================================\n");
310 EGS_Float normD = 1., normE=1.;
312 EGS_Float F =
app->getFluence();
316 normD = 1.602e-10*normE;
317 int irmax_digits = getDigits(max_dreg);
318 if (irmax_digits < 2) {
325 egsInformation(
"\n\n==> Summary of region dosimetry (per particle)\n\n");
327 "%*s %*s %12s %12s Edep (MeV) D (Gy) %n\n",
328 irmax_digits,
"ir",max_medl,
"medium",
"rho (g/cm^3)",
"Volume (cm^3)",&count);
331 egsInformation(
"\n==> Summary of region dosimetry (per fluence)\n");
333 "%*s %*s %12s %12s Edep (MeV*cm^2) D (Gy*cm^2) %n\n",
334 irmax_digits,
"ir",max_medl,
"medium",
"rho (g/cm^3)",
"Volume (cm^3)",&count);
336 line.append(count,
'-');
340 for (
int ireg = 0; ireg < nreg; ireg++) {
341 if (d_reg_index[ireg]>=0) {
346 EGS_Float rho =
app->getMediumRho(imed);
347 EGS_Float mass = vol[ireg]*rho;
355 egsInformation(
"%*d %-*s %11.8f %13.6f %12.4e +/- %-7.3f%% %10.4e +/- %-7.3f%%\n",
356 irmax_digits,ireg,max_medl,
app->getMediumName(imed),rho,vol[ireg],r*normE,dr*100.,r*normD/mass,dr*100.);
362 vector<EGS_Float> massM(nmedia,0);
364 for (
int ir=0; ir<nreg; ir++) {
371 EGS_Float volume = vol.size() > 1 ? vol[ir]:vol[0];
372 massM[imed] +=
app->getMediumRho(imed)*volume;
376 egsInformation(
"\n\n==> Summary of media dosimetry (per particle)\n");
378 "%*s %*s Edep/[MeV] D/[Gy] %n\n",
379 max_medl/2,
"medium",max_medl/2,
" ",&count);
382 egsInformation(
"\n\n==> Summary of media dosimetry (per fluence)\n");
384 "%*s %*s Edep/[MeV*cm2] D/[Gy*cm2] %n\n",
385 max_medl/2,
"medium",max_medl/2,
" ",&count);
388 line.append(count,
'-');
391 for (
int im=0; im<nmedia; im++) {
396 "%-*s %10.4e +/- %-7.3f%% %10.4e +/- %-7.3f%%\n",
397 max_medl,
app->getMediumName(im),r*normE,dr*100.,r*normD/massM[im],dr*100.);
404 egsInformation(
"\n EGS_DoseScoring: This is one of a number of parallel jobs. Will only output dose file on combining results.\n");
407 outputDoseFile(normD);
410 egsInformation(
"\n======================================================\n");
413 void EGS_DoseScoring::outputDoseFile(
const EGS_Float &normD) {
416 egsInformation(
"\n EGS_DoseScoring: Writing dose data for EGS_XYZGeometry %s... \n",dose_geom->
getName().c_str());
421 df_out.open(df_name.c_str());
423 egsFatal(
"\n EGS_DoseScoring: Error: Failed to open file %s\n",df_name.c_str());
431 df_out << nx <<
" " << ny <<
" " << nz <<
"\n";
433 float bound,
dose, doseun;
435 for (
int i=0; i<=nx; i++) {
437 df_out << bound <<
" ";
440 for (
int j=0; j<=ny; j++) {
442 df_out << bound <<
" ";
445 for (
int k=0; k<=nz; k++) {
447 df_out << bound <<
" ";
451 for (
int i=0; i<nx*ny*nz; i++) {
453 EGS_Float mass = dose_geom->
getVolume(i)*getRealRho(i);
455 df_out << dose <<
" ";
459 for (
int i=0; i<nx*ny*nz; i++) {
468 df_out << doseun <<
" ";
474 egsFatal(
"\n EGS_DoseScoring: Warning: Dose output file type not recognized.\n");
478 bool EGS_DoseScoring::storeState(ostream &data)
const {
484 if (dose && !dose->storeState(data)) {
496 bool EGS_DoseScoring::setState(istream &data) {
500 if (dose && !dose->setState(data)) {
512 bool EGS_DoseScoring::addState(istream &data) {
513 int err = addTheStates(data);
520 int EGS_DoseScoring::addTheStates(istream &data) {
528 if (!tmp.setState(data)) {
535 if (!tmpM.setState(data)) {
545 if (!tmpF.setState(data)) {
553 void EGS_DoseScoring::resetCounter() {
593 const static char *func =
"createAusgabObject(dose_scoring)";
598 vector <EGS_Float> v_in;
602 vector<string> allowed_mode;
603 allowed_mode.push_back(
"no");
604 allowed_mode.push_back(
"yes");
605 int d_in_medium = input->
getInput(
"medium dose",allowed_mode,0);
606 int d_in_region = input->
getInput(
"region dose",allowed_mode,1);
609 string d_regionsString;
610 vector <int> d_regions;
611 bool using_all_regions=
true;
612 vector <int> d_start, d_stop;
613 if (!input->
getInput(
"dose regions",d_regionsString)&& d_regionsString.length()>0) {
614 using_all_regions =
false;
617 int err1 = input->
getInput(
"dose start region",d_start);
618 int err2 = input->
getInput(
"dose stop region",d_stop);
619 if (!err1 && !err2) {
620 if (d_start.size()==d_stop.size()) {
621 for (
int i=0; i<d_start.size(); i++) {
622 int ir = d_start[i], fr = d_stop[i];
623 for (
int ireg=ir; ireg<=fr; ireg++) {
624 d_regions.push_back(ireg);
627 using_all_regions =
false;
630 "%s: Mismatch in start and stop dose region groups !!!\n"
631 " Calculating dose in ALL regions.\n",func);
634 EGS_Float norma = 1.0;
635 int err04 = input->
getInput(
"normalization",norma);
643 vector <EGS_Float> volin;
645 if (! using_all_regions && v_in.size()== d_start.size()) {
646 for (
int i=0; i<d_start.size(); i++) {
647 int ir = d_start[i], fr = d_stop[i];
648 for (
int ireg=ir; ireg<=fr; ireg++) {
649 volin.push_back(v_in[i]);
658 bool outputdosefile=
false;
665 int err05 = fileinp->
getInput(
"geometry name",gname);
667 egsFatal(
"EGS_DoseScoring: Output dose file: missing/incorrect input for name of geometry.\n");
672 egsFatal(
"EGS_DoseScoring: Output dose file: %s does not name an existing geometry\n",gname.c_str());
674 else if (dgeom->
getType()!=
"EGS_XYZGeometry") {
675 egsFatal(
"EGS_DoseScoring: Output dose file: %s is not an EGS_XYZGeometry.\n",gname.c_str());
679 if (fileinp->
getInput(
"file type", str) < 0) {
683 vector<string> allowed_ftype;
684 allowed_ftype.push_back(
"3ddose");
685 ftype = fileinp->
getInput(
"file type", allowed_ftype, -1);
687 egsFatal(
"EGS_DoseScoring: Output dose file: Invalid file type. Currently only 3ddose is supported.\n");
700 if (volin.size()==1) {
701 result->setVol(volin[0]);
703 else if (volin.size()) {
704 result->setVol(volin);
709 if (!using_all_regions) {
710 if (d_regions.size() > 0) {
711 result->setDoseRegions(d_regions);
714 result->setDoseRegions(d_regionsString);
718 result->setMediumScoring(
true);
721 result->setRegionScoring(
true);
723 if (outputdosefile) {
724 result->setOutputFile(
true,dgeom,ftype);
728 result->setUserNorm(norma);
virtual const string & getType() const =0
Get the geometry type.
A class for scoring an array of quantities (e.g. a dose distribution) in a Monte Carlo simulation...
bool EGS_EXPORT egsGetI64(istream &data, EGS_I64 &n)
Reads a 64 bit integer from the stream data and assigns it to n. Returns true on success, false on failure.
int regions() const
Returns the number of regions (or elements or bins, the most appropriate term depending on the way th...
EGS_ScoringArray * doseM
Scoring dose in each medium.
virtual int getNRegDir(int idir)
string description
A short ausgab object description.
void getNumberRegions(const string &str, vector< int > ®s)
Gets numbers out of str and pushes them onto regs.
void getLabelRegions(const string &str, vector< int > ®s)
Gets the regions for the labels in str and pushes onto regs.
const string & getObjectName() const
Get the object name.
A dose scoring object: header.
EGS_I64 m_lastCase
The event set via setCurrentCase()
A dose scoring ausgab object.
bool EGS_EXPORT egsStoreI64(ostream &data, EGS_I64 n)
Writes the 64 bit integer n to the output stream data and returns true on success, false on failure.
A class representing 3D vectors.
string egsJoinPath(const string &first, const string &second)
Join two path variables (or a path and a file name) using the platform specific directory separator a...
virtual EGS_Float getVolume(int ireg)
Calculates the volume of region ireg.
int getIparallel() const
Returns the job number in a parallel run.
Global egspp functions header file.
bool setState(istream &data)
Sets the state fof the scoring array object from the data in the input stream data.
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
virtual EGS_Float getBound(int idir, int ind)
Returns region boundaries in direction determined by idir.
bool storeState(ostream &data)
Stores the state of the scoring array object into the data stream data.
bool isRealRegion(int ireg)
Returns true if ireg is a real region, false otherwise.
static EGS_BaseGeometry * getGeometry(const string &Name)
Get a pointer to the geometry named Name.
EGS_ScoringArray * dose
Scoring in each dose scoring region.
EGS_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
EGS_ScoringArray * doseF
Scoring dose in each voxel in EGS_XYZGeometry.
EGS_InfoFunction EGS_EXPORT egsInformation
Always use this function for reporting the progress of a simulation and any other type of information...
void reset()
Reset the scoring array to a pristine state.
void currentResult(int ireg, double &r, double &dr)
Sets r to the result in region ireg and dr to its statistical uncertainty.
int getMedium(int ireg)
Returns the medium index in region ireg using C-style indexing.
void setName(EGS_Input *inp)
Set the name of the object from the information provided by inp.
string name
The object name.
const string & getName() const
Get the name of this geometry.
const string & getAppDir() const
Returns the absolute path to the user code directory.
string otype
The object type.
virtual void setApplication(EGS_Application *App)
Set the application this object belongs to.
Base class for advanced EGSnrc C++ applications.
EGS_Application * app
The application this object belongs to.
EGS_InfoFunction EGS_EXPORT egsWarning
Always use this function for reporting warnings.