EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_autoenvelope.cpp
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ auto envelope geometry
5 # Copyright (C) 2016 Randle E. P. Taylor, Rowan M. Thomson,
6 # Marc J. P. Chamberland, D. W. O. Rogers
7 #
8 # This file is part of EGSnrc.
9 #
10 # EGSnrc is free software: you can redistribute it and/or modify it under
11 # the terms of the GNU Affero General Public License as published by the
12 # Free Software Foundation, either version 3 of the License, or (at your
13 # option) any later version.
14 #
15 # EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
16 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 # FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for
18 # more details.
19 #
20 # You should have received a copy of the GNU Affero General Public License
21 # along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
22 #
23 ###############################################################################
24 #
25 # Author: Randle Taylor, 2016
26 #
27 # Contributors: Marc Chamberland
28 # Rowan Thomson
29 # Dave Rogers
30 #
31 ###############################################################################
32 #
33 # egs_autoenvelope was developed for the Carleton Laboratory for
34 # Radiotherapy Physics.
35 #
36 ###############################################################################
37 */
38 
39 
45 #include "egs_input.h"
46 #include "egs_autoenvelope.h"
47 #include "../egs_gtransformed/egs_gtransformed.h"
48 
49 #include <cstdlib>
50 #include <algorithm>
51 #include <numeric>
52 #include <iostream>
53 #include <fstream>
54 
55 #ifdef HAS_SOBOL
56  #include "sobol.h"
57 #endif
58 
59 
60 string EGS_AENVELOPE_LOCAL EGS_AEnvelope::type = "EGS_AEnvelope";
61 string EGS_AENVELOPE_LOCAL EGS_ASwitchedEnvelope::type = "EGS_ASwitchedEnvelope";
62 /* only geometries that support getting regional volume can be used as phantoms */
63 const string EGS_AENVELOPE_LOCAL EGS_AEnvelope::allowed_base_geom_types[] = {"EGS_cSpheres", "EGS_cSphericalShell", "EGS_XYZGeometry", "EGS_RZ"};
64 
65 static char EGS_AENVELOPE_LOCAL geom_class_msg[] = "createGeometry(AEnvelope): %s\n";
66 static char EGS_AENVELOPE_LOCAL base_geom_keyword[] = "base geometry";
67 static char EGS_AENVELOPE_LOCAL inscribed_geom_keyword[] = "inscribed geometry";
68 static char EGS_AENVELOPE_LOCAL inscribed_geom_name_keyword[] = "inscribed geometry name";
69 static char EGS_AENVELOPE_LOCAL transformations_keyword[] = "transformations";
70 static char EGS_AENVELOPE_LOCAL type_keyword[] = "type";
71 static char EGS_AENVELOPE_LOCAL transformation_keyword[] = "transformation";
72 
73 
74 EGS_AEnvelope::EGS_AEnvelope(EGS_BaseGeometry *base_geom,
75  const vector<AEnvelopeAux> inscribed, const string &Name, bool debug, string output_vc_file) :
76  EGS_BaseGeometry(Name), base_geom(base_geom), debug_info(debug), output_vc(output_vc_file) {
77 
78  bool volcor_available = allowedBaseGeomType(base_geom->getType());
79  bool volcor_requested = inscribed.size() > 0 && inscribed[0].vcopts->mode == CORRECT_VOLUME;
80 
81  if (!volcor_available && volcor_requested) {
82 
83  string msg(
84  "EGS_AEnvelope:: Volume correction is not available for geometry type '%s (%s)'. "
85  "Geometry types must implement getVolume. Valid choices are:\n\t"
86  );
87 
88  int end = (int)(sizeof(allowed_base_geom_types)/sizeof(string));
89  for (int i=0; i < end; i++) {
90  msg += allowed_base_geom_types[i] + " ";
91  }
92  msg += "\nPlease set `action = discover -or- correct and zero volume` or use a different base geometry type.\n";
93  egsFatal(msg.c_str(), base_geom->getType().c_str(), base_geom->getName().c_str());
94  }
95 
96 
97  base_geom->ref();
98  nregbase = base_geom->regions();
99  is_convex = base_geom->isConvex();
100  has_rho_scaling = getHasRhoScaling();
101 
102  ninscribed = inscribed.size();
103  if (ninscribed == 0) {
104  egsFatal("EGS_AEnvelope: no inscribed geometries!\n");
105  }
106 
107  // nreg is total regions in geometry = nregbase + sum(nreginscribed_j)
108  nreg = nregbase;
109 
110  // create a transformed copy of the inscribed geometry
111  int global = nreg;
112  for (int geom_idx=0; geom_idx < ninscribed; geom_idx++) {
113 
114  EGS_BaseGeometry *inscribed_geom = inscribed[geom_idx].geom;
115  EGS_AffineTransform *transform = inscribed[geom_idx].transform;
116  VCOptions *vcopt = inscribed[geom_idx].vcopts;
117  EGS_TransformedGeometry *transformed = new EGS_TransformedGeometry(inscribed_geom, *transform);
118 
119  inscribed_geoms.push_back(transformed);
120  transforms.push_back(transform);
121  opts.push_back(vcopt);
122 
123  int ninscribed_reg= transformed->regions();
124  nreg += ninscribed_reg;
125 
126  for (int lreg = 0; lreg < ninscribed_reg; lreg++) {
127  GeomRegPairT local(transformed, lreg);
128  local_to_global_reg[local] = global;
129  global_reg_to_local[global] = local;
130  global++;
131  }
132  }
133 
134  geoms_in_region = new vector<EGS_BaseGeometry *>[nregbase];
135 
136  if (inscribed[0].vcopts->vc_file != "") {
137  vc_results = loadFileResults(inscribed[0].vcopts, base_geom, inscribed_geoms, transforms);
138  egsInformation("loaded from %s\n", inscribed[0].vcopts->vc_file.c_str());
139  }
140  else {
141  // now run volume correction and figure out which regions have inscribed geometries
142  vc_results = findRegionsWithInscribed(inscribed[0].vcopts, base_geom, inscribed_geoms, transforms);
143  }
144 
145 
146  nreg_with_inscribed = 0;
147  for (volcor::RegionGeomSetT::iterator it=vc_results.regions_with_inscribed.begin();
148  it!=vc_results.regions_with_inscribed.end(); it++) {
149  nreg_with_inscribed++;
150  copy(it->second.begin(), it->second.end(), back_inserter(geoms_in_region[it->first]));
151  }
152 
153  if (getNRegWithInscribed() == 0) {
154  egsFatal("EGS_AEnvelope:: Failed to find any regions with inscribed geometries\n");
155  }
156 
157  if (output_vc=="yes" || output_vc=="text" || output_vc=="gzip") {
158  writeVolumeCorrection();
159  }
160 
161  if (debug_info) {
162  printInfo();
163  }
164 
165 }
166 
167 
168 EGS_AEnvelope::~EGS_AEnvelope() {
169  if (!base_geom->deref()) {
170  delete base_geom;
171  }
172 }
173 
174 bool EGS_AEnvelope::allowedBaseGeomType(const string &geom_type) {
175  // Check if the input geometry is one that aenvelope can handle
176 
177  int end = (int)(sizeof(allowed_base_geom_types)/sizeof(string));
178 
179  for (int i=0; i<end; i++) {
180  if (allowed_base_geom_types[i] == geom_type) {
181  return true;
182  }
183  }
184 
185  return false;
186 }
187 
188 int EGS_AEnvelope::getGlobalRegFromLocalReg(EGS_BaseGeometry *g, int local_reg) {
189  return getGlobalRegFromLocal(volcor::GeomRegPairT(g, local_reg));
190 }
191 
192 int EGS_AEnvelope::getGlobalRegFromLocal(const volcor::GeomRegPairT local) const {
193 
194  map<volcor::GeomRegPairT, int>::const_iterator glob_it = local_to_global_reg.find(local);
195  if (glob_it != local_to_global_reg.end()) {
196  return (*glob_it).second;
197  }
198  return -1;
199 }
200 
201 volcor::GeomRegPairT EGS_AEnvelope::getLocalFromGlobalReg(int ireg) const {
202 
203  map<int, volcor::GeomRegPairT>::const_iterator glob_it = global_reg_to_local.find(ireg);
204  if (glob_it != global_reg_to_local.end()) {
205  return (*glob_it).second;
206  }
207  return volcor::GeomRegPairT();
208 }
209 
210 int EGS_AEnvelope::getNRegWithInscribed() const {
211  return nreg_with_inscribed;
212 }
213 
214 bool EGS_AEnvelope::isRealRegion(int ireg) const {
215 
216  bool is_outside = ireg < 0 || ireg >= nreg;
217  if (is_outside) {
218  return false;
219  }
220 
221  bool in_base_geom = ireg < nregbase;
222  if (in_base_geom) {
223  return base_geom->isRealRegion(ireg);
224  }
225 
226  volcor::GeomRegPairT local;
227  local = getLocalFromGlobalReg(ireg);
228  return local.first->isRealRegion(local.second);
229 
230 };
231 
232 
233 bool EGS_AEnvelope::isInside(const EGS_Vector &x) {
234  return base_geom->isInside(x);
235 };
236 
237 
238 int EGS_AEnvelope::isWhere(const EGS_Vector &x) {
239 
240  int base_reg = base_geom->isWhere(x);
241 
242  // which inscribed geometries are in this region
243  vector<EGS_BaseGeometry *> geoms_in_reg = getGeomsInRegion(base_reg);
244 
245  if (base_reg < 0 || geoms_in_reg.size()==0) {
246  return base_reg;
247  }
248 
249  // loop over all inscribed geometries in this region and check
250  // if position is inside any of them
251  for (vector<EGS_BaseGeometry *>::iterator geom=geoms_in_reg.begin(); geom!= geoms_in_reg.end(); ++geom) {
252  int inscribed_reg = (*geom)->isWhere(x);
253  if (inscribed_reg >= 0) {
254  return getGlobalRegFromLocalReg(*geom, inscribed_reg);
255  }
256  }
257 
258  // not in any of the inscribed geometries so must be in envelope region
259  return base_reg;
260 };
261 
262 int EGS_AEnvelope::inside(const EGS_Vector &x) {
263  return isWhere(x);
264 };
265 
266 int EGS_AEnvelope::medium(int ireg) const {
267 
268  if (ireg < nregbase) {
269  return base_geom->medium(ireg);
270  }
271 
272  volcor::GeomRegPairT local = getLocalFromGlobalReg(ireg);
273 
274  return local.first->medium(local.second);
275 };
276 
277 vector<EGS_BaseGeometry *> EGS_AEnvelope::getGeomsInRegion(int ireg) {
278  if (ireg < 0 || ireg >= nregbase) {
279  vector<EGS_BaseGeometry *> empty;
280  return empty;
281  }
282  return geoms_in_region[ireg];
283 };
284 
285 int EGS_AEnvelope::computeIntersections(int ireg, int n, const EGS_Vector &X,
286  const EGS_Vector &u, EGS_GeometryIntersections *isections) {
287 
288  // For a given position x, direction u and region number ireg, this
289  // method returns the number of intersections with the geometry and
290  // distances, medium indeces and relative mass densities for all
291  // intersections, or, if this number is larger than the size n of
292  // isections, -1 (but the n intersections are still put in isections).
293  // If the position is outside, the method checks if the trajectory
294  // intersects the geometry and if yes, puts in the first element of
295  // isections the distance to the entry point and then finds all other
296  // intersections as in the case of x inside.
297 
298  if (n < 1) {
299  return -1;
300  }
301 
302  int isec_idx = 0; // current intersection index
303  EGS_Float t; // distance to next boundary
304  EGS_Float ttot = 0; // total distance along path
305  EGS_Vector x(X); // current position
306 
307  int imed;
308 
309  bool outside_envelope = ireg < 0;
310  bool in_base_geom = ireg >= 0 && ireg < nregbase;
311 
312  if (outside_envelope) {
313  t = 1e30;
314 
315  // region we would hit
316  ireg = howfar(ireg,x,u,t,&imed);
317  bool would_not_intersect = ireg < 0;
318 
319  if (would_not_intersect) {
320  return 0;
321  }
322 
323  isections[0].t = t; // distance to entry point
324  isections[0].rhof = 1;
325  isections[0].ireg = -1;
326  isections[0].imed = -1;
327  ttot = t;
328 
329  x += u*t;
330 
331  isec_idx = 1;
332  }
333  else {
334  // we're already inside the envelope
335  imed = medium(ireg);
336  }
337 
338  // keep looping until we exit the geometry!
339  while (1) {
340  //egsInformation("in loop: j=%d ireg=%d imed=%d x=(%g,%g,%g)\n",
341  // j,ireg,imed,x.x,x.y,x.z);
342  //
343 
344  // why are we setting these here?
345  // on first pass if we were outside envelope then:
346  // ireg is set to entry region
347  // imed is set to med of entry region (set in howfar above)
348  // else
349  // ireg is region of input ireg (current region?)
350  // imed is medium of input ireg (current region?)
351  isections[isec_idx].imed = imed; //
352  isections[isec_idx].ireg = ireg;
353  isections[isec_idx].rhof = getRelativeRho(ireg);
354 
355  if (in_base_geom) { // in one of the regions of the base geometry
356 
357  t = 1e30;
358 
359  // next region will be in base geom by default
360  ireg = base_geom->howfar(ireg,x,u,t,&imed);
361 
362  // if there's inscribed geoms in this region see if we will be intersecting one
363  vector<EGS_BaseGeometry *> geoms_in_reg = getGeomsInRegion(ireg);
364  for (vector<EGS_BaseGeometry *>::iterator geom=geoms_in_reg.begin(); geom!= geoms_in_reg.end(); ++geom) {
365 
366  int inscribed_reg = (*geom)->howfar(-1, x, u, t, &imed);
367 
368  bool hits_inscribed = inscribed_reg >= 0;
369 
370  if (hits_inscribed) {
371  ireg = getGlobalRegFromLocalReg(*geom, inscribed_reg);
372  }
373  }
374 
375  ttot += t;
376  isections[isec_idx++].t = ttot;
377 
378  // left the geometry (implies that no inscribed geoms hit
379  if (ireg < 0) {
380  return isec_idx;
381  }
382 
383 
384  // over limit of number of intersections
385  if (isec_idx >= n) {
386  return -1;
387  }
388 
389  // move to next boundary
390  x += u*t;
391 
392  }
393  else {
394 
395  // in inscribed geometry
396  volcor::GeomRegPairT local = getLocalFromGlobalReg(ireg);
397 
398  int max_isec = n - isec_idx;
399 
400  int ninter_sec = local.first->computeIntersections(local.second, max_isec,
401  x, u, &isections[isec_idx]);
402 
403  // convert intersection region indexes to global indexex
404  int nmax = ninter_sec >= 0 ? ninter_sec + isec_idx : n;
405  for (int i=isec_idx; i < nmax; i++) {
406  isections[i].ireg = getGlobalRegFromLocalReg(local.first, isections[i].ireg);
407  isections[i].t += ttot;
408  }
409 
410  //egsInformation("last intersection: %g\n",isections[nm-1].t);
411  if (ninter_sec < 0) {
412  return ninter_sec;
413  }
414 
415  isec_idx += ninter_sec;
416 
417  if (isec_idx >= n) {
418  return -1;
419  }
420 
421  t = isections[isec_idx-1].t - ttot;
422  x += u*t;
423  ttot = isections[isec_idx-1].t;
424  ireg = base_geom->isWhere(x);
425  //egsInformation("new region: %d\n",ireg);
426 
427  if (ireg < 0) {
428  return isec_idx;
429  }
430 
431  imed = base_geom->medium(ireg);
432  }
433  }
434  return -1;
435 }
436 
437 EGS_Float EGS_AEnvelope::howfarToOutside(int ireg, const EGS_Vector &x, const EGS_Vector &u) {
438 
439  if (ireg < 0) {
440  return 0;
441  }
442 
443  EGS_Float d;
444 
445  bool in_base_geom = ireg < nregbase;
446  if (in_base_geom) {
447  d = base_geom->howfarToOutside(ireg, x, u);
448  }
449  else if (base_geom->regions() == 1) {
450  d = base_geom->howfarToOutside(0, x, u);
451  }
452  else {
453  int ir = base_geom->isWhere(x);
454  d = base_geom->howfarToOutside(ir,x,u);
455  }
456  return d;
457 };
458 
459 
460 int EGS_AEnvelope::howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u,
461  EGS_Float &t, int *newmed, EGS_Vector *normal) {
462 
463  bool inside_geom = ireg >= 0;
464  if (inside_geom) {
465 
466  bool inside_base_geom = ireg < nregbase;
467 
468  if (inside_base_geom) {
469 
470  int base_reg = base_geom->howfar(ireg,x,u,t,newmed,normal);
471  vector<EGS_BaseGeometry *> geoms_in_reg = getGeomsInRegion(ireg);
472 
473  bool no_inscribed_in_reg = geoms_in_reg.size() == 0;
474  if (no_inscribed_in_reg) {
475  return base_reg;
476  }
477 
478  int inscribed_global = -1;
479  for (vector<EGS_BaseGeometry *>::iterator geom=geoms_in_reg.begin(); geom!= geoms_in_reg.end(); ++geom) {
480  int local_reg = (*geom)->howfar(-1, x, u, t, newmed, normal);
481  bool hits_inscribed = local_reg >= 0;
482  if (hits_inscribed) {
483  inscribed_global = getGlobalRegFromLocalReg(*geom, local_reg);
484  }
485  }
486 
487  bool hit_inscribed_first = inscribed_global >= 0;
488 
489  if (hit_inscribed_first) {
490  return inscribed_global;
491  }
492 
493  return base_reg;
494  }
495  else {
496 
497  // if here, we are in an inscribed geometry.
498  volcor::GeomRegPairT local = getLocalFromGlobalReg(ireg);
499 
500  // and then check if we will hit a boundary in this geometry.
501  int new_local = local.first->howfar(local.second, x, u, t, newmed, normal);
502  bool hit_boundary_in_inscribed = new_local >=0;
503  if (hit_boundary_in_inscribed) {
504  return getGlobalRegFromLocalReg(local.first, new_local);
505  }
506 
507  // new_local < 0 implies that we have exited the inscribed geometry
508  // => check to see in which base geometry region we are.
509  int inew = base_geom->isWhere(x + u*t);
510  if (inew >= 0 && newmed) {
511  *newmed = base_geom->medium(inew);
512  }
513  return inew;
514  }
515  }
516 
517  // if here, we are outside the base geometry.
518  // check to see if we will enter.
519  int new_base_reg = base_geom->howfar(ireg,x,u,t,newmed,normal);
520  vector<EGS_BaseGeometry *> geoms_in_reg = getGeomsInRegion(new_base_reg);
521  bool has_geoms_in_reg = geoms_in_reg.size() > 0;
522 
523  //cout << "new reg "<<new_base_reg << " has geoms " << has_geoms_in_reg << " t "<<t<<endl;
524  if (new_base_reg >= 0 && has_geoms_in_reg) {
525 
526  // check if we will be inside an inscribed geometry when we enter
527  for (vector<EGS_BaseGeometry *>::iterator geom=geoms_in_reg.begin(); geom!= geoms_in_reg.end(); ++geom) {
528  int new_local_reg = (*geom)->isWhere(x+u*t);
529  bool landed_in_inscribed = new_local_reg >= 0;
530  if (landed_in_inscribed) {
531  if (newmed) {
532  *newmed = (*geom)->medium(new_local_reg);
533  }
534  return getGlobalRegFromLocalReg(*geom, new_local_reg);
535  }
536 
537  }
538  }
539 
540  return new_base_reg;
541 };
542 
543 
544 EGS_Float EGS_AEnvelope::hownear(int ireg, const EGS_Vector &x) {
545  bool inside_envelope = ireg >= 0;
546 
547  if (!inside_envelope) {
548  return base_geom->hownear(ireg, x);
549  }
550 
551  bool in_base_geom = ireg < nregbase ;
552  if (!in_base_geom) {
553  volcor::GeomRegPairT local = getLocalFromGlobalReg(ireg);
554  return local.first->hownear(local.second, x);
555  }
556 
557  EGS_Float tmin = base_geom->hownear(ireg, x);
558  if (tmin <= 0) {
559  return tmin;
560  }
561  vector<EGS_BaseGeometry *> geoms_in_reg = getGeomsInRegion(ireg);
562  for (vector<EGS_BaseGeometry *>::iterator geom=geoms_in_reg.begin(); geom!= geoms_in_reg.end(); ++geom) {
563  EGS_Float local_t = (*geom)->hownear(-1, x);
564  if (local_t < tmin) {
565  tmin = local_t;
566  if (tmin < 0) {
567  return tmin;
568  }
569  }
570  }
571  return tmin;
572 };
573 
574 bool EGS_AEnvelope::hasBooleanProperty(int ireg, EGS_BPType prop) const {
575  if (ireg >= 0 && ireg < nreg) {
576  if (ireg < nregbase) {
577  return base_geom->hasBooleanProperty(ireg,prop);
578  }
579  volcor::GeomRegPairT local = getLocalFromGlobalReg(ireg);
580  return local.first->hasBooleanProperty(local.second, prop);
581  }
582  return false;
583 };
584 
585 void EGS_AEnvelope::setBooleanProperty(EGS_BPType prop) {
586  setPropertyError("setBooleanProperty()");
587 };
588 
589 void EGS_AEnvelope::addBooleanProperty(int bit) {
590  setPropertyError("addBooleanProperty()");
591 };
592 
593 void EGS_AEnvelope::setBooleanProperty(EGS_BPType prop, int start, int end,int step) {
594  setPropertyError("setBooleanProperty()");
595 };
596 
597 void EGS_AEnvelope::addBooleanProperty(int bit,int start,int end,int step) {
598  setPropertyError("addBooleanProperty()");
599 };
600 
601 int EGS_AEnvelope::getMaxStep() const {
602  int nstep = base_geom->getMaxStep();
603  for (size_t j=0; j< inscribed_geoms.size(); ++j) {
604  nstep += inscribed_geoms[j]->getMaxStep();
605  }
606  return nstep + inscribed_geoms.size();
607 };
608 
609 EGS_Float EGS_AEnvelope::getVolume(int ireg) {
610 
611  if (ireg < 0) {
612  return -1;
613  }
614  else if (ireg < nregbase) {
615  return vc_results.corrected_volumes[ireg];
616  }
617 
618  volcor::GeomRegPairT local = global_reg_to_local[ireg];
619 
620  return local.first->getVolume(local.second);
621 };
622 
623 EGS_Float EGS_AEnvelope::getCorrectionRatio(int ireg) {
624 
625  if (0 <= ireg && ireg < nregbase) {
626  return vc_results.corrected_volumes[ireg]/vc_results.uncorrected_volumes[ireg];
627  }
628 
629  return 1;
630 };
631 
632 
633 /* Print information about the geometry. If `print debug info = yes` is present
634  * in the geometry, extra information about which regions of the base geometry
635  * contain inscribed geometries and how their volumes were corrected.
636 */
637 void EGS_AEnvelope::printInfo() const {
638 
640 
641  if (debug_info) {
642 
643  for (int ir=0; ir < nregbase; ir++) {
644  if (fabs(vc_results.uncorrected_volumes[ir] - vc_results.corrected_volumes[ir]) > 1E-8) {
645  egsInformation(" volume of region %d was corrected from %.5E g to %.5E g\n",
646  ir, vc_results.uncorrected_volumes[ir], vc_results.corrected_volumes[ir]);
647  }
648  }
649 
650  for (int ir=0; ir < nregbase; ir++) {
651  if (geoms_in_region[ir].size() > 0) {
652  egsInformation(" region %d has %d inscribed geometries\n", ir, geoms_in_region[ir].size());
653  }
654  }
655  }
656 
657  vc_results.outputResults();
658 
659  egsInformation(" base geometry = %s (type %s) nreg=%d\n",base_geom->getName().c_str(),
660  base_geom->getType().c_str(), nregbase, base_geom->regions());
661  egsInformation(" # of inscribed geometries= %d in %d regions\n",inscribed_geoms.size(), getNRegWithInscribed());
662 
663 
664  egsInformation("=======================================================\n");
665 }
666 
667 void EGS_AEnvelope::writeVCToFile(ostream &out) {
668 
669  vector<int> to_write;
670 
671  for (int i=0; i < base_geom->regions(); i++) {
672  EGS_Float cor = vc_results.corrected_volumes[i];
673  EGS_Float uncor = vc_results.uncorrected_volumes[i];
674  bool has_correction = fabs(cor-uncor) > 1E-8;
675  if (has_correction) {
676  to_write.push_back(i);
677  }
678  }
679 
680  size_t nrecords = to_write.size();
681  out << nrecords << endl;
682 
683  for (size_t i=0; i<to_write.size(); i++) {
684  int ir = to_write[i];
685  EGS_Float cor = vc_results.corrected_volumes[ir];
686 
687  // find out the indexes of geometries inscribed in this region (if any)
688  vector<EGS_BaseGeometry *> geoms_in_reg = getGeomsInRegion(ir);
689  vector<int> geom_idxs_in_reg;
690  for (size_t i=0; i<geoms_in_reg.size(); i++) {
691  size_t pos = find(inscribed_geoms.begin(), inscribed_geoms.end(), geoms_in_reg[i]) - inscribed_geoms.begin();
692  if (pos < inscribed_geoms.size()) {
693  geom_idxs_in_reg.push_back(pos);
694  }
695  }
696 
697  int ninscribed = (int)geom_idxs_in_reg.size();
698 
699  out << ir << " " << cor;
700  // write number of and indexes of geometries inscribed in this region
701  out << " " << ninscribed;
702  for (size_t gidx = 0; gidx < geom_idxs_in_reg.size(); gidx++) {
703  out << " " <<geom_idxs_in_reg[gidx];
704  }
705  out << endl;
706  }
707 
708 }
709 
710 void EGS_AEnvelope::writeVolumeCorrection() {
711 
712  egsInformation("\nVolume Correction Output File \n%s\n", string(80,'=').c_str());
713 
714  string gname = base_geom->getName();
715  string fname = gname+".autoenv.volcor";
716  fname += (output_vc == "gzip" ? ".gz" : "");
717  string mode = (output_vc == "gzip" ? "gzip" : "text");
718 
720  "Writing volume correction file for %s to %s file %s\n",
721  gname.c_str(), mode.c_str(), fname.c_str());
722 
723  if (output_vc == "gzip") {
724 #ifdef HAS_GZSTREAM
725  ogzstream outgz(fname.c_str());
726  writeVCToFile(outgz);
727  outgz.close();
728 #endif
729  }
730  else {
731  ofstream out(fname.c_str());
732  writeVCToFile(out);
733  out.close();
734  }
735 
736 }
737 
738 
739 /* check if the base or any of inscribed geometries have density scaling on */
740 bool EGS_AEnvelope::getHasRhoScaling() {
741 
742  has_rho_scaling = base_geom->hasRhoScaling();
743  if (has_rho_scaling) {
744  return true;
745  }
746 
747  for (size_t geom_idx=0; geom_idx < inscribed_geoms.size(); geom_idx++) {
748  if (inscribed_geoms[geom_idx]->hasRhoScaling()) {
749  return true;
750  }
751  }
752  return false;
753 
754 }
755 
756 
757 void EGS_AEnvelope::setMedia(EGS_Input *,int,const int *) {
758  egsWarning("EGS_AEnvelope::setMedia: don't use this method. Use the\n"
759  " setMedia() methods of the geometry objects that make up this geometry\n");
760 }
761 
762 
763 void EGS_AEnvelope::setRelativeRho(int start, int end, EGS_Float rho) {
764  setRelativeRho(0);
765 }
766 
767 
768 void EGS_AEnvelope::setRelativeRho(EGS_Input *) {
769  egsWarning("EGS_AEnvelope::setRelativeRho(): don't use this method."
770  " Use the\n setRelativeRho methods of the geometry objects that make up"
771  " this geometry\n");
772 }
773 
774 
775 EGS_Float EGS_AEnvelope::getRelativeRho(int ireg) const {
776  if (ireg < 0 || ireg >= nreg) {
777  return 1;
778  }
779  if (ireg < nregbase) {
780  EGS_Float v = base_geom->getRelativeRho(ireg);
781  return v;
782  }
783  volcor::GeomRegPairT local = getLocalFromGlobalReg(ireg);
784  if (local.first) {
785  return local.first->getRelativeRho(local.second);
786  }
787  return 1;
788 };
789 
790 
791 /*************************************************************************/
792 /* Switched envelope *****************************************************/
793 /*************************************************************************/
794 
795 EGS_ASwitchedEnvelope::EGS_ASwitchedEnvelope(EGS_BaseGeometry *base_geom,
796  const vector<AEnvelopeAux> inscribed, const string &Name, bool debug, string output_vc_file):
797  EGS_AEnvelope(base_geom, inscribed, Name, debug, output_vc_file) {
798 
799  // activate a single geometry
800  cur_ptr = 0;
801  active_inscribed.push_back(inscribed_geoms[cur_ptr]);
802 
803 };
804 
805 
806 //TODO: this gets called a lot and is probably quite slow. Instead fo doing a
807 //set intersection on every call we can probably do it once when activated
808 //geometries change and cache it
809 vector<EGS_BaseGeometry *> EGS_ASwitchedEnvelope::getGeomsInRegion(int ireg) {
810  // find which geometries are present AND active in this region
811  vector<EGS_BaseGeometry *> active_in_reg;
812  if ((0 <= ireg) && (ireg < nregbase)) {
813  set_intersection(
814  geoms_in_region[ireg].begin(), geoms_in_region[ireg].end(),
815  active_inscribed.begin(), active_inscribed.end(),
816  back_inserter(active_in_reg)
817  );
818  }
819 
820  return active_in_reg;
821 
822 }
823 
824 
825 void EGS_ASwitchedEnvelope::setActiveGeometries(vector<EGS_BaseGeometry *> geoms) {
826  active_inscribed.clear();
827  active_inscribed = geoms;
828  cur_ptr = -1;
829  if (geoms.size() <= 0) {
830  return;
831  }
832 
833  for (size_t i=0; i < inscribed_geoms.size(); i++) {
834  if (inscribed_geoms[i] == geoms[0]) {
835  cur_ptr = i;
836  return;
837  }
838  }
839 }
840 
841 
842 void EGS_ASwitchedEnvelope::setActiveGeometries(vector<int> geom_indexes) {
843 
844  active_inscribed.clear();
845 
846  for (size_t i = 0; i < geom_indexes.size(); i++) {
847  int idx = geom_indexes[i];
848  if (idx < 0 || idx >= ninscribed) {
849  egsFatal("EGS_ASwitchedEnvelope:: %d is not a valid geometry index\n", idx);
850  }
851  active_inscribed.push_back(inscribed_geoms[i]);
852  }
853 
854  cur_ptr = geom_indexes[0];
855 
856 }
857 
858 
860 
861  if ((0 <= ireg) && (ireg < nregbase)) {
862  size_t nactive = getGeomsInRegion(ireg).size();
863  return nactive > 0;
864  }
865 
866  return false;
867 }
868 
869 
871 
872  if ((0 <= ireg) && (ireg < nregbase)) {
873  size_t nactive = getGeomsInRegion(ireg).size();
874  size_t ntotal = EGS_AEnvelope::getGeomsInRegion(ireg).size();
875  return nactive < ntotal;
876  }
877 
878  return false;
879 
880 }
881 
882 void EGS_ASwitchedEnvelope::setActiveByIndex(int inscribed_index) {
883  vector<EGS_BaseGeometry *> to_activate;
884  to_activate.push_back(inscribed_geoms[inscribed_index]);
885  setActiveGeometries(to_activate);
886  cur_ptr = inscribed_index;
887 }
888 
889 void EGS_ASwitchedEnvelope::activateByIndex(int inscribed_index) {
890  active_inscribed.push_back(inscribed_geoms[inscribed_index]);
891 }
892 
893 void EGS_ASwitchedEnvelope::deactivateByIndex(int inscribed_index) {
894  vector<EGS_BaseGeometry *>::iterator loc = find(
895  active_inscribed.begin(), active_inscribed.end(),
896  inscribed_geoms[inscribed_index]
897  );
898  if (loc != active_inscribed.end()) {
899  active_inscribed.erase(loc);
900  }
901 }
902 
904  cur_ptr = cur_ptr == ninscribed -1 ? 0 : cur_ptr + 1;
905  vector<EGS_BaseGeometry *> to_activate;
906  to_activate.push_back(inscribed_geoms[cur_ptr]);
907  setActiveGeometries(to_activate);
908 }
909 
910 
911 vector<EGS_AffineTransform *> EGS_AEnvelope::createTransforms(EGS_Input *input) {
912 
913  vector<EGS_AffineTransform *> transforms;
914  if (input) {
915  EGS_Input *i;
916 
917  while ((i = input->takeInputItem(transformation_keyword))) {
919  if (!transform) {
920  egsWarning("Invalid transform input given\n");
921 
922  }
923  transforms.push_back(transform);
924  delete i;
925 
926  }
927  }
928 
929  return transforms;
930 
931 }
932 
933 
934 extern "C" {
935 
936  EGS_AENVELOPE_EXPORT EGS_BaseGeometry *createGeometry(EGS_Input *input) {
937 
938  if (!input) {
939  egsWarning(geom_class_msg, "null input");
940  return 0;
941  }
942 
943  bool debug;
944  vector<string> debug_choices;
945  debug_choices.push_back("no");
946  debug_choices.push_back("yes");
947  debug = input->getInput("print debug info", debug_choices, 0);
948 
949  int output_vc_file_choice;
950  vector<string> vc_file_choices;
951  vc_file_choices.push_back("no");
952  vc_file_choices.push_back("yes");
953  vc_file_choices.push_back("text");
954  vc_file_choices.push_back("gzip");
955  output_vc_file_choice = input->getInput("output volume correction file", vc_file_choices, 0);
956  string output_vc_file = vc_file_choices[output_vc_file_choice];
957 
958 #ifndef HAS_GZSTREAM
959  if (output_vc_file == "gzip") {
960  egsWarning(
961  "GZip file output requested but not compiled with gzstream.\n"
962  "Please recompile with gzstream support.\n"
963  );
964  return 0;
965  }
966 #endif
967 
968 
969 
970  string base_geom_name;
971  int err = input->getInput(base_geom_keyword, base_geom_name);
972  if (err) {
973  egsWarning(geom_class_msg, ("'"+string(base_geom_keyword)+"' input not found").c_str());
974  return 0;
975  }
976 
977  string type;
978  err = input->getInput(type_keyword, type);
979  if (err) {
980  type = "AEnvelope";
981  }
982 
983  EGS_BaseGeometry *base_geom = EGS_BaseGeometry::getGeometry(base_geom_name);
984  if (!base_geom) {
985  egsWarning(geom_class_msg, ("Unable to find geometry '"+base_geom_name+"'").c_str());
986  return 0;
987  }
988 
989  EGS_Input *inscribed_input = input->takeInputItem(inscribed_geom_keyword);
990  if (!inscribed_input) {
991  egsWarning(geom_class_msg, ("Missing '"+string(inscribed_geom_keyword)+"' input item").c_str());
992  return 0;
993  }
994 
995  string inscribed_geom_name;
996  err = inscribed_input->getInput(inscribed_geom_name_keyword, inscribed_geom_name);
997  if (err) {
998  egsWarning(geom_class_msg, ("'"+string(inscribed_geom_name_keyword)+"' input not found").c_str());
999  return 0;
1000  }
1001 
1002  EGS_BaseGeometry *inscribed_geom = EGS_BaseGeometry::getGeometry(inscribed_geom_name);
1003  if (!inscribed_geom) {
1004  egsWarning(geom_class_msg, ("Unable to find geometry '"+inscribed_geom_name+"'").c_str());
1005  return 0;
1006  }
1007 
1008  vector<EGS_AffineTransform *> transforms;
1009  EGS_Input *trans_input = inscribed_input->takeInputItem(transformations_keyword);
1010  if (trans_input) {
1011  transforms = EGS_AEnvelope::createTransforms(trans_input);
1012  }
1013 
1014  delete trans_input;
1015 
1016  EGS_Input *volcor_input = inscribed_input->takeInputItem("region discovery");
1017  VCOptions *vcopts = new VCOptions(volcor_input);
1018  if (!vcopts->valid) {
1019  egsWarning(geom_class_msg, "Missing or invalid 'region discovery' input item");
1020  return 0;
1021  }
1022 
1023  delete inscribed_input;
1024 
1025 
1026  vector<AEnvelopeAux> inscribed;
1027  if (transforms.size()>0) {
1028  for (size_t i=0; i < transforms.size(); i++) {
1029  inscribed.push_back(AEnvelopeAux(inscribed_geom, transforms[i], vcopts));
1030  }
1031  }
1032  else {
1033  EGS_AffineTransform *unityt = new EGS_AffineTransform();
1034  inscribed.push_back(AEnvelopeAux(inscribed_geom, unityt, vcopts));
1035  }
1036 
1037  EGS_BaseGeometry *result;
1038  if (type == "EGS_ASwitchedEnvelope") {
1039  result = new EGS_ASwitchedEnvelope(base_geom, inscribed, "", (bool)debug, output_vc_file);
1040  }
1041  else {
1042  result = new EGS_AEnvelope(base_geom, inscribed, "", (bool)debug, output_vc_file);
1043  }
1044 
1045  result->setName(input);
1046  return result;
1047  }
1048 
1049 }
virtual const string & getType() const =0
Get the geometry type.
virtual bool isInside(const EGS_Vector &x)=0
Is the position x inside the geometry?
virtual bool hasBooleanProperty(int ireg, EGS_BPType prop) const
Is the boolean property prop set for region ireg ?
A fast envelope geometry with automatic region detection.
int regions() const
Returns the number of local regions in this geometry.
int deref()
Decrease the reference count to this geometry.
An envelope geometry with automatic inscribed region detection (inspired by EGS_FastEnvelope) ...
A helper class for initializing auto envelopes.
virtual int medium(int ireg) const
Returns the medium index in region ireg.
virtual void printInfo() const
Print information about this geometry.
EGS_Float t
distance to next region boundary
This geometry type allows you to activate and deactivate inscribed geometries in custom egspp user co...
void outputResults() const
Definition: volcor.h:475
virtual EGS_Float howfarToOutside(int ireg, const EGS_Vector &x, const EGS_Vector &u)
EGS_Float rhof
relative mass density in that region
A class providing affine transformations.
EGS_Input class header file.
virtual int getMaxStep() const
Returns the maximum number of steps through the geometry.
A class representing 3D vectors.
Definition: egs_vector.h:56
pair< EGS_BaseGeometry *, int > GeomRegPairT
Definition: volcor.h:88
bool hasRhoScaling() const
Does this geometry object have a mass density scaling feature?
int nregbase
Number of regions in the base geometry.
virtual int isWhere(const EGS_Vector &x)=0
In which region is poisition x?
vector< EGS_BaseGeometry * > inscribed_geoms
The inscribed geometries.
void activateByIndex(int inscribed_index)
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.
EGS_I32 imed
medium index
EGS_GLIB_EXPORT EGS_BaseGeometry * createGeometry(EGS_Input *input)
Definition: egs_glib.cpp:57
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
vector< EGS_AffineTransform * > transforms
The inscribed geometries.
bool hasActiveGeom(int ireg)
VCResults findRegionsWithInscribed(VCOptions *opts, EGS_BaseGeometry *base, vector< EGS_BaseGeometry * > inscribed, vector< EGS_AffineTransform * > transforms)
Run the MC simulation.
Definition: volcor.h:557
bool isConvex() const
Is the geometry convex?
void setActiveGeometries(vector< EGS_BaseGeometry * > geoms)
vector< EGS_Float > uncorrected_volumes
Definition: volcor.h:439
A transformed geometry.
static EGS_BaseGeometry * getGeometry(const string &Name)
Get a pointer to the geometry named Name.
VCResults loadFileResults(VCOptions *opts, EGS_BaseGeometry *base, vector< EGS_BaseGeometry * > inscribed, vector< EGS_AffineTransform * > transforms)
Load volume corrections from external file.
Definition: volcor.h:677
EGS_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
static bool allowedBaseGeomType(const string &geom_type)
function for checking whether a given geometry type is allowed to be used as a base geometry ...
EGS_InfoFunction EGS_EXPORT egsInformation
Always use this function for reporting the progress of a simulation and any other type of information...
EGS_BaseGeometry * base_geom
The envelope geometry.
void setName(EGS_Input *inp)
Set the name of the geometry from the input inp.
vector< EGS_Float > corrected_volumes
Definition: volcor.h:440
bool hasInactiveGeom(int ireg)
virtual EGS_Float hownear(int ireg, const EGS_Vector &x)=0
Calculate the distance to a boundary for position x in any direction.
virtual EGS_Float getRelativeRho(int ireg) const
Get the relative mass density in region ireg.
void deactivateByIndex(int inscribed_index)
void transform(EGS_Vector &v) const
Transforms the vector v.
static EGS_AffineTransform * getTransformation(EGS_Input *inp)
Constructs an affine transformation object from the input pointed to by inp and returns a pointer to ...
bool has_rho_scaling
Does this geometry have relative mass density scvaling?
int ref()
Increase the reference count to this geometry.
virtual bool isRealRegion(int ireg) const
Returnes true if ireg is a real region, false otherwise.
static string type
Geometry type.
EGS_I32 ireg
region index
Volume correction initialization helper class.
Definition: volcor.h:241
A class for storing information in a tree-like structure of key-value pairs. This class is used throu...
Definition: egs_input.h:182
const string & getName() const
Get the name of this geometry.
int nreg
Number of local regions in this geometry.
void setMedia(EGS_Input *, int, const int *)
Don&#39;t set media for an envelope geometry.
EGS_Input * takeInputItem(const string &key, bool self=true)
Get the property named key.
Definition: egs_input.cpp:226
static vector< EGS_AffineTransform * > createTransforms(EGS_Input *inpt)
Take a block of transformations and return vector of EGS_AffineTransforms.
int ninscribed
Number of regions in the base geometry.
void setActiveByIndex(int inscribed_index)
static string type
Geometry type.
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
EGS_InfoFunction EGS_EXPORT egsWarning
Always use this function for reporting warnings.