EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_lattice.h
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ lattice geometry headers
5 # Copyright (C) 2019 Rowan Thomson and Martin Martinov
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: Martin Martinov, 2019
25 #
26 # Contributors:
27 #
28 ###############################################################################
29 #
30 # When egs_lattice is used for publications, please cite the following paper:
31 #
32 # Martinov, Martin P., and Rowan M. Thomson. Taking EGSnrc to new lows:
33 # Development of egs++ lattice geometry and testing with microscopic
34 # geometries. Medical Physics 47, 3225-3232 (2020).
35 #
36 ###############################################################################
37 */
38 
39 
44 #include "egs_base_geometry.h"
45 #include "../egs_gtransformed/egs_gtransformed.h"
46 
47 #ifndef EGS_LATTICE_GEOMETRY_
48 #define EGS_LATTICE_GEOMETRY_
49 
50 #ifdef WIN32
51 
52  #ifdef BUILD_LATTICE_DLL
53  #define EGS_LATTICE_EXPORT __declspec(dllexport)
54  #else
55  #define EGS_LATTICE_EXPORT __declspec(dllimport)
56  #endif
57  #define EGS_LATTICE_LOCAL
58 
59 #else
60 
61  #ifdef HAVE_VISIBILITY
62  #define EGS_LATTICE_EXPORT __attribute__ ((visibility ("default")))
63  #define EGS_LATTICE_LOCAL __attribute__ ((visibility ("hidden")))
64  #else
65  #define EGS_LATTICE_EXPORT
66  #define EGS_LATTICE_LOCAL
67  #endif
68 
69 #endif
70 
145 class EGS_LATTICE_EXPORT EGS_Lattice : public EGS_BaseGeometry {
146 protected:
147 
150  int ind;
151  int maxStep;
152  EGS_Float a, b, c;
153  string type;
154 
155 public:
156 
157  EGS_Lattice(EGS_BaseGeometry *B, EGS_BaseGeometry *S, int i, EGS_Float x,
158  EGS_Float y, EGS_Float z, const string &Name = "");
159  ~EGS_Lattice();
160 
161  EGS_Vector closestPoint(const EGS_Vector &x) {
162  return EGS_Vector(int(round(x.x/a))*a,
163  int(round(x.y/b))*b,
164  int(round(x.z/c))*c);
165  };
166 
167  int computeIntersections(int ireg, int n, const EGS_Vector &x,
168  const EGS_Vector &u, EGS_GeometryIntersections *isections) {
169  return base->computeIntersections(ireg,n,x,u,isections);
170  };
171 
172  bool isRealRegion(int ireg) const {
173  if (ireg < base->regions()) { // If ireg is less than base regions
174  return base->isRealRegion(ireg); // check against base regions
175  }
176  return sub->isRealRegion(ireg - base->regions()); // then check sub regions
177  };
178 
179  bool isInside(const EGS_Vector &x) {
180  return base->isInside(x);
181  };
182 
183  int isWhere(const EGS_Vector &x) {
184  int temp = base->isWhere(x);
185  if (temp == ind) {// Are we in the subgeom?
186  sub->setTransformation(closestPoint(x));
187  if (sub->isInside(x)) {
188  return sub->isWhere(x) + base->regions();
189  }
190  }
191  return temp; // otherwise base geom
192  };
193 
194  int inside(const EGS_Vector &x) {
195  return isWhere(x);
196  };
197 
198  int medium(int ireg) const {
199  if (ireg >= base->regions()) { // If ireg is greater than base regions
200  return sub->medium(ireg-base->regions());
201  }
202  return base->medium(ireg); // If in base region
203  };
204 
205  EGS_Float howfarToOutside(int ireg, const EGS_Vector &x, const EGS_Vector &u) {
206  if (ireg < 0) { // Error catch, not inside
207  return 0;
208  }
209 
210  // Are we in the subgeom
211  if (ireg >= base->regions()) {
212  return base->howfarToOutside(ind,x,u);
213  }
214  return base->howfarToOutside(ireg,x,u);
215  };
216 
217  int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u,
218  EGS_Float &t,int *newmed=0, EGS_Vector *normal=0) {
219  // Catch t=0 exception, which sometimes causes issues when ind is the region neighboring -1
220  if (t < epsilon) {
221  t = 0.0;
222  return ireg;
223  }
224 
225  if (ireg >= base->regions()) { // Are we in the subgeom?
226  sub->setTransformation(closestPoint(x));
227 
228  // Do the howfar call
229  EGS_Float tempT = t;
230  int tempReg = sub->howfar(ireg-base->regions(),x,u,tempT,newmed,normal);
231 
232  // If we do leave the subgeom, we want to return the index
233  // of the region of the base geometry we would be entering,
234  // which is not necessarily ind if subgeom is at a boundary
235  EGS_Vector newX = u;
236  newX.normalize();
237  newX = x + (newX * tempT);
238 
239  if (base->isWhere(newX) != ind) {// If we leave ind region of base
240  t = tempT;
241 
242  tempReg = base->howfar(ind,x,u,t,0,normal);
243  if (newmed && tempReg >= 0) {
244  *newmed = base->medium(tempReg);
245  }
246 
247  //if (t<0) {
248  // egsWarning("Returning negative t from subgeom of %s\n",getName().c_str());
249  // egsWarning("howfar(%7d,[%e,%e,%e],[%e,%e,%e],%e)\n",ireg, x.x, x.y, x.z, u.x, u.y, u.z, t);
250  //}
251  return tempReg;
252  }
253 
254  if (!(tempReg+1)) {// If we leave sub geom back into ind
255  if (newmed) {
256  *newmed = base->medium(ind);
257  }
258  t = tempT;
259  return ind;
260  }
261 
262  // else we stay in sub geom
263  t = tempT;
264  //if (t<0) {
265  // egsWarning("Returning negative t from subgeom of %s\n",getName().c_str());
266  // egsWarning("howfar(%7d,[%e,%e,%e],[%e,%e,%e],%e)\n",ireg, x.x, x.y, x.z, u.x, u.y, u.z, t);
267  //}
268  return tempReg+base->regions();
269  }
270  else if (ireg == ind) {// If we are in the region that could contain subgeoms
271  // Determine the path travelled ------------------------------------------------ //
272  EGS_Float tempT = t;
273  base->howfar(ireg,x,u,tempT); // Get how far it is in temp
274  EGS_Float max = tempT;
275 
276  // Iterate through the lattice ------------------------------------------------- //
277  EGS_Float min = max, minX = max, minY = max, minZ = max; // Distance to closest subgeom for three different cases
278  EGS_Vector unit = u;
279  unit.normalize();
280  EGS_Vector xInt = unit*(a/4.0/fabs(unit.x)), yInt = unit*(b/4.0/fabs(unit.y)), zInt = unit*(c/4.0/fabs(unit.z));
281  EGS_Vector tempP; // Temporarily hold indices of subgeom we are testing against
282  EGS_Vector finalP, finalX, finalY, finalZ; // Current closest intersecting subgeom
283 
284  EGS_Vector x0;
285  EGS_Float max2 = max*max;
286 
287  x0 = x;
288  while ((x-x0).length2() < max2) {// Quit as soon as our testing position is beyond current closest intersection (which could be tempT)
289  tempT = minX;
290  tempP = closestPoint(x0);
291  sub->setTransformation(tempP);
292  if (sub->howfar(-1,x,u,tempT)+1) // Intersection!
293  if (tempT < max && tempT > 0) {
294  finalX = tempP;
295  minX = tempT;
296  break;
297  }
298  x0 += xInt;
299  }
300 
301  x0 = x;
302  while ((x-x0).length2() < max2) {// Quit as soon as our testing position is beyond current closest intersection (which could be tempT)
303  tempT = minY;
304  tempP = closestPoint(x0);
305  sub->setTransformation(tempP);
306  if (sub->howfar(-1,x,u,tempT)+1) // Intersection!
307  if (tempT < max && tempT > 0) {
308  finalY = tempP;
309  minY = tempT;
310  break;
311  }
312  x0 += yInt;
313  }
314 
315  x0 = x;
316  while ((x-x0).length2() < max2) {// Quit as soon as our testing position is beyond current closest intersection (which could be tempT)
317  tempT = minZ;
318  tempP = closestPoint(x0);
319  sub->setTransformation(tempP);
320  if (sub->howfar(-1,x,u,tempT)+1) // Intersection!
321  if (tempT < max && tempT > 0) {
322  finalZ = tempP;
323  minZ = tempT;
324  break;
325  }
326  x0 += zInt;
327  }
328 
329  finalP = finalX;
330  min = minX;
331  if (min > minY) {
332  min = minY;
333  finalP = finalY;
334  }
335  if (min > minZ) {
336  min = minZ;
337  finalP = finalZ;
338  }
339 
340  // Check last point
341  tempT = max;
342  sub->setTransformation(closestPoint(x+unit*tempT));
343  sub->howfar(-1,x,u,tempT);
344  if (tempT < min && tempT > 0) {
345  min = tempT;
346  finalP = closestPoint(x+unit*tempT);
347  }
348 
349  // We did intersect subgeom
350  if (min < max) {
351  tempT = t;
352  sub->setTransformation(finalP);
353  int tempReg = sub->howfar(-1,x,u,tempT,newmed,normal);
354  if (tempReg+1) {
355  if (newmed && tempReg >= 0) {
356  *newmed = sub->medium(tempReg);
357  }
358  t = tempT;
359 
360  //if (t<0) {
361  // egsWarning("Returning negative t from region %d in base geom of %s\n", ind, getName().c_str());
362  // egsWarning("howfar(%7d,[%e,%e,%e],[%e,%e,%e],%e)\n",ireg, x.x, x.y, x.z, u.x, u.y, u.z, t);
363  //}
364  return tempReg+base->regions();
365  }
366  }
367 
368  // We didn't intersect subgeom
369  int tempReg = base->howfar(ireg,x,u,t,newmed,normal);
370  //if (t<0) {
371  // egsWarning("Returning negative t from region %d in base geom of %s\n", ind, getName().c_str());
372  // egsWarning("howfar(%7d,[%e,%e,%e],[%e,%e,%e],%e)\n",ireg, x.x, x.y, x.z, u.x, u.y, u.z, t);
373  //}
374  if (newmed && tempReg >= 0) {
375  *newmed = base->medium(tempReg);
376  }
377 
378  return tempReg;
379  }
380  else {// Not in region containing subgeoms, then it is quite easy
381  // unless we enter directly into a subgeom when entering
382  // region ind
383  int tempReg = base->howfar(ireg,x,u,t,newmed,normal);
384 
385  // If we enter region ind
386  if (tempReg == ind) {
387  // newX is the point of entrance
388  EGS_Vector newX = u;
389  newX.normalize();
390  newX = x + (newX * t);
391 
392  sub->setTransformation(closestPoint(newX));
393  int newReg = sub->isWhere(newX);
394  if (newReg+1) { // see what region we are in
395  if (newmed && newReg >= 0) { // in subgeom, if its not -1
396  *newmed = sub->medium(newReg); // then return the proper
397  }
398  //if (t<0) { // media and region
399  // egsWarning("Returning negative t from region %d (not ind) in base geom of %s\n", ireg, getName().c_str());
400  // egsWarning("howfar(%7d,[%e,%e,%e],[%e,%e,%e],%e)\n",ireg, x.x, x.y, x.z, u.x, u.y, u.z, t);
401  //}
402  return newReg + base->regions();
403  }
404  }
405 
406  if (newmed && tempReg >= 0) {
407  *newmed = base->medium(tempReg);
408  }
409  //if (t<0) {
410  // egsWarning("Returning negative t from region %d (not ind) in base geom of %s\n", ireg, getName().c_str());
411  // egsWarning("howfar(%7d,[%e,%e,%e],[%e,%e,%e],%e)\n",ireg, x.x, x.y, x.z, u.x, u.y, u.z, t);
412  //}
413  return tempReg;
414  }
415  };
416 
417  EGS_Float hownear(int ireg, const EGS_Vector &x) {
418  EGS_Float temp, dist = base->hownear(ind,x);
419  if (ireg >= base->regions()) {
420  sub->setTransformation(closestPoint(x));
421  temp = sub->hownear(ireg-base->regions(),x);
422  dist=(temp<dist)?temp:dist;
423  }
424  else if (ireg == ind) {
425  // Check all nearby geometries
426  EGS_Float xa = floor(x.x/a)*a, yb = floor(x.y/b)*b, zc = floor(x.z/c)*c;
427  EGS_Vector x0[8] = {EGS_Vector(xa, yb, zc),
428  EGS_Vector(xa+a,yb, zc),
429  EGS_Vector(xa, yb+b,zc),
430  EGS_Vector(xa+a,yb+b,zc),
431  EGS_Vector(xa, yb, zc+c),
432  EGS_Vector(xa+a,yb, zc+c),
433  EGS_Vector(xa, yb+b,zc+c),
434  EGS_Vector(xa+a,yb+b,zc+c)
435  };
436  for (int i = 0; i < 8; i++) {
437  sub->setTransformation(x0[i]);
438  temp = sub->hownear(-1,x);
439  if (temp < dist) {
440  dist = temp;
441  }
442  }
443  }
444  else {
445  dist = base->hownear(ireg,x);
446  }
447  return dist;
448  };
449 
450  int regions() {
451  return base->regions() + sub->regions();
452  }
453 
454  int getMaxStep() const {
455  return maxStep;
456  };
457 
458  const string &getType() const {
459  return type;
460  };
461 
462  void printInfo() const;
463 
464  // I have no idea what this boolean stuff is for
465  bool hasBooleanProperty(int ireg, EGS_BPType prop) const {
466  if (ireg < base->regions()) { // If ireg is less than base regions
467  return base->hasBooleanProperty(ireg, prop); // check against base regions
468  }
469  return sub->hasBooleanProperty(ireg - base->regions(), prop); // then check sub regions
470  };
471  void setBooleanProperty(EGS_BPType prop) {
472  base->setBooleanProperty(prop);
473  };
474  void addBooleanProperty(int bit) {
475  base->addBooleanProperty(bit);
476  };
477  void setBooleanProperty(EGS_BPType prop, int start, int end, int step=1) {
478  base->setBooleanProperty(prop,start,end,step);
479  };
480  void addBooleanProperty(int bit, int start, int end, int step=1) {
481  base->addBooleanProperty(bit,start,end,step);
482  };
483 
484  EGS_Float getRelativeRho(int ireg) const {
485  if (ireg < base->regions()) { // If ireg is less than base regions
486  return base->getRelativeRho(ireg); // check against base regions
487  }
488  return sub->getRelativeRho(ireg - base->regions()); // then check sub regions
489  };
490  void setRelativeRho(int start, int end, EGS_Float rho);
491  void setRelativeRho(EGS_Input *);
492 
493  void setBScaling(int start, int end, EGS_Float rho);
494  void setBScaling(EGS_Input *);
495  EGS_Float getBScaling(int ireg) const {
496  if (ireg < base->regions()) { // If ireg is less than base regions
497  return base->getBScaling(ireg); // check against base regions
498  }
499  return sub->getBScaling(ireg - base->regions()); // then check sub regions
500  };
501 
502  virtual void getLabelRegions(const string &str, vector<int> &regs);
503 
504 protected:
505  void setMedia(EGS_Input *inp,int,const int *);
506 };
507 
508 class EGS_LATTICE_EXPORT EGS_Hexagonal_Lattice : public EGS_BaseGeometry {
509 protected:
510 
513  int ind;
514  int maxStep;
515  EGS_Float a;
516  vector<EGS_Float> d;
517  string type;
518  EGS_Float gap;
519  // y and z have a shorter recurrence scale than a
520 public:
521 
522  EGS_Hexagonal_Lattice(EGS_BaseGeometry *B, EGS_BaseGeometry *S, int i, EGS_Float x, const string &Name = "");
524 
525  // Based off the closestPoint() function in the egs_lattice definition, here we check against the four
526  // Bravais lattices that make up the hexagonal lattice to find the closest point
527  EGS_Vector closestPoint(const EGS_Vector &x) {
528  EGS_Float i1,j1,k1,i2,j2,j3,k3,j4;
529  EGS_Float r3a = 2.0*gap;
530 
531  i1 = a*(round(x.x/a));
532  j1 = r3a*(round(x.y/r3a));
533  k1 = r3a*(round(x.z/r3a));
534 
535  i2 = a*(0.5+round(x.x/a-0.5));
536  j2 = r3a*(0.5+round(x.y/r3a-0.5));
537 
538  j3 = r3a*(0.25+round(x.y/r3a-0.25));
539  k3 = r3a*(0.5+round(x.z/r3a-0.5));
540 
541  j4 = r3a*(-0.25+round(x.y/r3a+0.25));
542 
543  EGS_Vector p1(i1,j1,k1);
544  EGS_Vector p2(i2,j2,k1);
545  EGS_Vector p3(i1,j3,k3);
546  EGS_Vector p4(i2,j4,k3);
547 
548  d[0] = (x-p1).length2();
549  d[1] = (x-p2).length2();
550  d[2] = (x-p3).length2();
551  d[3] = (x-p4).length2();
552 
553  if (d[0] < d[1] && d[0] < d[2] && d[0] < d[3]) {
554  return p1;
555  }
556  if (d[1] < d[2] && d[1] < d[3]) {
557  return p2;
558  }
559  if (d[2] < d[3]) {
560  return p3;
561  }
562  return p4;
563  };
564 
565  int computeIntersections(int ireg, int n, const EGS_Vector &x,
566  const EGS_Vector &u, EGS_GeometryIntersections *isections) {
567  return base->computeIntersections(ireg,n,x,u,isections);
568  };
569 
570  bool isRealRegion(int ireg) const {
571  if (ireg < base->regions()) { // If ireg is less than base regions
572  return base->isRealRegion(ireg); // check against base regions
573  }
574  return sub->isRealRegion(ireg - base->regions()); // then check sub regions
575  };
576 
577  bool isInside(const EGS_Vector &x) {
578  return base->isInside(x);
579  };
580 
581  int isWhere(const EGS_Vector &x) {
582  int temp = base->isWhere(x);
583  if (temp == ind) {// Are we in the subgeom?
584  sub->setTransformation(closestPoint(x));
585  if (sub->isInside(x)) {
586  return sub->isWhere(x) + base->regions();
587  }
588  }
589  return temp; // otherwise base geom
590  };
591 
592  int inside(const EGS_Vector &x) {
593  return isWhere(x);
594  };
595 
596  int medium(int ireg) const {
597  if (ireg >= base->regions()) { // If ireg is greater than base regions
598  return sub->medium(ireg-base->regions());
599  }
600  return base->medium(ireg); // If in base region
601  };
602 
603  EGS_Float howfarToOutside(int ireg, const EGS_Vector &x, const EGS_Vector &u) {
604  if (ireg < 0) { // Error catch, not inside
605  return 0;
606  }
607 
608  // Are we in the subgeom
609  if (ireg >= base->regions()) {
610  return base->howfarToOutside(ind,x,u);
611  }
612  return base->howfarToOutside(ireg,x,u);
613  };
614 
615  int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u,
616  EGS_Float &t,int *newmed=0, EGS_Vector *normal=0) {
617  // Catch t=0 exception, which sometimes causes issues when ind is the region neighboring -1
618  if (t < epsilon) {
619  t = 0.0;
620  return ireg;
621  }
622 
623  if (ireg >= base->regions()) { // Are we in the subgeom?
624  sub->setTransformation(closestPoint(x));
625 
626  // Do the howfar call
627  EGS_Float tempT = t;
628  int tempReg = sub->howfar(ireg-base->regions(),x,u,tempT,newmed,normal);
629 
630  // If we do leave the subgeom, we want to return the index
631  // of the region of the base geometry we would be entering,
632  // which is not necessarily ind if subgeom is at a boundary
633  EGS_Vector newX = u;
634  newX.normalize();
635  newX = x + (newX * tempT);
636 
637  if (base->isWhere(newX) != ind) {// If we leave ind region of base
638  t = tempT;
639 
640  tempReg = base->howfar(ind,x,u,t,0,normal);
641  if (newmed && tempReg >= 0) {
642  *newmed = base->medium(tempReg);
643  }
644 
645  //if (t<0) {
646  // egsWarning("Returning negative t from subgeom of %s\n",getName().c_str());
647  // egsWarning("howfar(%7d,[%e,%e,%e],[%e,%e,%e],%e)\n",ireg, x.x, x.y, x.z, u.x, u.y, u.z, t);
648  //}
649  return tempReg;
650  }
651 
652  if (!(tempReg+1)) {// If we leave sub geom back into ind
653  if (newmed) {
654  *newmed = base->medium(ind);
655  }
656  t = tempT;
657  return ind;
658  }
659 
660  // else we stay in sub geom
661  t = tempT;
662  //if (t<0) {
663  // egsWarning("Returning negative t from subgeom of %s\n",getName().c_str());
664  // egsWarning("howfar(%7d,[%e,%e,%e],[%e,%e,%e],%e)\n",ireg, x.x, x.y, x.z, u.x, u.y, u.z, t);
665  //}
666  return tempReg+base->regions();
667  }
668  else if (ireg == ind) {// If we are in the region that could contain subgeoms
669  // Determine the path travelled ------------------------------------------------ //
670  EGS_Float tempT = t;
671  base->howfar(ireg,x,u,tempT); // Get how far it is in temp
672  EGS_Float max = tempT;
673 
674  // Iterate through the lattice ------------------------------------------------- //
675  EGS_Float min = max, minX = max, minY = max, minZ = max; // Distance to closest subgeom for three different cases
676  EGS_Vector unit = u;
677  unit.normalize();
678  EGS_Vector xInt = unit*(a/8.0/fabs(unit.x)), yInt = unit*(gap/8.0/fabs(unit.y)), zInt = unit*(gap/8.0/fabs(unit.z));
679  EGS_Vector tempP; // Temporarily hold indices of subgeom we are testing against
680  EGS_Vector finalP, finalX, finalY, finalZ; // Current closest intersecting subgeom
681 
682  EGS_Vector x0;
683  EGS_Float max2 = max*max;
684 
685  x0 = x;
686  while ((x-x0).length2() < max2) {// Quit as soon as our testing position is beyond current closest intersection (which could be tempT)
687  tempT = minX;
688  tempP = closestPoint(x0);
689  sub->setTransformation(tempP);
690  if (sub->howfar(-1,x,u,tempT)+1) // Intersection!
691  if (tempT < max && tempT > 0) {
692  finalX = tempP;
693  minX = tempT;
694  break;
695  }
696  x0 += xInt;
697  }
698 
699  x0 = x;
700  while ((x-x0).length2() < max2) {// Quit as soon as our testing position is beyond current closest intersection (which could be tempT)
701  tempT = minY;
702  tempP = closestPoint(x0);
703  sub->setTransformation(tempP);
704  if (sub->howfar(-1,x,u,tempT)+1) // Intersection!
705  if (tempT < max && tempT > 0) {
706  finalY = tempP;
707  minY = tempT;
708  break;
709  }
710  x0 += yInt;
711  }
712 
713  x0 = x;
714  while ((x-x0).length2() < max2) {// Quit as soon as our testing position is beyond current closest intersection (which could be tempT)
715  tempT = minZ;
716  tempP = closestPoint(x0);
717  sub->setTransformation(tempP);
718  if (sub->howfar(-1,x,u,tempT)+1) // Intersection!
719  if (tempT < max && tempT > 0) {
720  finalZ = tempP;
721  minZ = tempT;
722  break;
723  }
724  x0 += zInt;
725  }
726 
727  finalP = finalX;
728  min = minX;
729  if (min > minY) {
730  min = minY;
731  finalP = finalY;
732  }
733  if (min > minZ) {
734  min = minZ;
735  finalP = finalZ;
736  }
737 
738  // Check last point
739  tempT = max;
740  sub->setTransformation(closestPoint(x+unit*tempT));
741  sub->howfar(-1,x,u,tempT);
742  if (tempT < min && tempT > 0) {
743  min = tempT;
744  finalP = closestPoint(x+unit*tempT);
745  }
746 
747  // We did intersect subgeom
748  if (min < max) {
749  tempT = t;
750  sub->setTransformation(finalP);
751  int tempReg = sub->howfar(-1,x,u,tempT,newmed,normal);
752  if (tempReg+1) {
753  if (newmed && tempReg >= 0) {
754  *newmed = sub->medium(tempReg);
755  }
756  t = tempT;
757 
758  //if (t<0) {
759  // egsWarning("Returning negative t from region %d in base geom of %s\n", ind, getName().c_str());
760  // egsWarning("howfar(%7d,[%e,%e,%e],[%e,%e,%e],%e)\n",ireg, x.x, x.y, x.z, u.x, u.y, u.z, t);
761  //}
762  return tempReg+base->regions();
763  }
764  }
765 
766  // We didn't intersect subgeom
767  int tempReg = base->howfar(ireg,x,u,t,newmed,normal);
768  //if (t<0) {
769  // egsWarning("Returning negative t from region %d in base geom of %s\n", ind, getName().c_str());
770  // egsWarning("howfar(%7d,[%e,%e,%e],[%e,%e,%e],%e)\n",ireg, x.x, x.y, x.z, u.x, u.y, u.z, t);
771  //}
772  if (newmed && tempReg >= 0) {
773  *newmed = base->medium(tempReg);
774  }
775 
776  return tempReg;
777  }
778  else {// Not in region containing subgeoms, then it is quite easy
779  // unless we enter directly into a subgeom when entering
780  // region ind
781  int tempReg = base->howfar(ireg,x,u,t,newmed,normal);
782 
783  // If we enter region ind
784  if (tempReg == ind) {
785  // newX is the point of entrance
786  EGS_Vector newX = u;
787  newX.normalize();
788  newX = x + (newX * t);
789 
790  sub->setTransformation(closestPoint(newX));
791  int newReg = sub->isWhere(newX);
792  if (newReg+1) { // see what region we are in
793  if (newmed && newReg >= 0) { // in subgeom, if its not -1
794  *newmed = sub->medium(newReg); // then return the proper
795  }
796  //if (t<0) { // media and region
797  // egsWarning("Returning negative t from region %d (not ind) in base geom of %s\n", ireg, getName().c_str());
798  // egsWarning("howfar(%7d,[%e,%e,%e],[%e,%e,%e],%e)\n", ireg, x.x, x.y, x.z, u.x, u.y, u.z, t);
799  //}
800  return newReg + base->regions();
801  }
802  }
803 
804  if (newmed && tempReg >= 0) {
805  *newmed = base->medium(tempReg);
806  }
807  //if (t<0) {
808  // egsWarning("Returning negative t from region %d (not ind) in base geom of %s\n", ireg, getName().c_str());
809  // egsWarning("howfar(%7d,[%e,%e,%e],[%e,%e,%e],%e)\n", ireg, x.x, x.y, x.z, u.x, u.y, u.z, t);
810  //}
811  return tempReg;
812  }
813  };
814 
815  EGS_Float hownear(int ireg, const EGS_Vector &x) {
816  EGS_Float temp, dist = base->hownear(ireg>=base->regions()?ind:ireg,x);
817  if (ireg >= base->regions()) {
818  sub->setTransformation(closestPoint(x));
819  temp = sub->hownear(ireg-base->regions(),x);
820  dist=(temp<dist)?temp:dist;
821  }
822  else if (ireg == ind) {
823  // Check all nearby geometries, making use of the 4 Bravais lattice substructure of the geometry
824  EGS_Float i,j,k;
825  EGS_Float r3a = 2.0*gap, hgap = gap/2.0, ha = a/2.0;
826  EGS_Vector x0[8];
827 
828  // Calculate the edge of an oblique rhombic prism formed around x
829  i = a*(floor(x.x/a));
830  j = r3a*(round(x.y/r3a));
831  k = r3a*(round(x.z/r3a));
832 
833  // One of the following 4 trigonal trapezohedrons must enclose x, they all share an edge defined
834  // as (i,j,k) to (i+a,j,k) and can recurse through all space
835  if (x.y<=j && x.z<=k)
836  EGS_Vector x0[8] = {EGS_Vector(i,j,k),
837  EGS_Vector(i-ha,j-gap,k),
838  EGS_Vector(i,j-gap-hgap,k-gap),
839  EGS_Vector(i+ha,j-hgap,k-gap),
840  EGS_Vector(i+a,j,k),
841  EGS_Vector(i+ha,j-gap,k),
842  EGS_Vector(i+a,j-gap-hgap,k-gap),
843  EGS_Vector(i+a+ha,j-hgap,k-gap)
844  };
845  else if (x.y>j && x.z<=k)
846  EGS_Vector x0[8] = {EGS_Vector(i,j,k),
847  EGS_Vector(i+ha,j+gap,k),
848  EGS_Vector(i+a,j+hgap,k-gap),
849  EGS_Vector(i+ha,j-hgap,k-gap),
850  EGS_Vector(i+a,j,k),
851  EGS_Vector(i+a+ha,j+gap,k),
852  EGS_Vector(i+a+a,j+hgap,k-gap),
853  EGS_Vector(i+a+ha,j-hgap,k-gap)
854  };
855  else if (x.y<j && x.z>k)
856  EGS_Vector x0[8] = {EGS_Vector(i,j,k),
857  EGS_Vector(i+ha,j-gap,k),
858  EGS_Vector(i,j-gap-hgap,k+gap),
859  EGS_Vector(i-ha,j-hgap,k+gap),
860  EGS_Vector(i+a,j,k),
861  EGS_Vector(i+a+ha,j-gap,k),
862  EGS_Vector(i+a,j-gap-hgap,k+gap),
863  EGS_Vector(i+ha,j-hgap,k+gap)
864  };
865  else //(x.y>j && x.z>k)
866  EGS_Vector x0[8] = {EGS_Vector(i,j,k),
867  EGS_Vector(i-ha,j+gap,k),
868  EGS_Vector(i-a,j+hgap,k+gap),
869  EGS_Vector(i-ha,j-hgap,k+gap),
870  EGS_Vector(i+a,j,k),
871  EGS_Vector(i+ha,j+gap,k),
872  EGS_Vector(i,j+hgap,k+gap),
873  EGS_Vector(i+ha,j-hgap,k+gap)
874  };
875  for (int i = 0; i < 8; i++) {
876  sub->setTransformation(x0[i]);
877  temp = sub->hownear(-1,x);
878  if (temp < dist) {
879  dist = temp;
880  }
881  }
882  }
883  return dist;
884  };
885 
886  int regions() {
887  return base->regions() + sub->regions();
888  }
889 
890  int getMaxStep() const {
891  return maxStep;
892  };
893 
894  const string &getType() const {
895  return type;
896  };
897 
898  void printInfo() const;
899 
900  // I have no idea what this boolean stuff is for
901  bool hasBooleanProperty(int ireg, EGS_BPType prop) const {
902  if (ireg < base->regions()) { // If ireg is less than base regions
903  return base->hasBooleanProperty(ireg, prop); // check against base regions
904  }
905  return sub->hasBooleanProperty(ireg - base->regions(), prop); // then check sub regions
906  };
907  void setBooleanProperty(EGS_BPType prop) {
908  base->setBooleanProperty(prop);
909  };
910  void addBooleanProperty(int bit) {
911  base->addBooleanProperty(bit);
912  };
913  void setBooleanProperty(EGS_BPType prop, int start, int end, int step=1) {
914  base->setBooleanProperty(prop,start,end,step);
915  };
916  void addBooleanProperty(int bit, int start, int end, int step=1) {
917  base->addBooleanProperty(bit,start,end,step);
918  };
919 
920  EGS_Float getRelativeRho(int ireg) const {
921  if (ireg < base->regions()) { // If ireg is less than base regions
922  return base->getRelativeRho(ireg); // check against base regions
923  }
924  return sub->getRelativeRho(ireg - base->regions()); // then check sub regions
925  };
926  void setRelativeRho(int start, int end, EGS_Float rho);
927  void setRelativeRho(EGS_Input *);
928 
929  void setBScaling(int start, int end, EGS_Float rho);
930  void setBScaling(EGS_Input *);
931  EGS_Float getBScaling(int ireg) const {
932  if (ireg < base->regions()) { // If ireg is less than base regions
933  return base->getBScaling(ireg); // check against base regions
934  }
935  return sub->getBScaling(ireg - base->regions()); // then check sub regions
936  };
937 
938  virtual void getLabelRegions(const string &str, vector<int> &regs);
939 
940 protected:
941  void setMedia(EGS_Input *inp,int,const int *);
942 };
943 
944 #endif
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 ?
int regions() const
Returns the number of local regions in this geometry.
virtual int medium(int ireg) const
Returns the medium index in region ireg.
virtual void printInfo() const
Print information about this geometry.
vector< EGS_Float > d
Don&#39;t redefine 4 dists in closestPoint each invocation.
Definition: egs_lattice.h:516
EGS_Float c
The center-to-center distance along x, y, and z.
Definition: egs_lattice.h:152
virtual EGS_Float howfarToOutside(int ireg, const EGS_Vector &x, const EGS_Vector &u)
EGS_Float y
y-component
Definition: egs_vector.h:61
virtual int getMaxStep() const
Returns the maximum number of steps through the geometry.
EGS_TransformedGeometry * sub
The sub geometry that could appear within base.
Definition: egs_lattice.h:149
A class representing 3D vectors.
Definition: egs_vector.h:56
virtual int isWhere(const EGS_Vector &x)=0
In which region is poisition x?
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.
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
virtual void setRelativeRho(int start, int end, EGS_Float rho)
Set the relative mass density in regions.
A Bravais, cubic, and hexagonal lattice geometryA geometry which embeds a lattice of one geometry (na...
Definition: egs_lattice.h:145
EGS_BaseGeometry * base
The geometry within which the sub geometry appears.
Definition: egs_lattice.h:148
virtual EGS_Float getBScaling(int ireg) const
Get the B field scaling factor in region ireg.
virtual void setBScaling(int start, int end, EGS_Float bf)
Set the B field scaling factor in regions.
EGS_Float z
z-component
Definition: egs_vector.h:62
A transformed geometry.
virtual void getLabelRegions(const string &str, vector< int > &regs)
Get the list of all regions labeled with str.
void setMedia(EGS_Input *inp)
Set the media in the geometry from the input pointed to by inp.
virtual int inside(const EGS_Vector &x)=0
Returns the region index, if inside, or -1 if outside (obsolete)
virtual int computeIntersections(int ireg, int n, const EGS_Vector &x, const EGS_Vector &u, EGS_GeometryIntersections *isections)
Calculates intersection distances to region boundaries.
EGS_Float a
The center-to-center distance to the nearest 12 neighbours.
Definition: egs_lattice.h:515
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.
int ind
The region in base geom where we could encounter sub geom.
Definition: egs_lattice.h:513
EGS_Float x
x-component
Definition: egs_vector.h:60
virtual bool isRealRegion(int ireg) const
Returnes true if ireg is a real region, false otherwise.
EGS_Float gap
Translating the hexagonal lattice to xyz coordinate,.
Definition: egs_lattice.h:518
int ind
The region in base geom where we could encounter sub geom.
Definition: egs_lattice.h:150
const EGS_Float epsilon
The epsilon constant for floating point comparisons.
Definition: egs_functions.h:61
string type
The geometry type.
Definition: egs_lattice.h:517
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_TransformedGeometry * sub
The sub geometry that could appear within base.
Definition: egs_lattice.h:512
string type
The geometry type.
Definition: egs_lattice.h:153
EGS_BaseGeometry class header file.
EGS_BaseGeometry * base
The geometry within which the sub geometry appears.
Definition: egs_lattice.h:511
int maxStep
The maximum number of steps.
Definition: egs_lattice.h:151
virtual void addBooleanProperty(int bit)
Add a boolean property for the entire geometry by setting the bit&#39;th bit.
virtual void setBooleanProperty(EGS_BPType prop)
Set the boolean properties of the entire geometry to prop.
int maxStep
The maximum number of steps.
Definition: egs_lattice.h:514