EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_base_geometry.cpp
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ base geometry
5 # Copyright (C) 2015 National Research Council Canada
6 #
7 # This file is part of EGSnrc.
8 #
9 # EGSnrc is free software: you can redistribute it and/or modify it under
10 # the terms of the GNU Affero General Public License as published by the
11 # Free Software Foundation, either version 3 of the License, or (at your
12 # option) any later version.
13 #
14 # EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
15 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
16 # FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for
17 # more details.
18 #
19 # You should have received a copy of the GNU Affero General Public License
20 # along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
21 #
22 ###############################################################################
23 #
24 # Author: Iwan Kawrakow, 2005
25 #
26 # Contributors: Frederic Tessier
27 # Marc Chamberland
28 # Reid Townson
29 # Ernesto Mainegra-Hing
30 # Hugo Bouchard
31 # Hubert Ho
32 #
33 ###############################################################################
34 */
35 
36 
44 #include "egs_base_geometry.h"
45 #include "egs_functions.h"
46 #include "egs_library.h"
47 #include "egs_input.h"
48 #include "egs_application.h"
49 
50 #include <algorithm>
51 #include <vector>
52 #include <cstdio>
53 #include <cstdarg>
54 #include <cstdlib>
55 
56 using namespace std;
57 
58 typedef EGS_BaseGeometry *(*EGS_GeometryCreationFunction)(EGS_Input *);
59 
60 #ifndef SKIP_DOXYGEN
66 class EGS_LOCAL EGS_GeometryPrivate {
67 public:
68 
69  int nnow, ntot;
70  EGS_BaseGeometry **geoms;
71  vector<string> media;
72  vector<EGS_Library *> glibs;
73  static string geom_delimeter;
74  static string libkey;
75  static string create_key;
76  string dso_path;
77  EGS_Application *app;
78 
79  EGS_GeometryPrivate() : nnow(0), ntot(0), geoms(0), app(0) {
80  //egsInformation("EGS_GeometryPrivate() at 0x%x\n",this);
81  setUp();
82  };
83 
84  EGS_GeometryPrivate(const EGS_GeometryPrivate &p) :
85  nnow(0), ntot(0), geoms(0), app(0) {
86  //egsInformation("EGS_GeometryPrivate(0x%x) at 0x%x\n",&p,this);
87  setUp();
88  };
89 
90  void setUp() {
91  char *hhouse = getenv("HEN_HOUSE");
92  if (!hhouse) {
93  egsFatal("Environment variable HEN_HOUSE must be defined\n");
94  }
95  dso_path = hhouse;
96 #if defined WIN32 && !defined CYGWIN
97  char c = '\\';
98 #else
99  char c = '/';
100 #endif
101  if (dso_path[dso_path.size()-1] != c) {
102  dso_path += c;
103  }
104  //dso_path += "geometry"; dso_path += c;
105  dso_path += "egs++";
106  dso_path += c;
107  dso_path += "dso";
108  dso_path += c;
109  dso_path += CONFIG_NAME;
110  };
111 
112  ~EGS_GeometryPrivate() {
113  //egsInformation("Destructing EGS_GeometryPrivate at 0x%x, "
114  // "ntot=%d nnow=%d\n",this,ntot,nnow);
115  if (!ntot) {
116  return;
117  }
118  clearGeometries();
119  delete [] geoms;
120  //egsInformation("Deleting geometry libs\n");
121  {
122  for (unsigned int j=0; j<glibs.size(); j++) {
123  delete glibs[j];
124  }
125  }
126  };
127 
128  void clearGeometries() {
129  media.clear();
130  if (!ntot) {
131  return;
132  }
133  int j = 0, iloop = 0;
134  while (nnow > 0) {
135  if (geoms[j]->deref() == -1) {
136  removeGeometry(geoms[j]);
137  }
138  else {
139  geoms[j++]->ref();
140  }
141  if (j >= nnow && nnow) {
142  j = 0;
143  ++iloop;
144 
145  for (int i=0; i<nnow; ++i) {
146  if (!geoms[i]->deref()) {
147  removeGeometry(geoms[i]);
148  }
149  }
150 
151  if (iloop > 20) {
152 // egsWarning("~EGS_GeometryPrivate(): failed "
153 // "to delete all geometries after 20 loops!\n");
154 
155  break;
156  }
157  }
158  }
159  };
160 
161  void grow(int ngrow) {
162  ntot += ngrow;
163  EGS_BaseGeometry **tmp = new EGS_BaseGeometry* [ntot];
164  for (int j=0; j<nnow; j++) {
165  tmp[j] = geoms[j];
166  }
167  if (geoms) {
168  delete [] geoms;
169  }
170  geoms = tmp;
171  };
172 
173  int addGeometry(EGS_BaseGeometry *g) {
174  if (!g) {
175  return -1;
176  }
177  for (int j=0; j<nnow; j++) {
178  if (geoms[j]->getName() == g->getName()) {
179  return -1;
180  }
181  }
182  if (nnow >= ntot) {
183  grow(10);
184  }
185  geoms[nnow++] = g;
186  return nnow-1;
187 
188  };
189 
190  void removeGeometry(EGS_BaseGeometry *g) {
191  ntot--;
192  EGS_BaseGeometry **tmp = new EGS_BaseGeometry* [ntot];
193  int i=0;
194  for (int j=0; j<nnow; j++) {
195  if (geoms[j] != g) {
196  tmp[i++] = geoms[j];
197  }
198  }
199  if (geoms) {
200  delete [] geoms;
201  }
202  nnow--;
203  geoms = tmp;
204  };
205 
206  EGS_BaseGeometry *getGeometry(const string &name) {
207  for (int j=0; j<nnow; j++)
208  if (geoms[j]->getName() == name) {
209  return geoms[j];
210  }
211  return 0;
212  };
213 
214  int addMedium(const string &Name) {
215  if (EGS_Input::compare(Name,"vacuum")) {
216  return -1;
217  }
218  for (unsigned int j=0; j<media.size(); j++)
219  if (media[j] == Name) {
220  return j;
221  }
222  media.push_back(Name);
223  return media.size()-1;
224  };
225 
226  int getMediumIndex(const string &Name) {
227  for (unsigned int j=0; j<media.size(); j++)
228  if (media[j] == Name) {
229  return j;
230  }
231  return -1;
232  };
233 
234  int nMedia() const {
235  return media.size();
236  };
237 
238  const char *getMediumName(int ind) const {
239  if (ind < 0 || ind > media.size()-1) {
240  return 0;
241  }
242  return media[ind].c_str();
243  };
244 
245  EGS_Float getMediumRho(int ind) const {
246  if (ind==-1) {
247  return -1;
248  }
249  else {
250  if (app) {
251  return app->getMediumRho(ind);
252  }
253  else {
254  return -1;
255  }
256  }
257  };
258 
259  void setApplication(EGS_Application *App) {
260  app = App;
261  };
262 
263  EGS_BaseGeometry *createSingleGeometry(EGS_Input *i);
264 
265 };
266 
267 class EGS_LOCAL EGS_PrivateGeometryLists {
268 public:
269  EGS_PrivateGeometryLists() : nnow(0), ntot(0) {
270  addList(new EGS_GeometryPrivate);
271  };
272  ~EGS_PrivateGeometryLists() {
273  //egsInformation("Deleting geometry lists\n");
274  if (ntot > 0) {
275  for (int j=0; j<nnow; j++) {
276  //egsInformation("Deleting list %d\n",j);
277  delete lists[j];
278  }
279  delete [] lists;
280  }
281  };
282  int size() const {
283  return nnow;
284  };
285  EGS_GeometryPrivate &operator[](int j) {
286  return *lists[j];
287  };
288  void addList(EGS_GeometryPrivate *l) {
289  if (!l) {
290  return;
291  }
292  if (nnow >= ntot) {
293  EGS_GeometryPrivate **tmp = new EGS_GeometryPrivate* [ntot+10];
294  for (int j=0; j<nnow; j++) {
295  tmp[j] = lists[j];
296  }
297  if (ntot > 0) {
298  delete [] lists;
299  }
300  lists = tmp;
301  ntot += 10;
302  }
303  lists[nnow++] = l;
304  };
305  int nnow, ntot;
306  EGS_GeometryPrivate **lists;
307 };
308 
309 string EGS_GeometryPrivate::geom_delimeter = "geometry definition";
310 string EGS_GeometryPrivate::libkey = "library";
311 string EGS_GeometryPrivate::create_key = "createGeometry";
312 #endif
313 
314 int EGS_BaseGeometry::active_glist = 0;
315 
316 static char buf_unique[32];
317 
319 
320 #ifndef SKIP_DOXYGEN
321 EGS_BaseGeometry *EGS_GeometryPrivate::createSingleGeometry(EGS_Input *i) {
322  string libname;
323  if (!i) {
324  egsWarning("createSingleGeometry: null input?\n");
325  return 0;
326  }
327  int error = i->getInput(EGS_GeometryPrivate::libkey,libname);
328  if (error) {
329  egsWarning("createSingleGeometry: input item %s does not define the"
330  " geometry library\n",i->name());
331  return 0;
332  }
333  EGS_Library *lib = 0;
334  for (unsigned int j=0; j<glibs.size(); j++) {
335  if (libname == glibs[j]->libraryName()) {
336  lib = glibs[j];
337  break;
338  }
339  }
340  if (!lib) {
341  lib = new EGS_Library(libname.c_str(),dso_path.c_str());
342  lib->load();
343  if (!lib->isLoaded()) {
344  egsWarning("createSingleGeometry: Failed to load library '%s' from"
345  " %s\n",libname.c_str(),dso_path.c_str());
346  return 0;
347  }
348  glibs.push_back(lib);
349  }
350  EGS_GeometryCreationFunction gcreate = (EGS_GeometryCreationFunction)
351  lib->resolve(EGS_GeometryPrivate::create_key.c_str());
352  if (!gcreate) {
353  egsWarning("createSingleGeometry: failed to resolve the %s function\n"
354  " in geometry library %s\n",
355  EGS_GeometryPrivate::create_key.c_str(),lib->libraryName());
356  return 0;
357  }
358  EGS_BaseGeometry *g = gcreate(i);
359  if (!g) {
360  egsWarning("createSingleGeometry: got null geometry\n");
361  egsWarning(" library: %s\n",lib->libraryName());
362  egsWarning(" input:\n");
363  i->print(4,cerr);
364  return 0;
365  }
366  if (!addGeometry(g)) {
367  egsWarning("createSingleGeometry: failed to add the geometry %s\n"
368  " to the list of geometries. This implies that a geometry with"
369  " this name already exists\n",g->getName().c_str());
370  delete g;
371  return 0;
372  }
373  return g;
374 
375 }
376 #endif
377 
378 //static EGS_LOCAL vector<EGS_GeometryPrivate> egs_geometries;
379 static EGS_LOCAL EGS_PrivateGeometryLists egs_geometries;
380 
381 EGS_Float EGS_BaseGeometry::howfarToOutside(int ireg, const EGS_Vector &x,
382  const EGS_Vector &u) {
383  if (ireg < 0) {
384  return 0;
385  }
386  EGS_Vector xx(x);
387  EGS_Float ttot = 0;
388  for (EGS_I64 loopCount=0; loopCount<=loopMax; ++loopCount) {
389  if (loopCount == loopMax) {
390  egsFatal("EGS_BaseGeometry::howfarToOutside: Too many iterations were required! Input may be invalid, or consider increasing loopMax.");
391  return 0;
392  }
393  EGS_Float t = veryFar;
394  int inew = howfar(ireg,xx,u,t);
395  ttot += t;
396  if (inew < 0) {
397  break;
398  }
399  xx += u*t;
400  ireg = inew;
401  }
402  return ttot;
403 }
404 
406  int n = egs_geometries.size();
407  //egsInformation("EGS_BaseGeometry::setActiveGeometryList: size=%d list=%d\n",
408  // n,list);
409  for (int j=n; j<=list; j++)
410  //egs_geometries.push_back(EGS_GeometryPrivate());
411  {
412  egs_geometries.addList(new EGS_GeometryPrivate);
413  }
414  active_glist = list;
415 }
416 
417 /*
418 extern "C" void __list_geometries() {
419  egsInformation("The following %d geometries are defined:\n",
420  egs_geometries[active_glist].nnow);
421  for(int j=0; j<egs_geometries[active_glist].nnow; j++) {
422  egsInformation("%d %s %s\n",j,
423  egs_geometries[active_glist].geoms[j]->getType().c_str(),
424  egs_geometries[active_glist].geoms[j]->getName().c_str());
425  }
426 }
427 */
428 
429 EGS_BaseGeometry::EGS_BaseGeometry(const string &Name) : nreg(0), name(Name),
430  region_media(0), med(-1), has_rho_scaling(false), rhor(0),
431  has_B_scaling(false), has_Ref_rho(false), bfactor(0), rhoRef(1.0),
432  nref(0), debug(false), is_convex(true), bproperty(0), bp_array(0),
433  boundaryTolerance(distanceEpsilon) {
434 
435  halfBoundaryTolerance = boundaryTolerance/2.;
436  if (!egs_geometries.size()) {
437  egs_geometries.addList(new EGS_GeometryPrivate);
438  }
439  if (!name.size()) {
440  name = getUniqueName();
441  }
442  if (egs_geometries[active_glist].addGeometry(this) < 0)
443  egsFatal("EGS_BaseGeometry::EGS_BaseGeometry:\n"
444  " a geometry with name %s already exists\n",name.c_str());
445 }
446 
448  if (region_media) {
449  delete [] region_media;
450  }
451  if (rhor && has_rho_scaling) {
452  delete [] rhor;
453  }
454  if (bp_array) {
455  delete [] bp_array;
456  }
457  if (bfactor && has_B_scaling) {
458  delete [] bfactor;
459  }
460  //egsInformation("Deleting geometry at 0x%x, list=%d\n",this,active_glist);
461  egs_geometries[active_glist].removeGeometry(this);
462 }
463 
465  egs_geometries[active_glist].clearGeometries();
466 }
467 
469  return egs_geometries[active_glist].getGeometry(Name);
470 }
471 
472 EGS_BaseGeometry **EGS_BaseGeometry::getGeometries() {
473  return egs_geometries[active_glist].geoms;
474 }
475 
476 int EGS_BaseGeometry::getNGeometries() {
477  return egs_geometries[active_glist].nnow;
478 }
479 
480 void EGS_BaseGeometry::setMedium(const string &Name) {
481  med = egs_geometries[active_glist].addMedium(Name);
482  if (region_media)
483  for (int j=0; j<nreg; j++) {
484  region_media[j] = med;
485  }
486 }
487 
488 int EGS_BaseGeometry::addMedium(const string &medname) {
489  return egs_geometries[active_glist].addMedium(medname);
490 }
491 
492 int EGS_BaseGeometry::getMediumIndex(const string &medname) {
493  return egs_geometries[active_glist].getMediumIndex(medname);
494 }
495 
496 void EGS_BaseGeometry::setMedium(int istart, int iend, const string &Name,
497  int delta) {
498  int imed = egs_geometries[active_glist].addMedium(Name);
499  setMedium(istart,iend,imed,delta);
500 }
501 
502 void EGS_BaseGeometry::setMedium(int istart, int iend, int imed, int delta) {
503  if (nreg <= 1) {
504  med = imed;
505  return;
506  }
507  if (delta <= 0) {
508  return;
509  }
510  if (istart < 0) {
511  istart = 0;
512  }
513  if (iend > nreg-1) {
514  iend = nreg-1;
515  }
516  if (!region_media) {
517  region_media = new short [nreg];
518  for (int j=0; j<nreg; j++) {
519  region_media[j] = med;
520  }
521  }
522  for (int j=istart; j<=iend; j+=delta) {
523  region_media[j] = imed;
524  }
525 }
526 
528  return egs_geometries[active_glist].nMedia();
529 }
530 
531 const char *EGS_BaseGeometry::getMediumName(int ind) {
532  return egs_geometries[active_glist].getMediumName(ind);
533 }
534 
535 EGS_Float EGS_BaseGeometry::getMediumRho(int ind) const {
536  return egs_geometries[active_glist].getMediumRho(ind);
537 }
538 
539 void EGS_BaseGeometry::setApplication(EGS_Application *App) {
540  return egs_geometries[active_glist].setApplication(App);
541 }
542 
544  return egs_geometries[active_glist].createSingleGeometry(input);
545 }
546 
548  EGS_Input *ginput = input;
549  bool delete_it = false;
550  if (!input->isA(egs_geometries[active_glist].geom_delimeter)) {
551  ginput = input->takeInputItem(egs_geometries[active_glist].geom_delimeter);
552  delete_it = true;
553  }
554  if (!ginput) {
555  egsWarning("EGS_BaseGeometry::createGeometry: no geometry specification"
556  " in this input\n");
557  return 0;
558  }
559  EGS_Input *ij;
560  bool error = false;
561  while ((ij = ginput->takeInputItem("geometry")) != 0) {
562  EGS_BaseGeometry *g = egs_geometries[active_glist].createSingleGeometry(ij);
563  if (!g) {
564  error = true;
565  }
566  delete ij;
567  }
568  // Check to make sure that geometries have unique names
569  for (int j=0; j<egs_geometries[active_glist].nnow; j++) {
570  string gname = egs_geometries[active_glist].geoms[j]->getName();
571  for (int k=0; k<egs_geometries[active_glist].nnow; k++) {
572  if (k == j) {
573  continue;
574  }
575  if (gname == egs_geometries[active_glist].geoms[k]->getName()) {
576  egsFatal("\ncreateGeometry: Error: multiple geometries with"
577  " the same name exist: %s\n\n", gname.c_str());
578  return 0;
579  }
580  }
581  }
582  if (error) {
583  egsFatal("EGS_BaseGeometry::createGeometry: errors during geometry"
584  " definition\n");
585  return 0;
586  }
587  string sim_geom;
588  int err = ginput->getInput("simulation geometry",sim_geom);
589  if (err) {
590  egsWarning("EGS_BaseGeometry::createGeometry: missing/wrong keyword"
591  " 'simulation geometry'\n");
592  return 0;
593  }
594  EGS_BaseGeometry *g = egs_geometries[active_glist].getGeometry(sim_geom);
595  if (!g) egsWarning("EGS_BaseGeometry::createGeometry: a geometry with "
596  "the name %s does not exist\n",sim_geom.c_str());
597  if (delete_it) {
598  delete ginput;
599  }
600  return g;
601 }
602 
604  sprintf(buf_unique,"geometry%d",egs_geometries[active_glist].nnow);
605  string result(buf_unique);
606  return result;
607 }
608 
610  int err = i->getInput("name",name);
611  if (err) {
612  name = getUniqueName();
613  }
614  EGS_Input *inp;
615  int irep=0;
616  while ((inp = i->takeInputItem("replica"))) {
617  string typ;
618  int ncopy;
619  vector<EGS_Float> trans, trans_o;
620  vector<EGS_Float> rot_axis;
621  EGS_Float rot_angle, rot_angle_o;
622  int err1 = inp->getInput("type",typ);
623  int err2 = inp->getInput("number of copies",ncopy);
624  int err3 = inp->getInput("translation delta",trans);
625  int err4 = inp->getInput("first translation",trans_o);
626  int err5 = inp->getInput("rotation axis",rot_axis);
627  int err6 = inp->getInput("rotation delta",rot_angle);
628  int err7 = inp->getInput("first rotation",rot_angle_o);
629  bool do_it = true;
630  int ttype;
631  if (err1 || err2) {
632  if (err1) egsWarning("geometry replication: 'type' not defined ->"
633  " ignoring input\n");
634  if (err2) egsWarning("geometry replication: 'number of copies' "
635  "not defined -> ignoring input\n");
636  do_it = false;
637  }
638  else {
639  if (ncopy < 1) {
640  egsWarning("geometry replication: %d copies?\n",ncopy);
641  do_it = false;
642  }
643  if (typ == "line") {
644  if (err3 || trans.size() != 3) {
645  egsWarning("geometry replication: got %d inputs for "
646  "'translation', need 3\n",trans.size());
647  do_it = false;
648  }
649  else {
650  if (err4) {
651  trans_o.push_back(0);
652  trans_o.push_back(0);
653  trans_o.push_back(0);
654  }
655  }
656  ttype = 0;
657  }
658  else if (typ == "rotation") {
659  if (err5 || rot_axis.size() != 3) {
660  egsWarning("geometry replication: got %d inputs for "
661  "'rotation axis', need 3\n",rot_axis.size());
662  do_it = false;
663  }
664  if (err6) {
665  egsWarning("geometry replication: missing 'rotation delta'"
666  " input\n");
667  do_it = false;
668  }
669  if (err7) {
670  rot_angle_o = 0;
671  }
672  ttype = 1;
673  }
674  else {
675  egsWarning("geometry replication: unknown replica type %s\n",
676  typ.c_str());
677  do_it = false;
678  }
679  }
680  if (do_it) {
681  ++irep;
682  char buf[1024];
683  for (int icopy=1; icopy<=ncopy; icopy++) {
684  //string content(":start geometry:\n");
685  string content;
686  content += " library = egs_gtransformed\n";
687  sprintf(buf," name = %s_rep%d_%d\n",name.c_str(),irep,icopy);
688  content += buf;
689  sprintf(buf," my geometry = %s\n",name.c_str());
690  content += buf;
691  content += " :start transformation:\n";
692  if (ttype == 0)
693  sprintf(buf," translation = %g %g %g\n",
694  trans_o[0]+trans[0]*icopy,
695  trans_o[1]+trans[1]*icopy,
696  trans_o[2]+trans[2]*icopy);
697  else
698  sprintf(buf," rotation = %g %g %g %g\n",
699  rot_axis[0],rot_axis[1],rot_axis[2],
700  rot_angle_o+rot_angle*icopy);
701  content += buf;
702  content += " :stop transformation:\n";
703  //content += ":stop geometry:\n";
704  EGS_Input aux;
705  aux.setContentFromString(content);
706  EGS_BaseGeometry *g =
708  if (!g) egsWarning("geometry replication: failed to create"
709  " replica %d of %s\n",icopy,name.c_str());
710  //else egsInformation("geometry replication: created replica"
711  // " %s\n",g->getName().c_str());
712  }
713  }
714  }
715 }
716 
718  int err = i->getInput("boundary tolerance", boundaryTolerance);
719  if (err > 0) {
720  egsWarning("EGS_BaseGeometry::setBoundaryTolerance(): error while reading 'boundary tolerance' input\n");
721  return;
722  }
723  halfBoundaryTolerance = boundaryTolerance/2.;
724 }
725 
727  egsInformation("======================== geometry =====================\n");
728  egsInformation(" type = %s\n",getType().c_str());
729  egsInformation(" name = %s\n",getName().c_str());
730  egsInformation(" number of regions = %d\n",nreg);
731  if (hasBScaling()) {
732  egsInformation("\nB scaling ON\n");
733  }
734 }
735 
737  egsInformation("\nThe following geometries are defined:\n\n");
738  for (int j=0; j<egs_geometries[active_glist].nnow; j++) {
739  egs_geometries[active_glist].geoms[j]->printInfo();
740  }
741 }
742 
744  EGS_Input *input = inp;
745  bool delete_it = false;
746  if (!input->isA("media input")) {
747  input = inp->takeInputItem("media input");
748  if (!input) {
749  return;
750  }
751  // i.e., if there is no media related input for this geometry,
752  // we don't warn as we assume that media will be set from the outside
753  delete_it = true;
754  }
755  vector<string> media_names;
756  int err = input->getInput("media",media_names);
757  int *med_ind = 0;
758  int nmed = media_names.size();
759  if (!err && nmed > 0) {
760  med_ind = new int [nmed];
761  for (int j=0; j<nmed; j++) {
762  med_ind[j] = egs_geometries[active_glist].addMedium(media_names[j]);
763  }
764  }
765  else {
766  nmed = nMedia();
767  if (nmed < 1) {
768  if (delete_it) {
769  delete input;
770  }
771  return;
772  }
773  med_ind = new int [nmed];
774  for (int j=0; j<nmed; j++) {
775  med_ind[j] = j;
776  }
777  }
778  setMedia(input,nmed,med_ind);
779  setRelativeRho(input);
780  setBScaling(input);
781  delete [] med_ind;
782  if (delete_it) {
783  delete input;
784  }
785 }
786 
787 void EGS_BaseGeometry::setMedia(EGS_Input *input, int nmed, const int *mind) {
788  EGS_Input *i;
789  med = mind[0];
790  while ((i = input->takeInputItem("set medium"))) {
791  vector<int> inp;
792  int err = i->getInput("set medium",inp);
793  delete i;
794  if (!err) {
795  //if( inp.size() == 2 ) setMedium(inp[0],inp[0]+1,mind[inp[1]]);
796  if (inp.size() == 2) {
797  setMedium(inp[0],inp[0],mind[inp[1]]);
798  }
799  else if (inp.size() == 3) {
800  setMedium(inp[0],inp[1],mind[inp[2]]);
801  }
802  else if (inp.size() == 4) {
803  setMedium(inp[0],inp[1],mind[inp[2]],inp[3]);
804  }
805  else egsWarning("EGS_BaseGeometry::setMedia(): found %d inputs\n"
806  "in a 'set medium' input. 2 or 3 are allowed\n",inp.size());
807  }
808  else egsWarning("EGS_BaseGeometry::setMedia(): wrong 'set medium'"
809  " input\n");
810  }
811 }
812 
813 void EGS_BaseGeometry::setRelativeRho(int start, int end, EGS_Float rho) {
814  if (start < 0) {
815  start = 0;
816  }
817  if (end >= nreg) {
818  end = nreg-1;
819  }
820  if (end >= start) {
821  int j;
822  if (!rhor) {
823  rhor = new EGS_Float [nreg];
824  for (j=0; j<nreg; j++) {
825  rhor[j] = 1;
826  }
827  }
828  for (j=start; j<=end; j++) {
829  rhor[j] = rho;
830  }
831  has_rho_scaling = true;
832  }
833 }
834 
836  EGS_Input *i;
837  while ((i = input->takeInputItem("set relative density"))) {
838  vector<EGS_Float> tmp;
839  int err = i->getInput("set relative density",tmp);
840  if (!err) {
841  if (tmp.size() == 2) {
842  int start = (int)(tmp[0]+0.1);
843  int end = start;
844  setRelativeRho(start,end,tmp[1]);
845  }
846  else if (tmp.size() == 3) {
847  int start = (int)(tmp[0]+0.1);
848  int end = (int)(tmp[1]+0.1);
849  setRelativeRho(start,end,tmp[2]);
850  }
851  else {
852  egsWarning("EGS_BaseGeometry::setRelativeRho(): found %d "
853  "inputs in a 'set relative density' input.\n",tmp.size());
854  egsWarning(" 2 or 3 are allowed => input ignored\n");
855  }
856  }
857  delete i;
858  }
859 }
860 
861 void EGS_BaseGeometry::setBScaling(int start, int end, EGS_Float bf) {
862  if (start < 0) {
863  start = 0;
864  }
865  if (end >= nreg) {
866  end = nreg-1;
867  }
868  if (end >= start) {
869  int j;
870  if (!bfactor) {
871  bfactor = new EGS_Float [nreg];
872  for (j=0; j<nreg; j++) {
873  bfactor[j] = 1.0;
874  }
875  }
876  for (j=start; j<=end; j++) {
877  bfactor[j] = bf;
878  }
879  has_B_scaling = true;
880  }
881 }
882 
884  EGS_Input *i;
885  /* Check whether scaling of the B field requested */
886  while ((i = input->takeInputItem("set B scaling"))) {
887  vector<EGS_Float> tmp;
888  int err = i->getInput("set B scaling",tmp);
889  if (!err) {
890  if (tmp.size() == 2) {
891  int start = (int)(tmp[0]+0.1);
892  int end = start;
893  setBScaling(start, end, tmp[1]);
894  }
895  else if (tmp.size() == 3) {
896  int start = (int)(tmp[0]+0.1);
897  int end = (int)(tmp[1]+0.1);
898  setBScaling(start, end, tmp[2]);
899  }
900  else {
901  egsWarning("EGS_BaseGeometry::setBScaling(): found %d "
902  "inputs in a 'set B scaling' input.\n", tmp.size());
903  egsWarning(" 2 or 3 are allowed => input ignored\n");
904  }
905  }
906  delete i;
907  }
908  /* Check whether mass density scaling of the B field requested */
909  EGS_Float refD;
910  int err0 = input->getInput("B scaling reference density", refD);
911  if (!err0) {
912  rhoRef = refD;
913  has_Ref_rho = true;
914  }
915 
916 }
917 
919  const EGS_Vector &u, EGS_GeometryIntersections *isections) {
920  if (n < 1) {
921  return -1;
922  }
923  int ifirst = 0;
924  EGS_Float t, ttot = 0;
925  EGS_Vector x(X);
926  int imed;
927  if (ireg < 0) {
928  t = veryFar;
929  ireg = howfar(ireg,x,u,t,&imed);
930  if (ireg < 0) {
931  return 0;
932  }
933  isections[0].t = t;
934  isections[0].rhof = 1;
935  isections[0].ireg = -1;
936  isections[0].imed = -1;
937  ttot = t;
938  ++ifirst;
939  x += u*t;
940  }
941  else {
942  imed = medium(ireg);
943  }
944 
945  for (int j=ifirst; j<n; j++) {
946  isections[j].imed = imed;
947  isections[j].rhof = getRelativeRho(ireg);
948  isections[j].ireg = ireg;
949  t = veryFar;
950  int inew = howfar(ireg,x,u,t,&imed);
951  ttot += t;
952  isections[j].t = ttot;
953  if (inew < 0 || inew == ireg) {
954  return j+1;
955  }
956  ireg = inew;
957  x += u*t;
958  }
959 
960  return ireg >= 0 ? -1 : n;
961 
962 }
963 
965  bproperty = prop;
966  if (bp_array) {
967  delete [] bp_array;
968  bp_array = 0;
969  }
970 }
971 
973  if (bit < 0 || bit >= 8*sizeof(EGS_BPType)) {
974  egsWarning("EGS_BaseGeometry::addBooleanProperty: attempt to set the "
975  "%d'th bith!\n",bit);
976  return;
977  }
978  EGS_BPType prop = 1 << bit;
979  bproperty |= prop;
980  if (bp_array) {
981  for (int j=0; j<nreg; j++) {
982  bp_array[j] |= prop;
983  }
984  }
985 }
986 
987 void EGS_BaseGeometry::setBooleanProperty(EGS_BPType prop, int start, int end,
988  int step) {
989  if (start < 0) {
990  start = 0;
991  }
992  if (end >= nreg) {
993  end = nreg-1;
994  }
995  if (start == 0 && end == nreg-1 && step==1) {
996  setBooleanProperty(prop);
997  }
998  else {
999  if (!bp_array) {
1000  bp_array = new EGS_BPType [nreg];
1001  for (int j=0; j<nreg; j++) {
1002  bp_array[j] = 0;
1003  }
1004  }
1005  for (int j=start; j<=end; j+=step) {
1006  bp_array[j] = prop;
1007  }
1008  }
1009 }
1010 
1011 void EGS_BaseGeometry::addBooleanProperty(int bit, int start, int end,
1012  int step) {
1013  if (bit < 0 || bit >= 8*sizeof(EGS_BPType)) {
1014  egsWarning("EGS_BaseGeometry::addBooleanProperty: attempt to set the "
1015  "%d'th bith!\n",bit);
1016  return;
1017  }
1018  if (start < 0) {
1019  start = 0;
1020  }
1021  if (end >= nreg) {
1022  end = nreg-1;
1023  }
1024  if (start == 0 && end == nreg-1 && step==1) {
1025  addBooleanProperty(bit);
1026  }
1027  else {
1028  EGS_BPType prop = 1 << bit;
1029  if (!bp_array) {
1030  bp_array = new EGS_BPType [nreg];
1031  for (int j=0; j<nreg; j++) {
1032  bp_array[j] = 0;
1033  }
1034  }
1035  for (int j=start; j<=end; j+=step) {
1036  bp_array[j] |= prop;
1037  }
1038  }
1039 }
1040 
1041 // Gets region numbers from a string
1042 // Pushes the regions onto the array regs
1043 void EGS_BaseGeometry::getNumberRegions(const string &str, vector<int> &regs) {
1044 
1045  if (!str.empty()) {
1046 
1047  // Tokenize the input string
1048  vector<string> tokens;
1049  const char *ptr = str.c_str();
1050  do {
1051  const char *begin = ptr;
1052  while (*ptr != ' ' && *ptr) {
1053  ptr++;
1054  }
1055  tokens.push_back(string(begin, ptr));
1056  }
1057  while (*ptr++ != '\0');
1058 
1059  for (int i=0; i<tokens.size(); i++) {
1060  // Search for tokens that are numbers, not strings
1061  // Push the region numbers onto the regions array
1062  if (tokens[i].find_first_not_of(" -0123456789") == std::string::npos) {
1063  regs.push_back(atoi(tokens[i].c_str()));
1064  }
1065  }
1066  }
1067 }
1068 
1069 void EGS_BaseGeometry::getLabelRegions(const string &str, vector<int> &regs) {
1070 
1071  // Tokenize the input string - this allows for multiple labels
1072  vector<string> tokens;
1073  const char *ptr = str.c_str();
1074  do {
1075  const char *begin = ptr;
1076  while (*ptr != ' ' && *ptr) {
1077  ptr++;
1078  }
1079  tokens.push_back(string(begin, ptr));
1080  }
1081  while (*ptr++ != '\0');
1082 
1083  // Get all regions lists for this named label
1084  for (int j=0; j<tokens.size(); j++) {
1085  for (int i=0; i<labels.size(); i++) {
1086  if (labels[i].name.compare(tokens[j]) == 0) {
1087  regs.insert(regs.end(), labels[i].regions.begin(), labels[i].regions.end());
1088  }
1089  }
1090  }
1091 
1092  // Sort region list and remove duplicates
1093  sort(regs.begin(), regs.end());
1094  regs.erase(unique(regs.begin(), regs.end()), regs.end());
1095 }
1096 
1097 
1099  EGS_Input *i;
1100  int labelCount=0;
1101  while ((i = input->takeInputItem("set label"))) {
1102 
1103  // get input string
1104  string inp;
1105  int err = i->getInput("set label",inp);
1106  delete i;
1107 
1108  // bail out on read error
1109  if (err) {
1110  egsWarning("EGS_BaseGeometry::setLabels(): error while reading 'set label' input\n");
1111  return 0;
1112  }
1113 
1114  // set the labels from input string
1115  int nLabel = setLabels(inp);
1116 
1117  // increase label count
1118  labelCount += nLabel;
1119  }
1120 
1121  return labelCount;
1122 }
1123 
1124 
1125 int EGS_BaseGeometry::setLabels(const string &inp) {
1126 
1127  // tokenize input string
1128  vector<string> tokens;
1129  const char *ptr = inp.c_str();
1130  do {
1131  const char *begin = ptr;
1132  while (*ptr != ' ' && *ptr) {
1133  ptr++;
1134  }
1135  tokens.push_back(string(begin, ptr));
1136  }
1137  while (*ptr++ != '\0');
1138 
1139  // bail out if there are no label tokens
1140  if (tokens.size() < 1) {
1141  egsWarning("EGS_BaseGeometry::setLabels(): no label name\n");
1142  return 0;
1143  }
1144 
1145  // parse label into a label class
1146  label lab;
1147  lab.name = tokens[0];
1148  for (int i=1; i<tokens.size(); i++) {
1149  int reg = atoi(tokens[i].c_str());
1150  if (reg < nreg && reg >= 0) {
1151  lab.regions.push_back(reg);
1152  }
1153  else {
1154  egsWarning("EGS_BaseGeometry::setLabels(): label \"%s\": region %d is beyond the number " \
1155  "of regions in this geometry\n", lab.name.c_str(), reg);
1156  }
1157  }
1158 
1159  // continue if there is no region
1160  if (lab.regions.size() <= 0) {
1161  return 0;
1162  }
1163 
1164  // sort region list and remove duplicates
1165  sort(lab.regions.begin(), lab.regions.end());
1166  lab.regions.erase(unique(lab.regions.begin(), lab.regions.end()), lab.regions.end());
1167 
1168  // push current label onto vector of labels
1169  labels.push_back(lab);
1170 
1171  return 1;
1172 }
Base class for advanced EGSnrc C++ applications.
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
virtual int computeIntersections(int ireg, int n, const EGS_Vector &x, const EGS_Vector &u, EGS_GeometryIntersections *isections)
Calculates intersection distances to region boundaries.
static int nMedia()
Get the number of media registered so far by all geometries.
EGS_BPType bproperty
A bit mask of boolean properties for the entire geometry.
virtual EGS_Float howfarToOutside(int ireg, const EGS_Vector &x, const EGS_Vector &u)
virtual void setBScaling(int start, int end, EGS_Float bf)
Set the B field scaling factor in regions.
void setMedia(EGS_Input *inp)
Set the media in the geometry from the input pointed to by inp.
static EGS_BaseGeometry * createSingleGeometry(EGS_Input *inp)
Create a single geometry from the input inp.
virtual const string & getType() const =0
Get the geometry type.
virtual void getLabelRegions(const string &str, vector< int > &regs)
Get the list of all regions labeled with str.
static void describeGeometries()
Describes all existing geometries.
int nreg
Number of local regions in this geometry.
virtual int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u, EGS_Float &t, int *newmed=0, EGS_Vector *normal=0)=0
Calculate the distance to a boundary from x along the direction u.
bool has_B_scaling
Does this geometry has B field scaling factor?
bool has_rho_scaling
Does this geometry have relative mass density scvaling?
EGS_Float boundaryTolerance
Boundary tolerance for geometries that need it.
bool hasBScaling() const
Does this geometry object have a B field scaling feature?
short * region_media
Array of media indeces.
void setName(EGS_Input *inp)
Set the name of the geometry from the input inp.
static void setActiveGeometryList(int list)
Set the currently active geometry list.
void setMedium(const string &Name)
Set all regions to a medium with name Name.
static EGS_BaseGeometry * createGeometry(EGS_Input *)
Create a geometry (or geometries) from a given input.
const string & getName() const
Get the name of this geometry.
EGS_Float * rhor
Array with relative mass densities.
virtual EGS_Float getRelativeRho(int ireg) const
Get the relative mass density in region ireg.
EGS_Float * bfactor
Array with B field scaling factors.
int med
Medium index.
virtual void setBooleanProperty(EGS_BPType prop)
Set the boolean properties of the entire geometry to prop.
static string getUniqueName()
Get a unique geometry name.
static void clearGeometries()
Clears (deletes) all geometries in the currently active geometry list.
virtual int medium(int ireg) const
Returns the medium index in region ireg.
static int error_flag
Set to non-zero status if a geometry problem is encountered.
string name
Name of this geometry.
int setLabels(EGS_Input *input)
Set the labels from an input block.
virtual void setRelativeRho(int start, int end, EGS_Float rho)
Set the relative mass density in regions.
vector< label > labels
Labels.
virtual void printInfo() const
Print information about this geometry.
int ref()
Increase the reference count to this geometry.
static EGS_BaseGeometry * getGeometry(const string &Name)
Get a pointer to the geometry named Name.
virtual void addBooleanProperty(int bit)
Add a boolean property for the entire geometry by setting the bit'th bit.
static int addMedium(const string &medname)
Add a medium or get the index of an existing medium.
EGS_Float rhoRef
Reference density for B field scaling.
static int getMediumIndex(const string &medname)
Get the index of a medium named medname.
virtual ~EGS_BaseGeometry()
Destructor.
EGS_BaseGeometry(const string &Name)
Construct a geometry named Name.
static const char * getMediumName(int ind)
Get the name of medium with index ind.
void setBoundaryTolerance(EGS_Input *inp)
Set the value of the boundary tolerance from the input inp.
EGS_BPType * bp_array
An array of boolean properties on a region by region basis.
virtual void getNumberRegions(const string &str, vector< int > &regs)
Get a list of all the regions labeled with a number.
A class for storing information in a tree-like structure of key-value pairs. This class is used throu...
Definition: egs_input.h:182
EGS_Input * takeInputItem(const string &key, bool self=true)
Get the property named key.
Definition: egs_input.cpp:226
void print(int nind, ostream &)
Used for debugging purposes.
Definition: egs_input.cpp:1174
const char * name() const
Get the name of this property.
Definition: egs_input.cpp:271
int setContentFromString(string &input)
Definition: egs_input.cpp:206
static bool compare(const string &s1, const string &s2)
Definition: egs_input.cpp:1170
bool isA(const string &key) const
Definition: egs_input.cpp:278
int getInput(const string &key, vector< string > &values) const
Assign values to an array of strings from an input identified by key.
Definition: egs_input.cpp:338
A class for dynamically loading shared libraries.
Definition: egs_library.h:52
const char * libraryName() const
Returns the name of the library object as given in the constructor.
bool isLoaded() const
Returns true if the library is loaded, false otherwise.
bool load()
Loads the library.
void * resolve(const char *func)
Returns the address of the exported symbol func.
A class representing 3D vectors.
Definition: egs_vector.h:57
EGS_Application class header file.
EGS_BaseGeometry class header file.
Global egspp functions header file.
EGS_Input class header file.
EGS_Library class header file.
EGS_InfoFunction EGS_EXPORT egsInformation
Always use this function for reporting the progress of a simulation and any other type of information...
EGS_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
const EGS_Float distanceEpsilon
The distanceEpsilon constant for physical distance comparisons.
Definition: egs_functions.h:86
const EGS_I64 loopMax
The maximum number of iterations for near-infinite loops.
Definition: egs_functions.h:96
EGS_InfoFunction EGS_EXPORT egsWarning
Always use this function for reporting warnings.
const EGS_Float veryFar
A very large float.
EGS_I32 ireg
region index
EGS_Float t
distance to next region boundary
EGS_Float rhof
relative mass density in that region
EGS_I32 imed
medium index