EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_nd_geometry.cpp
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ nd 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 # Reid Townson
28 # Ernesto Mainegra-Hing
29 # Hubert Ho
30 # Randle Taylor
31 # Manuel Stoeckl
32 # Marc Chamberland
33 # Martin Martinov
34 # Hannah Gallop
35 #
36 ###############################################################################
37 */
38 
39 
45 #include "egs_nd_geometry.h"
46 #include "egs_input.h"
47 #include "egs_functions.h"
48 #include "egs_transformations.h"
49 
50 #ifdef HAS_GZSTREAM
51  #include "gzstream.h"
52 #endif
53 
54 #include <vector>
55 #include <fstream>
56 
57 using namespace std;
58 
59 EGS_NDGeometry::EGS_NDGeometry(int ng, EGS_BaseGeometry **G,
60  const string &Name, bool O) : EGS_BaseGeometry(Name), N(ng), ortho(O) {
61  nreg = 1;
62  setup();
63  for (int j=0; j<N; j++) {
64  g[j] = G[j];
65  n[j] = nreg;
66  nreg *= g[j]->regions();
67  if (!g[j]->isConvex()) {
68  is_convex = false;
69  }
70  }
71  n[N] = nreg;
72 }
73 
74 EGS_NDGeometry::EGS_NDGeometry(vector<EGS_BaseGeometry *> &G,
75  const string &Name, bool O) : EGS_BaseGeometry(Name), N(G.size()), ortho(O) {
76  nreg = 1;
77  setup();
78  for (int j=0; j<N; j++) {
79  g[j] = G[j];
80  n[j] = nreg;
81  nreg *= g[j]->regions();
82  if (!g[j]->isConvex()) {
83  is_convex = false;
84  }
85  }
86  n[N] = nreg;
87 }
88 
89 EGS_NDGeometry::~EGS_NDGeometry() {
90  for (int j=0; j<N; j++) {
91  int ref = g[j]->deref();
92  if (!ref) {
93  delete g[j];
94  }
95  }
96  delete [] n;
97  delete [] g;
98 }
99 
100 void EGS_NDGeometry::setup() {
101  n = new int [N+1];
102  g = new EGS_BaseGeometry* [N];
103 }
104 
105 void EGS_NDGeometry::printInfo() const {
106  EGS_BaseGeometry::printInfo();
107  egsInformation(" number of dimensions = %d\n",N);
108  for (int j=0; j<N; j++)
109  egsInformation(" dimension %d = %s (type %s)\n",
110  j+1,g[j]->getName().c_str(),g[j]->getType().c_str());
111  egsInformation(
112  "=======================================================\n");
113 }
114 
115 void EGS_NDGeometry::setMedia(EGS_Input *input, int nmed, const int *mind) {
116  EGS_Input *i;
117  med = mind[0];
118  while ((i = input->takeInputItem("set medium"))) {
119  vector<int> inp;
120  int err = i->getInput("set medium",inp);
121  delete i;
122  if (!err) {
123  if (inp.size() == 2) {
124  setMedium(inp[0],inp[0],mind[inp[1]]);
125  }
126  else if (inp.size() == 3) {
127  setMedium(inp[0],inp[1],mind[inp[2]]);
128  }
129  else if (inp.size() == 4) setMedium(inp[0],inp[1],mind[inp[2]],
130  inp[3]);
131  else if (inp.size() == 2*N+1) {
132  setM(0,0,inp,mind[inp[2*N]]);
133  }
134  else egsWarning("EGS_NDGeometry::setMedia(): found %dinputs\n"
135  "in a 'set medium' input. 2, 3, 4, or %d are allowed\n",2*N+1);
136  }
137  else egsWarning("EGS_NDGeometry::setMedia(): wrong 'set medium'"
138  " input\n");
139  }
140 }
141 
142 void EGS_NDGeometry::setM(int ibase, int idim,
143  const vector<int> &ranges, int medium) {
144  int istart = ranges[2*idim], iend = ranges[2*idim+1];
145  if (istart < 0) {
146  istart = 0;
147  }
148  int ndim = n[idim+1] / n[idim];
149  if (iend > ndim) {
150  iend = ndim;
151  }
152  if (idim < N-1) for (int j=istart; j<=iend; j++) {
153  setM(ibase+j*n[idim],idim+1,ranges,medium);
154  }
155  else for (int j=istart; j<=iend; j++) {
156  setMedium(ibase+j*n[idim],ibase+j*n[idim],medium);
157  }
158 }
159 
160 #ifdef EXPLICIT_XYZ
161 
162 string EGS_XYZGeometry::type = "EGS_XYZGeometry";
163 
164 EGS_XYZGeometry::EGS_XYZGeometry(EGS_PlanesX *Xp, EGS_PlanesY *Yp,
165  EGS_PlanesZ *Zp, const string &Name) : EGS_BaseGeometry(Name),
166  xp(Xp), yp(Yp), zp(Zp) {
167  nx = xp->regions();
168  ny = yp->regions();
169  nz = zp->regions();
170  nxy = nx*ny;
171  nreg = nxy*nz;
172  xp->ref();
173  yp->ref();
174  zp->ref();
175  xpos = xp->getPositions();
176  ypos = yp->getPositions();
177  zpos = zp->getPositions();
178 }
179 
180 EGS_XYZGeometry::~EGS_XYZGeometry() {
181  if (!xp->deref()) {
182  delete xp;
183  }
184  if (!yp->deref()) {
185  delete yp;
186  }
187  if (!zp->deref()) {
188  delete zp;
189  }
190 }
191 
192 void EGS_XYZGeometry::printInfo() const {
193  EGS_BaseGeometry::printInfo();
194  xp->printInfo();
195  yp->printInfo();
196  zp->printInfo();
197  egsInformation(
198  "=======================================================\n");
199 }
200 
201 int EGS_XYZGeometry::medIndex(char medium) {
202  std::string medIdx = "123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz";
203  size_t pos = medIdx.find(medium);
204  if (pos == std::string::npos) {
205  pos = -1; // Assume vacuum if not found
206  }
207  return static_cast<int>(pos);
208 }
209 
210 void EGS_XYZGeometry::setMedia(EGS_Input *input, int nmed, const int *mind) {
211  EGS_Input *i;
212  med = mind[0];
213  while ((i = input->takeInputItem("set medium"))) {
214  vector<int> inp;
215  int err = i->getInput("set medium",inp);
216  delete i;
217  if (!err) {
218  if (inp.size() == 2) {
219  setMedium(inp[0],inp[0],mind[inp[1]]);
220  }
221  else if (inp.size() == 3) {
222  setMedium(inp[0],inp[1],mind[inp[2]]);
223  }
224  else if (inp.size() == 4) setMedium(inp[0],inp[1],mind[inp[2]],
225  inp[3]);
226  else if (inp.size() == 7) {
227  int is = inp[0], ie = inp[1];
228  if (is < 0) {
229  is = 0;
230  }
231  if (ie > nx-1) {
232  ie = nx-1;
233  }
234  int js = inp[2], je = inp[3];
235  if (js < 0) {
236  js = 0;
237  }
238  if (je > nxy/nx-1) {
239  je = nxy/nx-1;
240  }
241  int ks = inp[4], ke = inp[5];
242  if (ks < 0) {
243  ks = 0;
244  }
245  if (ke > nreg/nxy-1) {
246  ke = nreg/nxy-1;
247  }
248  //egsInformation("setting medium to %d between %d %d %d %d "
249  // "%d %d\n",mind[inp[6]],is,ie,js,je,ks,ke);
250  for (int i=is; i<=ie; i++) {
251  for (int j=js; j<=je; j++) {
252  for (int k=ks; k<=ke; k++) {
253  int ireg = i + j*nx + k*nxy;
254  //egsWarning("set: %d %d %d %d %d\n",
255  // i,j,k,ireg,mind[inp[6]]);
256  setMedium(ireg,ireg,mind[inp[6]]);
257  }
258  }
259  }
260  }
261  else if (inp.size() == 10) {
262  int is = inp[0], ie = inp[1], idelta = inp[2];
263  if (is < 0) {
264  is = 0;
265  }
266  if (ie > nx-1) {
267  ie = nx-1;
268  }
269  if (idelta <= 0) {
270  idelta = 1;
271  }
272  int js = inp[3], je = inp[4], jdelta = inp[5];
273  if (js < 0) {
274  js = 0;
275  }
276  if (je > nxy/nx-1) {
277  je = nxy/nx-1;
278  }
279  if (jdelta <= 0) {
280  jdelta = 1;
281  }
282  int ks = inp[6], ke = inp[7], kdelta = inp[8];
283  if (ks < 0) {
284  ks = 0;
285  }
286  if (ke > nreg/nxy-1) {
287  ke = nreg/nxy-1;
288  }
289  if (kdelta <= 0) {
290  kdelta = 1;
291  }
292  //egsInformation("setting medium to %d between %d %d %d %d "
293  // "%d %d\n",mind[inp[6]],is,ie,js,je,ks,ke);
294  for (int i=is; i<=ie; i+=idelta) {
295  for (int j=js; j<=je; j+=jdelta) {
296  for (int k=ks; k<=ke; k+=kdelta) {
297  int ireg = i + j*nx + k*nxy;
298  //egsWarning("set: %d %d %d %d %d\n",
299  // i,j,k,ireg,mind[inp[6]]);
300  setMedium(ireg,ireg,mind[inp[9]]);
301  }
302  }
303  }
304  }
305  else egsWarning("EGS_XYZGeometry::setMedia(): found %dinputs\n"
306  "in a 'set medium' input. 2, 3, 4, 7, or 10 are allowed\n");
307  }
308  else egsWarning("EGS_XYZGeometry::setMedia(): wrong 'set medium'"
309  " input\n");
310  }
311 }
312 
313 EGS_XYZGeometry *EGS_XYZGeometry::constructGeometry(const char *dens_file,
314  const char *ramp_file, int dens_or_egsphant_or_interfile) {
315  const static char *func = "EGS_XYZGeometry::constructGeometry";
316  // A ramp file is only needed for density matrix or interfile
317  if (!dens_file || (dens_or_egsphant_or_interfile != 1 && !ramp_file)) {
318  return 0;
319  }
320 
321  ifstream ramp(ramp_file);
322  if (dens_or_egsphant_or_interfile != 1 && !ramp) {
323  egsWarning("%s: failed to open CT ramp file %s\n",func,ramp_file);
324  return 0;
325  }
326  vector<string> med_names;
327  vector<EGS_Float> rho_min,rho_max,rho_def;
328 
329  if (dens_or_egsphant_or_interfile == 2) {
330  // interfile format of data
331  int meds;
332  ramp >> meds;
333  for (int i=0; i < meds; i++) {
334  string medname;
335  EGS_Float rmin,rmax,rdef;
336  ramp >> rmin >> rmax >> medname;
337  rdef=(rmin+rmax) / 2.0; // is this correct definition of rdef ??
338  if (ramp.eof() || ramp.fail() || !ramp.good()) {
339  break;
340  }
341  egsInformation("Using medium %s for rho=%g...%g\n",medname.c_str(),
342  rmin,rmax);
343  med_names.push_back(medname);
344  rho_min.push_back(rmin);
345  rho_max.push_back(rmax+1);
346  rho_def.push_back(rdef);
347  }
348  if (med_names.size() < 1) {
349  egsWarning("%s: no media defined in the CT ramp "
350  "file %s!\n",func,ramp_file);
351  return 0;
352  }
353  }
354  else if (dens_or_egsphant_or_interfile == 0) {
355  // Density matrix format uses a ramp
356  for (EGS_I64 loopCount=0; loopCount<=loopMax; ++loopCount) {
357  if (loopCount == loopMax) {
358  egsFatal("EGS_XYZGeometry::constructGeometry: Too many iterations were required! Input may be invalid, or consider increasing loopMax.");
359  return 0;
360  }
361  // other format of data
362  string medname;
363  EGS_Float rmin,rmax,rdef;
364  ramp >> medname >> rmin >> rmax >> rdef;
365  if (ramp.eof() || ramp.fail() || !ramp.good()) {
366  break;
367  }
368  egsInformation("Using medium %s for rho=%g...%g\n",medname.c_str(),
369  rmin,rmax);
370  med_names.push_back(medname);
371  rho_min.push_back(rmin);
372  rho_max.push_back(rmax);
373  rho_def.push_back(rdef);
374  }
375  if (med_names.size() < 1) {
376  egsWarning("%s: no media defined in the CT ramp "
377  "file %s!\n",func,ramp_file);
378  return 0;
379  }
380  }
381 
382  int Nx, Ny, Nz;
383  EGS_Float *xx, *yy, *zz;
384  float *rho;
385  if (dens_or_egsphant_or_interfile == 0) {
386  int my_endian = egsGetEndian();
387  if (my_endian < 0) egsFatal("%s: machine has an unknown"
388  " endianess!\n",func);
389  ifstream dens(dens_file,ios::binary);
390  if (!dens) {
391  egsWarning("%s: failed to open density matrix file "
392  "%s\n",func,dens_file);
393  return 0;
394  }
395  char their_endian;
396  dens.read(&their_endian,1);
397  if (their_endian != 0 && their_endian != 1)
398  egsFatal("%s: density data created on a machine with"
399  " unknown endianess!\n",func);
400  bool swap = my_endian != their_endian;
401  dens.read((char *) &Nx,sizeof(int));
402  if (swap) {
403  egsSwapBytes(&Nx);
404  }
405  dens.read((char *) &Ny,sizeof(int));
406  if (swap) {
407  egsSwapBytes(&Ny);
408  }
409  dens.read((char *) &Nz,sizeof(int));
410  if (swap) {
411  egsSwapBytes(&Nz);
412  }
413  if (Nx < 1 || Ny < 1 || Nz < 1) {
414  egsWarning("%s: invalid density matrix file: "
415  " Nx=%d Ny=%d Nz=%d\n",func,Nx,Ny,Nz);
416  return 0;
417  }
418  float *x = new float [Nx+1];
419  float *y = new float [Ny+1];
420  float *z = new float [Nz+1];
421  dens.read((char *)x, (Nx+1)*sizeof(float));
422  if (dens.eof() || dens.fail() || !dens.good()) {
423  egsWarning("%s: error while reading %d x-planes from"
424  " file %s\n",func,Nx+1,dens_file);
425  delete [] x;
426  delete [] y;
427  delete [] z;
428  return 0;
429  }
430  dens.read((char *)y, (Ny+1)*sizeof(float));
431  if (dens.eof() || dens.fail() || !dens.good()) {
432  egsWarning("%s: error while reading %d y-planes from"
433  " file %s\n",func,Ny+1,dens_file);
434  delete [] x;
435  delete [] y;
436  delete [] z;
437  return 0;
438  }
439  dens.read((char *)z, (Nz+1)*sizeof(float));
440  if (dens.eof() || dens.fail() || !dens.good()) {
441  egsWarning("%s: error while reading %d z-planes from"
442  " file %s\n",func,Nz+1,dens_file);
443  delete [] x;
444  delete [] y;
445  delete [] z;
446  return 0;
447  }
448  if (swap) {
449  int j;
450  for (j=0; j<Nx; j++) {
451  egsSwapBytes(&x[j]);
452  }
453  for (j=0; j<Ny; j++) {
454  egsSwapBytes(&y[j]);
455  }
456  for (j=0; j<Nz; j++) {
457  egsSwapBytes(&z[j]);
458  }
459  }
460  {
461  int j;
462  xx = new EGS_Float [Nx+1];
463  for (j=0; j<=Nx; j++) {
464  xx[j] = x[j];
465  }
466  yy = new EGS_Float [Ny+1];
467  for (j=0; j<=Ny; j++) {
468  yy[j] = y[j];
469  }
470  zz = new EGS_Float [Nz+1];
471  for (j=0; j<=Nz; j++) {
472  zz[j] = z[j];
473  }
474  delete [] x;
475  delete [] y;
476  delete [] z;
477  }
478  rho = new float [Nx*Ny*Nz];
479  dens.read((char *)rho, Nx*Ny*Nz*sizeof(float));
480  if (dens.fail()) {
481  egsWarning("%s: failed reading mass densities from density matrux file\n",func);
482  delete [] rho;
483  delete [] xx;
484  delete [] yy;
485  delete [] zz;
486  return 0;
487  }
488  if (swap) {
489  for (int j=0; j<Nx*Ny*Nz; j++) {
490  egsSwapBytes(&rho[j]);
491  }
492  }
493  }
494  else if (dens_or_egsphant_or_interfile == 2) {
495  ifstream h_file(dens_file);
496  if (!h_file) {
497  egsWarning("%s: failed to open the interfile header %s\n",func,dens_file);
498  return 0;
499  }
500  string data_file("");
501  int data_type = -1;
502  float scale_x=0.0, scale_y=0.0, scale_z=0.0;
503  for (EGS_I64 loopCount=0; loopCount<=loopMax; ++loopCount) {
504  if (loopCount == loopMax) {
505  egsFatal("EGS_XYZGeometry::constructGeometry: Too many iterations were required! Input may be invalid, or consider increasing loopMax.");
506  return 0;
507  }
508  string line, key, value;
509  size_t pos;
510  getline(h_file, line);
511  if (h_file.eof() || h_file.fail() || !h_file.good()) {
512  break;
513  }
514 
515  pos = line.find(":=");
516  if (pos == string::npos) {
517  continue;
518  }
519  key = line.substr(0, int(pos));
520  value = line.substr(int(pos)+2);
521  while ((key[0] == '!') || (key[0] == ' ')) {
522  key.erase(0,1);
523  }
524  while (key[key.length()-1] == ' ') {
525  key.erase(key.length()-1, 1);
526  }
527  while (value[0] == ' ') {
528  value.erase(0,1);
529  }
530  while (value[value.length()-1] == ' ') {
531  value.erase(value.length()-1, 1);
532  }
533 
534  if (key == "matrix size [1]") {
535  sscanf(value.c_str(), "%u", &Nx);
536  }
537  else if (key == "matrix size [2]") {
538  sscanf(value.c_str(), "%u", &Ny);
539  }
540  else if ((key == "number of slices") || (key == "number of images")) {
541  sscanf(value.c_str(), "%u", &Nz);
542  }
543  else if (key == "scaling factor (mm/pixel) [1]") {
544  sscanf(value.c_str(), "%f", &scale_x);
545  }
546  else if (key == "scaling factor (mm/pixel) [2]") {
547  sscanf(value.c_str(), "%f", &scale_y);
548  }
549  else if (key == "slice thickness (pixels)") {
550  sscanf(value.c_str(), "%f", &scale_z); // ????????????????
551  }
552  else if (key == "name of data file") {
553  data_file = value;
554  }
555  else if (key == "number format") {
556  if ((value == "float") || (value == "FLOAT")) {
557  data_type = 0;
558  }
559  else if ((value == "unsigned integer") || (value == "UNSIGNED INTEGER")) {
560  data_type = 1;
561  }
562  else {
563  egsWarning("%s: unrecognised 'number format' type: %s \n",func,value.c_str());
564  }
565  }
566  }
567  if (Nx < 1 || Ny < 1 || Nz < 1 || scale_x <= 0.0 || scale_y <= 0.0 || scale_z <= 0.0 || data_file == "" || data_type == -1) {
568  egsWarning("%s: invalid interfile header information: "
569  "Nx=%d Ny=%d Nz=%d scale_x=%f scale_y=%f scale_z=%f "
570  "data_file='%s' number_format=%d\n",func,Nx,Ny,Nz,scale_x,scale_y,scale_z,data_file.c_str(),data_type);
571  return 0;
572  }
573  ifstream i_file(data_file.c_str(),ios::binary);
574  if (!i_file) {
575  egsWarning("%s: failed to open interfile data "
576  "%s\n",func,data_file.c_str());
577  return 0;
578  }
579  scale_x /= 10.0;
580  scale_y /= 10.0;
581  scale_z /= 10.0; // convert to cm
582  {
583  int j;
584  xx = new EGS_Float [Nx+1];
585  for (j=0; j<=Nx; j++) {
586  xx[j] = -(float)Nx*scale_x/2.0 + j*scale_x;
587  }
588  yy = new EGS_Float [Ny+1];
589  for (j=0; j<=Ny; j++) {
590  yy[j] = -(float)Ny*scale_y/2.0 + j*scale_y;
591  }
592  zz = new EGS_Float [Nz+1];
593  for (j=0; j<=Nz; j++) {
594  zz[j] = -(float)Nz*scale_z/2.0 + j*scale_z;
595  }
596  }
597  rho = new float [Nx*Ny*Nz];
598  if (data_type == 0) {
599  i_file.read((char *)rho, Nx*Ny*Nz*sizeof(float));
600  }
601  else {
602  unsigned short int *rho_tmp = new unsigned short int [Nx*Ny*Nz];
603  i_file.read((char *)rho_tmp, (Nx*Ny*Nz)*sizeof(unsigned short int));
604  for (int cc = 0; cc<Nx*Ny*Nz; cc++) {
605  rho[cc] = (float)(rho_tmp[cc]);
606  }
607  delete [] rho_tmp;
608  }
609  if (i_file.fail()) {
610  egsWarning("%s: failed reading interfile data\n",func);
611  delete [] rho;
612  delete [] xx;
613  delete [] yy;
614  delete [] zz;
615  return 0;
616  }
617  }
618  else {
619 
620  ifstream tempf(dens_file, ios::binary);
621  istream *data;
622  ifstream textf;
623 #ifdef HAS_GZSTREAM
624  igzstream binf;
625 #endif
626 
627  if (!tempf) {
628  egsWarning("%s: failed to open .egsphant file %s\n",func,dens_file);
629  return 0;
630  }
631 
632  bool is_gzip = (tempf.get() == 0x1f && tempf.get() == 0x8b);
633  tempf.close();
634 
635  if (is_gzip) {
636 #ifdef HAS_GZSTREAM
637  binf.open(dens_file);
638  data = &binf;
639 #else
640  egsWarning("Tried to read gzipped egsphant but egs_ndgeometry was not compiled with gzip support\n");
641  return 0;
642 #endif
643  }
644  else {
645  textf.open(dens_file);
646  data = &textf;
647  }
648 
649 
650  int nmed;
651  (*data) >> nmed;
652  if ((*data).fail()) {
653  egsWarning("%s: failed reading number of media\n",func);
654  return 0;
655  }
656  char buf [1024];
657  int imed;
658  (*data).getline(buf,1023);
659  for (imed=0; imed<nmed; ++imed) {
660  (*data).getline(buf,1023);
661 
662  if ((*data).fail()) {
663  egsWarning("%s: failed reading medium name for %d'th medium\n",func,imed+1);
664  }
665 
666  med_names.push_back(string(buf));
667  addMedium(med_names.back());
668  egsInformation("Using medium %s as mednum %d\n",
669  med_names.back().c_str(),med_names.size());
670  }
671  // ignore estepe
672  EGS_Float dum;
673  for (imed=0; imed<nmed; ++imed) {
674  (*data) >> dum;
675  }
676  (*data) >> Nx >> Ny >> Nz;
677  if ((*data).fail()) {
678  egsWarning("%s: failed reading number of voxels\n",func);
679  return 0;
680  }
681  xx = new EGS_Float [Nx+1];
682  yy = new EGS_Float [Ny+1];
683  zz = new EGS_Float [Nz+1];
684  int j;
685  for (j=0; j<=Nx; ++j) {
686  (*data) >> xx[j];
687  }
688  if ((*data).fail()) {
689  egsWarning("%s: failed reading x-planes\n",func);
690  delete [] xx;
691  delete [] yy;
692  delete [] zz;
693  return 0;
694  }
695  for (j=0; j<=Ny; ++j) {
696  (*data) >> yy[j];
697  }
698  if ((*data).fail()) {
699  egsWarning("%s: failed reading y-planes\n",func);
700  delete [] xx;
701  delete [] yy;
702  delete [] zz;
703  return 0;
704  }
705  for (j=0; j<=Nz; ++j) {
706  (*data) >> zz[j];
707  }
708  if ((*data).fail()) {
709  egsWarning("%s: failed reading z-planes\n",func);
710  delete [] xx;
711  delete [] yy;
712  delete [] zz;
713  return 0;
714  }
715  int nr = Nx*Ny*Nz;
716 
717  EGS_PlanesX *xp = new EGS_PlanesX(Nx+1,xx,"",EGS_XProjector("x-planes"));
718  EGS_PlanesY *yp = new EGS_PlanesY(Ny+1,yy,"",EGS_YProjector("y-planes"));
719  EGS_PlanesZ *zp = new EGS_PlanesZ(Nz+1,zz,"",EGS_ZProjector("z-planes"));
720  EGS_XYZGeometry *result = new EGS_XYZGeometry(xp,yp,zp);
721 
722  for (int i=0; i<Nz; ++i) {
723  for (j=0; j<Ny; ++j) {
724  string medLine;
725  //(*data).getline(buf,1023);
726  (*data) >> medLine;
727 
728  if (medLine.size() < Nx) {
729  egsWarning("%s: Number of media per line in correct.\n",func);
730  delete [] xx;
731  delete [] yy;
732  delete [] zz;
733  delete result;
734  return 0;
735  }
736  for (int k=0; k<Nx; ++k) {
737  int index = i * (Ny * Nx) + j * Nx + k;
738 
739  string medName;
740  int medIndex = result->medIndex(medLine[k]);
741  if (medIndex > -1) {
742  medName = med_names[medIndex];
743  }
744  else {
745  medName = "vacuum";
746  }
747 
748  result->setMedium(index,index,medName);
749  //egsInformation("%d",result->medIndex(medLine[k]));
750  }
751  //egsInformation("\n");
752  }
753  //egsInformation("\n");
754  (*data).getline(buf,1023); // Skip a blank line
755  }
756 
757  //int med;
758  //for(j=0; j<nr; ++j) data >> med;
759  /*(*data).getline(buf,1023);
760  for (int iz=0; iz<Nz; ++iz) {
761  for (int iy=0; iy<Ny; ++iy) {
762  (*data).getline(buf,1023);
763  }
764  (*data).getline(buf,1023);
765  }*/
766  if ((*data).fail()) {
767  egsWarning("%s: failed reading media indeces matrix\n",func);
768  delete [] xx;
769  delete [] yy;
770  delete [] zz;
771  delete result;
772  return 0;
773  }
774 
775  return result;
776  }
777 
778  // Only get here for density matrix or interfile
779  EGS_PlanesX *xp = new EGS_PlanesX(Nx+1,xx,"",EGS_XProjector("x-planes"));
780  EGS_PlanesY *yp = new EGS_PlanesY(Ny+1,yy,"",EGS_YProjector("y-planes"));
781  EGS_PlanesZ *zp = new EGS_PlanesZ(Nz+1,zz,"",EGS_ZProjector("z-planes"));
782  EGS_XYZGeometry *result = new EGS_XYZGeometry(xp,yp,zp);
783 
784  EGS_Float rhomin=1e30,rhomax=0;
785  int j;
786  for (j=0; j<Nx*Ny*Nz; j++) {
787  if (rho[j] > rhomax) {
788  rhomax = rho[j];
789  }
790  if (rho[j] < rhomin) {
791  rhomin = rho[j];
792  }
793  }
794  egsInformation("Min. density found: %g\n",rhomin);
795  egsInformation("Max. density found: %g\n",rhomax);
796  int nmed = med_names.size();
797  int *imed = new int [nmed];
798  for (j=0; j<nmed; j++) {
799  imed[j] = -2;
800  }
801  for (j=0; j<Nx*Ny*Nz; j++) {
802  int i;
803  for (i=0; i<nmed; i++) {
804  if (rho[j] >= rho_min[i] && rho[j] < rho_max[i]) {
805  break;
806  }
807  }
808  if (i >= nmed) {
809  egsWarning("%s: no material defined for density %g, leaving voxel %d"
810  " as vacuum\n",func,rho[j],j);
811  result->setMedium(j,j,-1);
812  }
813  else {
814  if (imed[i] < -1) {
815  imed[i] = addMedium(med_names[i]);
816  egsInformation("Using medium %s as mednum %d\n",
817  med_names[i].c_str(),imed[i]);
818  }
819  result->setMedium(j,j,imed[i]);
820  if (rho_def[i] > 0) {
821  EGS_Float rrho = rho[j]/rho_def[i];
822  if (fabs(rrho-1) > epsilon) {
823  result->setRelativeRho(j,j,rrho);
824  }
825  }
826  }
827  }
828  delete [] rho;
829  delete [] imed;
830  return result;
831 }
832 
833 void EGS_XYZGeometry::finishInitialization() {
834 
835  // Return if this isn't egsphant input
836  if (dens_or_egsphant_or_interfile != 1 || useEgsphantDensities == false) {
837  return;
838  }
839  const static char *func = "EGS_XYZGeometry::finishInitialization";
840 
841  ifstream tempf(dens_file, ios::binary);
842  istream *data;
843  ifstream textf;
844 #ifdef HAS_GZSTREAM
845  igzstream binf;
846 #endif
847 
848  if (!tempf) {
849  egsWarning("%s: failed to open .egsphant file %s\n",func,dens_file.c_str());
850  return;
851  }
852 
853  bool is_gzip = (tempf.get() == 0x1f && tempf.get() == 0x8b);
854  tempf.close();
855 
856  if (is_gzip) {
857 #ifdef HAS_GZSTREAM
858  binf.open(dens_file);
859  data = &binf;
860 #else
861  egsWarning("Tried to read gzipped egsphant but egs_ndgeometry was not compiled with gzip support\n");
862  return;
863 #endif
864  }
865  else {
866  textf.open(dens_file);
867  data = &textf;
868  }
869 
870  int nmed;
871  (*data) >> nmed;
872  if ((*data).fail()) {
873  egsWarning("%s: failed reading number of media\n",func);
874  return;
875  }
876 
877  char buf [1024];
878  int imed;
879  (*data).getline(buf,1023);
880  for (imed=0; imed<nmed; ++imed) {
881  (*data).getline(buf,1023);
882 
883  if ((*data).fail()) {
884  egsWarning("%s: failed reading medium name for %d'th medium\n",func,imed+1);
885  }
886  }
887 
888  // ignore estepe
889  EGS_Float dum;
890  for (imed=0; imed<nmed; ++imed) {
891  (*data) >> dum;
892  }
893  int Nx, Ny, Nz;
894  (*data) >> Nx >> Ny >> Nz;
895  if ((*data).fail()) {
896  egsWarning("%s: failed reading number of voxels\n",func);
897  return;
898  }
899 
900  int j;
901  for (j=0; j<=Nx; ++j) {
902  (*data) >> dum;
903  }
904  if ((*data).fail()) {
905  egsWarning("%s: failed reading x-planes\n",func);
906  return;
907  }
908  for (j=0; j<=Ny; ++j) {
909  (*data) >> dum;
910  }
911  if ((*data).fail()) {
912  egsWarning("%s: failed reading y-planes\n",func);
913  return;
914  }
915  for (j=0; j<=Nz; ++j) {
916  (*data) >> dum;
917  }
918  if ((*data).fail()) {
919  egsWarning("%s: failed reading z-planes\n",func);
920  return;
921  }
922  int nr = Nx*Ny*Nz;
923 
924  for (int i=0; i<Nz; ++i) {
925  for (j=0; j<Ny; ++j) {
926  //string medLine;
927  (*data).getline(buf,1023);
928  //(*data) >> medLine;
929  }
930  //egsInformation("\n");
931  (*data).getline(buf,1023); // Skip a blank line
932  }
933 
934  for (j=0; j<nr; ++j) {
935  float rho;
936  (*data) >> rho;
937 
938  EGS_Float rrho = rho / getMediumRho(medium(j));
939  if (fabs(rrho-1) > epsilon) {
940  setRelativeRho(j,j,rrho);
941  }
942  }
943  if ((*data).fail()) {
944  egsWarning("%s: failed reading mass density matrix\n",func);
945  return;
946  }
947 
948  return;
949 }
950 
951 string EGS_DeformedXYZ::def_type = "EGS_DeformedXYZ";
952 
953 char EGS_DeformedXYZ::tetrahedra[] = {0,2,3,6, 0,6,3,7, 0,4,6,7,
954  7,0,1,3, 7,1,0,4, 7,5,1,4
955  };
956 
957 char EGS_DeformedXYZ::plane_order[] = {0,1,2, 0,2,3, 0,3,1, 1,3,2};
958 
959 char EGS_DeformedXYZ::enter_tetra1[] = {0, 3, 4, 0, 0, 2};
960 char EGS_DeformedXYZ::enter_plane1[] = {1, 1, 1, 3, 0, 3};
961 char EGS_DeformedXYZ::enter_tetra2[] = {2, 5, 5, 1, 3, 5};
962 char EGS_DeformedXYZ::enter_plane2[] = {0, 0, 2, 3, 2, 2};
963 
964 
965 char EGS_DeformedXYZ::tetra_data[] = {
966  2,-1, 2, // plane 0 in tetrahedron 0
967  -1, 0, 1, // plane 1 in tetrahedron 0
968  0,-1, 3, // plane 2 in tetrahedron 0
969  1, 1, 4, // plane 3 in tetrahedron 0
970 //------------------------------------------------------------------
971  -1, 0, 0, // plane 0 in tetrahedron 1
972  -1, 0, 3, // plane 1 in tetrahedron 1
973  -1, 0, 2, // plane 2 in tetrahedron 1
974  1, 1, 5, // plane 3 in tetrahedron 1
975 //------------------------------------------------------------------
976  0,-1, 5, // plane 0 in tetrahedron 2
977  -1, 0, 1, // plane 1 in tetrahedron 2
978  -1, 0, 4, // plane 2 in tetrahedron 2
979  2, 1, 0, // plane 3 in tetrahedron 2
980 //------------------------------------------------------------------
981  -1, 0, 4, // plane 0 in tetrahedron 3
982  0, 1, 0, // plane 1 in tetrahedron 3
983  -1, 0, 1, // plane 2 in tetrahedron 3
984  2,-1, 5, // plane 3 in tetrahedron 3
985 //------------------------------------------------------------------
986  -1, 0, 3, // plane 0 in tetrahedron 4
987  -1, 0, 2, // plane 1 in tetrahedron 4
988  -1, 0, 5, // plane 2 in tetrahedron 4
989  1,-1, 0, // plane 3 in tetrahedron 4
990 //------------------------------------------------------------------
991  0, 1, 2, // plane 0 in tetrahedron 5
992  -1, 0, 4, // plane 1 in tetrahedron 5
993  2, 1, 3, // plane 2 in tetrahedron 5
994  1,-1, 1 // plane 3 in tetrahedron 5
995 };
996 
997 
998 EGS_DeformedXYZ::EGS_DeformedXYZ(EGS_PlanesX *Xp, EGS_PlanesY *Yp,
999  EGS_PlanesZ *Zp, const char *defFile, const string &Name) :
1000  EGS_XYZGeometry(Xp,Yp,Zp,Name), vectors(0) {
1001  setDeformations(defFile);
1002  np[0] = nx;
1003  np[1] = ny;
1004  np[2] = nz;
1005  nxyp1 = nx + ny + 1;
1006 }
1007 
1008 int EGS_DeformedXYZ::setDeformations(const char *defFile) {
1009  ifstream data(defFile,ios::binary);
1010  if (!data) {
1011  egsWarning("Failed to open deformations file %s\n",defFile);
1012  return 1;
1013  }
1014  int nvec;
1015  data.read((char *)&nvec,sizeof(int));
1016  int nneed = (nx+1)*(ny+1)*(nz+1);
1017  if (nvec != nneed) {
1018  egsWarning("Inconsistent deformations file: %d instead of %d vectors\n",
1019  nvec,nneed);
1020  return 2;
1021  }
1022  if (!vectors) {
1023  vectors = new EGS_Vector [nvec];
1024  }
1025  float tmp[3];
1026  int i,j,k;
1027  for (j=0; j<nvec; ++j) {
1028  data.read((char *)tmp,3*sizeof(float));
1029  if (data.fail()) {
1030  egsWarning("Error while readinf vector %d from file\n",j+1);
1031  return 3;
1032  }
1033  vectors[j] = EGS_Vector(tmp[0],tmp[1],tmp[2]);
1034  }
1035 
1036  // add plane positions to vectors
1037  for (k=0; k<=nz; ++k) for (j=0; j<=ny; ++j) for (i=0; i<=nx; ++i) {
1038  int index = i + j*(nx+1) + k*(nx+1)*(ny+1);
1039  vectors[index] += EGS_Vector(xpos[i],ypos[j],zpos[k]);
1040  }
1041 
1042  for (j=0; j<24; ++j) {
1043  int i = tetrahedra[j];
1044  switch (i) {
1045  case 0:
1046  tnodes[j] = 0;
1047  break;
1048  case 1:
1049  tnodes[j] = 1;
1050  break;
1051  case 2:
1052  tnodes[j] = nx+1;
1053  break;
1054  case 3:
1055  tnodes[j] = nx+2;
1056  break;
1057  case 4:
1058  tnodes[j] = (nx+1)*(ny+1);
1059  break;
1060  case 5:
1061  tnodes[j] = (nx+1)*(ny+1)+1;
1062  break;
1063  case 6:
1064  tnodes[j] = (nx+1)*(ny+1)+nx+1;
1065  break;
1066  case 7:
1067  tnodes[j] = (nx+1)*(ny+1)+nx+2;
1068  break;
1069  }
1070  }
1071 
1072  nreg = 6*nx*ny*nz;
1073 
1074  return 0;
1075 }
1076 
1077 EGS_DeformedXYZ::~EGS_DeformedXYZ() {
1078  if (vectors) {
1079  delete [] vectors;
1080  }
1081 }
1082 
1083 string EGS_XYZRepeater::type = "EGS_XYZRepeater";
1084 
1085 EGS_XYZRepeater::EGS_XYZRepeater(EGS_Float xmin, EGS_Float xmax,
1086  EGS_Float ymin, EGS_Float ymax,
1087  EGS_Float zmin, EGS_Float zmax,
1088  int Nx, int Ny, int Nz,
1089  EGS_BaseGeometry *G,
1090  const string &Name) :
1091  EGS_BaseGeometry(Name), g(G), nx(Nx), ny(Ny), nz(Nz) {
1092  dx = (xmax - xmin)/nx;
1093  dxi = 1/dx;
1094  ax = -xmin*dxi;
1095  dy = (ymax - ymin)/ny;
1096  dyi = 1/dy;
1097  ay = -ymin*dyi;
1098  dz = (zmax - zmin)/nz;
1099  dzi = 1/dz;
1100  az = -zmin*dzi;
1101  EGS_PlanesX *xp;
1102  EGS_PlanesY *yp;
1103  EGS_PlanesZ *zp;
1104  xp = new EGS_PlanesX(xmin,dx,nx,"",EGS_XProjector("x-planes"));
1105  yp = new EGS_PlanesY(ymin,dy,ny,"",EGS_YProjector("y-planes"));
1106  zp = new EGS_PlanesZ(zmin,dz,nz,"",EGS_ZProjector("z-planes"));
1107  xyz = new EGS_XYZGeometry(xp,yp,zp);
1108  xyz->ref();
1109  g->ref();
1110  nxy = nx*ny;
1111  nxyz = nx*ny*nz;
1112  ng = g->regions();
1113  nreg = nxyz*ng + 1;
1114  translation = new EGS_Vector [nxyz];
1115  for (int ix=0; ix<nx; ix++)
1116  for (int iy=0; iy<ny; iy++)
1117  for (int iz=0; iz<nz; iz++)
1118  translation[ix+iy*nx+iz*nxy] = EGS_Vector(xmin + dx*(0.5+ix),
1119  ymin + dy*(0.5+iy),
1120  zmin + dz*(0.5+iz));
1121 
1122 }
1123 
1124 EGS_XYZRepeater::~EGS_XYZRepeater() {
1125  if (!g->deref()) {
1126  delete g;
1127  }
1128  if (!xyz->deref()) {
1129  delete xyz;
1130  }
1131  delete [] translation;
1132 }
1133 
1134 void EGS_XYZRepeater::printInfo() const {
1135  EGS_BaseGeometry::printInfo();
1136  egsInformation("%d repetitions between xmin=%g and xmax=%g\n",
1137  nx,xyz->getXPositions()[0],xyz->getXPositions()[nx]);
1138  egsInformation("%d repetitions between ymin=%g and ymax=%g\n",
1139  ny,xyz->getYPositions()[0],xyz->getYPositions()[ny]);
1140  egsInformation("%d repetitions between zmin=%g and zmax=%g\n",
1141  nz,xyz->getZPositions()[0],xyz->getZPositions()[nz]);
1142  const char *med_name = getMediumName(med);
1143  if (med_name) {
1144  egsInformation("box is filled with %s\n",med_name);
1145  }
1146  else {
1147  egsInformation("box is filled with vacuum\n");
1148  }
1149  egsInformation("Repeated geometry:\n");
1150  g->printInfo();
1151 }
1152 
1153 void EGS_XYZGeometry::voxelizeGeometry(EGS_Input *input) {
1154  string gname;
1155  int err = input->getInput("voxelize geometry",gname);
1156  if (err) {
1157  return;
1158  }
1159  egsInformation("\nEGS_XYZGeometry(%s)\n",name.c_str());
1160  egsInformation(" setting media from geometry %s\n",gname.c_str());
1161  EGS_BaseGeometry *geometry = EGS_BaseGeometry::getGeometry(gname);
1162  if (!geometry) {
1163  egsInformation(" this geometry does not exist -> media will be not set.\n");
1164  return;
1165  }
1166  EGS_AffineTransform *T = EGS_AffineTransform::getTransformation(input);
1167  if (region_media) {
1168  delete [] region_media;
1169  }
1170  if (rhor) {
1171  delete [] rhor;
1172  rhor = 0;
1173  }
1174  if (bfactor) {
1175  delete [] bfactor;
1176  bfactor = 0;
1177  }
1178  region_media = new short [nreg];
1179 
1180  bool hrs = geometry->hasRhoScaling();
1181  if (hrs) {
1182  has_rho_scaling = true;
1183  rhor = new EGS_Float [nreg];
1184  }
1185  else {
1186  has_rho_scaling = false;
1187  }
1188 
1189  bool hbs = geometry->hasBScaling();
1190  if (hbs) {
1191  has_B_scaling = true;
1192  bfactor = new EGS_Float [nreg];
1193  }
1194  else {
1195  has_B_scaling = false;
1196  }
1197 
1198  EGS_Vector v1(xp->position(0),yp->position(0),zp->position(0));
1199  EGS_Vector v2(xp->position(nx),yp->position(ny),zp->position(nz));
1200  egsInformation(" top/left/front corner : (%g,%g,%g)\n",v1.x,v1.y,v1.z);
1201  egsInformation(" bottom/right/back corner: (%g,%g,%g)\n",v2.x,v2.y,v2.z);
1202  if (T) {
1203  egsInformation(" the above are transformed as follows before looking up media:\n");
1204  T->transform(v1);
1205  T->transform(v2);
1206  egsInformation(" top/left/front corner : (%g,%g,%g)\n",v1.x,v1.y,v1.z);
1207  egsInformation(" bottom/right/back corner: (%g,%g,%g)\n",v2.x,v2.y,v2.z);
1208  }
1209  for (int ix=0; ix<nx; ++ix) for (int iy=0; iy<ny; iy++) for (int iz=0; iz<nz; ++iz) {
1210  int ir = ix + iy*nx + iz*nxy;
1211  EGS_Vector v(0.5*(xp->position(ix)+xp->position(ix+1)),
1212  0.5*(yp->position(iy)+yp->position(iy+1)),
1213  0.5*(zp->position(iz)+zp->position(iz+1)));
1214  if (T) {
1215  T->transform(v);
1216  }
1217  int ireg = geometry->isWhere(v);
1218  if (ireg >= 0) {
1219  region_media[ir] = geometry->medium(ireg);
1220  if (hrs) {
1221  rhor[ir] = geometry->getRelativeRho(ir);
1222  }
1223  if (hbs) {
1224  bfactor[ir] = geometry->getBScaling(ir);
1225  }
1226  }
1227  else {
1228  region_media[ir] = -1;
1229  if (hrs) {
1230  rhor[ir] = 1;
1231  }
1232  if (hbs) {
1233  bfactor[ir] = 1;
1234  }
1235  }
1236  }
1237  vector<string> options;
1238  options.push_back("no");
1239  options.push_back("yes");
1240  bool delete_geometry = input->getInput("delete geometry",options,0);
1241  if (delete_geometry) {
1242  if (geometry->deref() == -1) {
1243  egsInformation(" deleting geometry %s as requested\n",geometry->getName().c_str());
1244  delete geometry;
1245  }
1246  else {
1247  egsInformation(" not deleting geometry %s as requested due to remaining references\n",
1248  geometry->getName().c_str());
1249  geometry->ref();
1250  }
1251  }
1252  egsInformation("\n");
1253 }
1254 
1255 const char *err_msg1 = "createGeometry(EGS_XYZRepeater)";
1256 
1257 #endif
1258 
1259 static string EGS_NDG_LOCAL typeStr("EGS_NDGeometry");
1260 string EGS_NDGeometry::type(typeStr);
1261 
1262 static bool EGS_NDG_LOCAL inputSet = false;
1263 
1264 extern "C" {
1265  static void setInputs() {
1266  inputSet = true;
1267 
1268  setBaseGeometryInputs();
1269 
1270  geomBlockInput->getSingleInput("library")->setValues({"egs_ndgeometry"});
1271 
1272  // Format: name, isRequired, description, vector string of allowed values
1273  auto typePtr = geomBlockInput->addSingleInput("type", false, "Alternative types of nd_geometry. Neglect this input to use a standard egs_nd_geometry.", {"EGS_XYZGeometry", "EGS_XYZRepeater"});
1274 
1275  auto dimPtr = geomBlockInput->addSingleInput("dimensions", true, "A list of previously defined geometries.");
1276  dimPtr->addDependency(typePtr, "", true);
1277  auto hownPtr = geomBlockInput->addSingleInput("hownear method", false, "0(for orthogonal constituent geometries) or 1");
1278  hownPtr->addDependency(typePtr, "", true);
1279 
1280  // EGS_XYZGeometry
1281  // First method
1282  auto xPtr = geomBlockInput->addSingleInput("x-planes", false, "A list of the x-plane positions");
1283  xPtr->addDependency(typePtr, "EGS_XYZGeometry");
1284  auto yPtr = geomBlockInput->addSingleInput("y-planes", false, "A list of the y-plane positions");
1285  yPtr->addDependency(typePtr, "EGS_XYZGeometry");
1286  auto zPtr = geomBlockInput->addSingleInput("z-planes", false, "A list of the z-plane positions");
1287  zPtr->addDependency(typePtr, "EGS_XYZGeometry");
1288 
1289  // Second method
1290  auto densityPtr = geomBlockInput->addSingleInput("density matrix", false, "A binary file containing the density for every region. See the documentation for format details.");
1291  densityPtr->addDependency(typePtr, "EGS_XYZGeometry");
1292  auto ctPtr = geomBlockInput->addSingleInput("ct ramp", false, "A ramp file for converting density to medium.");
1293  ctPtr->addDependency(typePtr, "EGS_XYZGeometry");
1294  auto phantPtr = geomBlockInput->addSingleInput("egsphant file", false, "An egsphant file, as defined in the DOSXYZnrc documentation.");
1295  phantPtr->addDependency(typePtr, "EGS_XYZGeometry");
1296 
1297  // For second method, must either use "density matrix" or "egsphant file"
1298  densityPtr->addDependency(phantPtr, "", true);
1299  phantPtr->addDependency(densityPtr, "", true);
1300  ctPtr->addDependency(densityPtr);
1301 
1302  // Third method
1303  auto xslabPtr = geomBlockInput->addSingleInput("x-slabs", false, "Xo Dx Nx");
1304  xslabPtr->addDependency(typePtr, "EGS_XYZGeometry");
1305  auto yslabPtr = geomBlockInput->addSingleInput("y-slabs", false, "Yo Dy Ny");
1306  yslabPtr->addDependency(typePtr, "EGS_XYZGeometry");
1307  auto zslabPtr = geomBlockInput->addSingleInput("z-slabs", false, "Zo Dz Nz");
1308  zslabPtr->addDependency(typePtr, "EGS_XYZGeometry");
1309 
1310  // Can only use one method
1311  xPtr->addDependency(densityPtr, "", true);
1312  yPtr->addDependency(densityPtr, "", true);
1313  zPtr->addDependency(densityPtr, "", true);
1314  xPtr->addDependency(ctPtr, "", true);
1315  yPtr->addDependency(ctPtr, "", true);
1316  zPtr->addDependency(ctPtr, "", true);
1317  xPtr->addDependency(phantPtr, "", true);
1318  yPtr->addDependency(phantPtr, "", true);
1319  zPtr->addDependency(phantPtr, "", true);
1320  xPtr->addDependency(xslabPtr, "", true);
1321  yPtr->addDependency(xslabPtr, "", true);
1322  zPtr->addDependency(xslabPtr, "", true);
1323  xPtr->addDependency(yslabPtr, "", true);
1324  yPtr->addDependency(yslabPtr, "", true);
1325  zPtr->addDependency(yslabPtr, "", true);
1326  xPtr->addDependency(zslabPtr, "", true);
1327  yPtr->addDependency(zslabPtr, "", true);
1328  zPtr->addDependency(zslabPtr, "", true);
1329  densityPtr->addDependency(xPtr, "", true);
1330  ctPtr->addDependency(xPtr, "", true);
1331  phantPtr->addDependency(xPtr, "", true);
1332  densityPtr->addDependency(yPtr, "", true);
1333  ctPtr->addDependency(yPtr, "", true);
1334  phantPtr->addDependency(yPtr, "", true);
1335  densityPtr->addDependency(zPtr, "", true);
1336  ctPtr->addDependency(zPtr, "", true);
1337  phantPtr->addDependency(zPtr, "", true);
1338  densityPtr->addDependency(xslabPtr, "", true);
1339  ctPtr->addDependency(xslabPtr, "", true);
1340  phantPtr->addDependency(xslabPtr, "", true);
1341  densityPtr->addDependency(yslabPtr, "", true);
1342  ctPtr->addDependency(yslabPtr, "", true);
1343  phantPtr->addDependency(yslabPtr, "", true);
1344  densityPtr->addDependency(zslabPtr, "", true);
1345  ctPtr->addDependency(zslabPtr, "", true);
1346  phantPtr->addDependency(zslabPtr, "", true);
1347  xslabPtr->addDependency(xPtr, "", true);
1348  xslabPtr->addDependency(yPtr, "", true);
1349  xslabPtr->addDependency(zPtr, "", true);
1350  xslabPtr->addDependency(densityPtr, "", true);
1351  xslabPtr->addDependency(ctPtr, "", true);
1352  xslabPtr->addDependency(phantPtr, "", true);
1353  yslabPtr->addDependency(xPtr, "", true);
1354  yslabPtr->addDependency(yPtr, "", true);
1355  yslabPtr->addDependency(zPtr, "", true);
1356  yslabPtr->addDependency(densityPtr, "", true);
1357  yslabPtr->addDependency(ctPtr, "", true);
1358  yslabPtr->addDependency(phantPtr, "", true);
1359  zslabPtr->addDependency(xPtr, "", true);
1360  zslabPtr->addDependency(yPtr, "", true);
1361  zslabPtr->addDependency(zPtr, "", true);
1362  zslabPtr->addDependency(densityPtr, "", true);
1363  zslabPtr->addDependency(ctPtr, "", true);
1364  zslabPtr->addDependency(phantPtr, "", true);
1365 
1366 
1367  // EGS_XYZRepeater
1368  auto regeomPtr = geomBlockInput->addSingleInput("repeated geometry", true, "The name of a previously defined geometry");
1369  regeomPtr->addDependency(typePtr, "EGS_XYZRepeater");
1370  auto medPtr = geomBlockInput->addSingleInput("medium", false, "The medium the space between xmin..xmax, ymin..ymax, and zmin..zmax is filled with");
1371  medPtr->addDependency(typePtr, "EGS_XYZRepeater");
1372  auto rexPtr = geomBlockInput->addSingleInput("repeat x", true, "xmin xmax Nx");
1373  rexPtr->addDependency(typePtr, "EGS_XYZRepeater");
1374  auto reyPtr = geomBlockInput->addSingleInput("repeat y", true, "ymin ymax Ny");
1375  reyPtr->addDependency(typePtr, "EGS_XYZRepeater");
1376  auto rezPtr = geomBlockInput->addSingleInput("repeat z", true, "zmin zmax Nz");
1377  rezPtr->addDependency(typePtr, "EGS_XYZRepeater");
1378  }
1379 
1380  EGS_NDG_EXPORT string getExample() {
1381  string example;
1382  example = {
1383  R"(
1384  # Example of egs_ndgeometry
1385  #:start geometry:
1386  library = EGS_NDGeometry
1387  name = my_ndgeometry
1388  dimensions = geom1 geom2
1389  :start media input:
1390  media = water
1391  :stop media input:
1392  :stop geometry:
1393 )"};
1394  return example;
1395  }
1396 
1397  EGS_NDG_EXPORT shared_ptr<EGS_BlockInput> getInputs() {
1398  if(!inputSet) {
1399  setInputs();
1400  }
1401  return geomBlockInput;
1402  }
1403 
1404  EGS_NDG_EXPORT EGS_BaseGeometry *createGeometry(EGS_Input *input) {
1405 #ifdef EXPLICIT_XYZ
1406  const static char *func = "createGeometry(XYZ)";
1407  string type;
1408  int is_xyz = input->getInput("type",type);
1409  if (!is_xyz && input->compare("EGS_XYZRepeater",type)) {
1410  string base,medium;
1411  vector<EGS_Float> xr, yr, zr;
1412  int err1 = input->getInput("repeated geometry",base);
1413  int err2 = input->getInput("repeat x",xr);
1414  int err3 = input->getInput("repeat y",yr);
1415  int err4 = input->getInput("repeat z",zr);
1416  int err5 = input->getInput("medium",medium);
1417  EGS_BaseGeometry *g = 0;
1418  if (err1) {
1419  egsWarning("%s: missing 'repeated geometry' input\n",err_msg1);
1420  }
1421  else {
1422  g = EGS_BaseGeometry::getGeometry(base);
1423  if (!g) {
1424  egsWarning("%s: no geometry named %s exists\n",err_msg1,
1425  base.c_str());
1426  err1 = 1;
1427  }
1428  }
1429  if (err2 || xr.size() != 3) {
1430  err1 = 1;
1431  egsWarning("%s: wrong/missing 'repeat x' input\n",err_msg1);
1432  }
1433  if (err3 || yr.size() != 3) {
1434  err1 = 1;
1435  egsWarning("%s: wrong/missing 'repeat y' input\n",err_msg1);
1436  }
1437  if (err4 || zr.size() != 3) {
1438  err1 = 1;
1439  egsWarning("%s: wrong/missing 'repeat z' input\n",err_msg1);
1440  }
1441  if (err1) {
1442  return 0;
1443  }
1444  EGS_Float xmin = xr[0], xmax = xr[1];
1445  int nx = (int)(xr[2]+0.1);
1446  if (xmin >= xmax || nx < 1) {
1447  egsWarning("%s: wrong 'repeat x' input xmin=%g xmax=%g nx=%d\n",
1448  xmin,xmax,nx);
1449  err1 = 1;
1450  }
1451  EGS_Float ymin = yr[0], ymax = yr[1];
1452  int ny = (int)(yr[2]+0.1);
1453  if (ymin >= ymax || ny < 1) {
1454  egsWarning("%s: wrong 'repeat y' input ymin=%g ymax=%g ny=%d\n",
1455  ymin,ymax,ny);
1456  err1 = 1;
1457  }
1458  EGS_Float zmin = zr[0], zmax = zr[1];
1459  int nz = (int)(zr[2]+0.1);
1460  if (zmin >= zmax || nz < 1) {
1461  egsWarning("%s: wrong 'repeat z' input zmin=%g zmax=%g nz=%d\n",
1462  zmin,zmax,nz);
1463  err1 = 1;
1464  }
1465  if (err1) {
1466  return 0;
1467  }
1468 
1469  EGS_XYZRepeater *result =
1470  new EGS_XYZRepeater(xmin,xmax,ymin,ymax,zmin,zmax,nx,ny,nz,g);
1471  result->setName(input);
1472  result->setBoundaryTolerance(input);
1473  if (!err5) {
1474  result->setMedium(medium);
1475  }
1476  result->setXYZLabels(input);
1477  result->setLabels(input);
1478  return result;
1479  }
1480  else if (!is_xyz && input->compare("EGS_XYZGeometry",type)) {
1481  string dens_file, ramp_file, egsphant_file, interfile_file;
1482  int ierr1 = input->getInput("density matrix",dens_file);
1483  int ierr2 = input->getInput("ct ramp",ramp_file);
1484  int ierr3 = input->getInput("egsphant file",egsphant_file);
1485  int ierr4 = input->getInput("interfile header",interfile_file);
1486  int dens_or_egsphant_or_interfile = -1;
1487  if (!ierr1) {
1488  dens_or_egsphant_or_interfile = 0;
1489  }
1490  else if (!ierr3) {
1491  dens_or_egsphant_or_interfile = 1;
1492  dens_file = egsphant_file;
1493  }
1494  else if (!ierr4) {
1495  dens_or_egsphant_or_interfile = 2;
1496  dens_file = interfile_file;
1497  }
1498  if (dens_or_egsphant_or_interfile >= 0 || !ierr2) {
1499  if (dens_or_egsphant_or_interfile < 0) {
1500  egsWarning("%s: no 'density matrix', 'egsphant file' or 'interfile header' input\n",func);
1501  return 0;
1502  }
1503  if (dens_or_egsphant_or_interfile != 1 && ierr2) {
1504  egsWarning("%s: no 'ct ramp' input\n",func);
1505  return 0;
1506  }
1507  EGS_XYZGeometry *result =
1508  EGS_XYZGeometry::constructGeometry(dens_file.c_str(),ramp_file.c_str(),dens_or_egsphant_or_interfile);
1509 
1510  if (result) {
1511  result->dens_file = dens_file;
1512  result->dens_or_egsphant_or_interfile = dens_or_egsphant_or_interfile;
1513  result->setName(input);
1514  result->setBoundaryTolerance(input);
1515  result->setBScaling(input);
1516 
1517  if (dens_or_egsphant_or_interfile == 1) {
1518  vector<string> options;
1519  options.push_back("no");
1520  options.push_back("yes");
1521  bool useEgsphantDensities = input->getInput("use egsphant densities",options,1);
1522  result->useEgsphantDensities = useEgsphantDensities;
1523  }
1524  }
1525 
1526  return result;
1527  }
1528  vector<EGS_Float> xpos, ypos, zpos, xslab, yslab, zslab;
1529  int ix = input->getInput("x-planes",xpos);
1530  int iy = input->getInput("y-planes",ypos);
1531  int iz = input->getInput("z-planes",zpos);
1532  int ix1 = input->getInput("x-slabs",xslab);
1533  int iy1 = input->getInput("y-slabs",yslab);
1534  int iz1 = input->getInput("z-slabs",zslab);
1535  int nx=0, ny=0, nz=0;
1536  if (!ix1) {
1537  if (xslab.size() != 3) {
1538  egsWarning("createGeometry(XYZ): exactly 3 inputs are required"
1539  " when using 'x-slabs' input method\n");
1540  ix1 = 1;
1541  }
1542  else {
1543  nx = (int)(xslab[2]+0.1);
1544  if (nx < 1) {
1545  egsWarning("createGeometry(XYZ): number of slabs must be"
1546  " positive!\n");
1547  ix1 = 1;
1548  }
1549  if (xslab[1] <= 0) {
1550  egsWarning("createGeometry(XYZ): slab thickness must be"
1551  " positive!\n");
1552  ix1 = 1;
1553  }
1554  }
1555  }
1556  if (ix && ix1) {
1557  egsWarning("createGeometry(XYZ): wrong/missing 'x-planes' "
1558  "and 'x-slabs' input\n");
1559  return 0;
1560  }
1561  if (!iy1) {
1562  if (yslab.size() != 3) {
1563  egsWarning("createGeometry(XYZ): exactly 3 inputs are required"
1564  " when using 'y-slabs' input method\n");
1565  iy1 = 1;
1566  }
1567  else {
1568  ny = (int)(yslab[2]+0.1);
1569  if (ny < 1) {
1570  egsWarning("createGeometry(XYZ): number of slabs must be"
1571  " positive!\n");
1572  iy1 = 1;
1573  }
1574  if (yslab[1] <= 0) {
1575  egsWarning("createGeometry(XYZ): slab thickness must be"
1576  " positive!\n");
1577  iy1 = 1;
1578  }
1579  }
1580  }
1581  if (iy && iy1) {
1582  egsWarning("createGeometry(XYZ): wrong/missing 'y-planes' "
1583  "and 'y-slabs' input\n");
1584  return 0;
1585  }
1586  if (!iz1) {
1587  if (zslab.size() != 3) {
1588  egsWarning("createGeometry(XYZ): exactly 3 inputs are required"
1589  " when using 'z-slabs' input method\n");
1590  iz1 = 1;
1591  }
1592  else {
1593  nz = (int)(zslab[2]+0.1);
1594  if (nz < 1) {
1595  egsWarning("createGeometry(XYZ): number of slabs must be"
1596  " positive!\n");
1597  iz1 = 1;
1598  }
1599  if (zslab[1] <= 0) {
1600  egsWarning("createGeometry(XYZ): slab thickness must be"
1601  " positive!\n");
1602  iz1 = 1;
1603  }
1604  }
1605  }
1606  if (iz && iz1) {
1607  egsWarning("createGeometry(XYZ): wrong/missing 'z-planes' "
1608  "and 'z-slabs' input\n");
1609  return 0;
1610  }
1611  EGS_PlanesX *xp = !ix1 ?
1612  new EGS_PlanesX(xslab[0],xslab[1],nx,"",EGS_XProjector("x-planes")) :
1613  new EGS_PlanesX(xpos,"",EGS_XProjector("x-planes"));
1614  EGS_PlanesY *yp = !iy1 ?
1615  new EGS_PlanesY(yslab[0],yslab[1],ny,"",EGS_YProjector("y-planes")) :
1616  new EGS_PlanesY(ypos,"",EGS_YProjector("y-planes"));
1617  EGS_PlanesZ *zp = !iz1 ?
1618  new EGS_PlanesZ(zslab[0],zslab[1],nz,"",EGS_ZProjector("z-planes")) :
1619  new EGS_PlanesZ(zpos,"",EGS_ZProjector("z-planes"));
1620  EGS_XYZGeometry *result = new EGS_XYZGeometry(xp,yp,zp);
1621 
1622  if (result) {
1623  egsWarning("**********************************************\n");
1624  EGS_BaseGeometry *g = result;
1625  result->setName(input);
1626  result->setBoundaryTolerance(input);
1627  g->setMedia(input);
1628  result->voxelizeGeometry(input);
1629 
1630  // labels
1631  result->setXYZLabels(input);
1632  g->setLabels(input);
1633  }
1634  return result;
1635  }
1636 #endif
1637  vector<EGS_BaseGeometry *> dims;
1638  EGS_Input *ij;
1639  int error = 0;
1640  while ((ij = input->takeInputItem("geometry",false)) != 0) {
1641  EGS_BaseGeometry *g = EGS_BaseGeometry::createSingleGeometry(ij);
1642  if (g) {
1643  dims.push_back(g);
1644  g->ref();
1645  }
1646  else {
1647  error++;
1648  }
1649  delete ij;
1650  }
1651  vector<string> gnames;
1652  int err1 = input->getInput("dimensions",gnames);
1653  if (!err1) {
1654  for (unsigned int j=0; j<gnames.size(); j++) {
1655  EGS_BaseGeometry *g = EGS_BaseGeometry::getGeometry(gnames[j]);
1656  if (g) {
1657  dims.push_back(g);
1658  g->ref();
1659  }
1660  else {
1661  egsWarning("Geometry %s does not exist\n",gnames[j].c_str());
1662  error++;
1663  }
1664  }
1665  }
1666  if (error) {
1667  egsWarning("createGeometry(ND_Geometry): %d errors while "
1668  "creating/geting geometries defining individual dimensions\n",error);
1669  return 0;
1670  }
1671  if (dims.size() < 2) {
1672  egsWarning("createGeometry(ND_Geometry): why do you want to "
1673  "construct a ND geometry with a single dimension?\n");
1674  input->print(0,cerr);
1675  for (int j=0; j<gnames.size(); j++) egsWarning("dimension %d: %s\n",
1676  j+1,gnames[j].c_str());
1677  }
1678  int n_concav = 0;
1679  for (int j=0; j<dims.size(); j++) if (!dims[j]->isConvex()) {
1680  n_concav++;
1681  }
1682  bool is_ok = true;
1683  if (n_concav > 1) {
1684  egsWarning("createGeometry(ND_Geometry): a ND geometry can not have "
1685  " more than one non-convex dimension, yours has %d\n",n_concav);
1686  is_ok = false;
1687  }
1688  else if (n_concav == 1) {
1689  if (dims[dims.size()-1]->isConvex()) {
1690  egsWarning("createGeometry(ND_Geometry): the non-convex "
1691  "dimension must be the last dimension\n");
1692  is_ok = false;
1693  }
1694  }
1695  if (!is_ok) {
1696  for (int j=0; j<dims.size(); j++)
1697  if (dims[j]->deref()) {
1698  delete dims[j];
1699  }
1700  return 0;
1701  }
1702  int hn_method=0;
1703  err1 = input->getInput("hownear method",hn_method);
1704  EGS_BaseGeometry *result;
1705  if (!err1 && hn_method == 1) {
1706  result = new EGS_NDGeometry(dims,"",false);
1707  }
1708  else {
1709  result = new EGS_NDGeometry(dims);
1710  }
1711  result->setName(input);
1712  result->setBoundaryTolerance(input);
1713  result->setMedia(input);
1714  result->setLabels(input);
1715  return result;
1716  }
1717 
1718 
1719  void EGS_XYZGeometry::setXYZLabels(EGS_Input *input) {
1720 
1721  // x,y,z labels
1722  string inp;
1723  int err;
1724 
1725  err = input->getInput("set x label", inp);
1726  if (!err) {
1727  xp->setLabels(inp);
1728  }
1729 
1730  err = input->getInput("set y label", inp);
1731  if (!err) {
1732  yp->setLabels(inp);
1733  }
1734 
1735  err = input->getInput("set z label", inp);
1736  if (!err) {
1737  zp->setLabels(inp);
1738  }
1739  }
1740 
1741 
1742  void EGS_NDGeometry::ndRegions(int r, int dim, int dimk, int k, vector<int> &regs) {
1743 
1744  // skip looping over selected dimension
1745  if (dim == dimk) {
1746  r += k*n[dimk];
1747  if (dim == N-1) {
1748  regs.push_back(r);
1749  }
1750  else {
1751  ndRegions(r, dim+1, dimk, k, regs);
1752  }
1753  }
1754 
1755  // last dimension: end recursion and record global region number
1756  else if (dim == N-1) {
1757  for (int j=0; j<g[dim]->regions(); j++) {
1758  regs.push_back(r+j*n[dim]);
1759  }
1760  }
1761 
1762  // keep collecting by recursion
1763  else {
1764  for (int j=0; j<g[dim]->regions(); j++) {
1765  ndRegions(r + j*n[dim], dim+1, dimk, k, regs);
1766  }
1767  }
1768  }
1769 
1770 
1771  void EGS_NDGeometry::getLabelRegions(const string &str, vector<int> &regs, bool sanitize) {
1772 
1773  // label defined in the sub-geometries
1774  vector<int> local_regs;
1775  for (int i=0; i<N; i++) {
1776  local_regs.clear();
1777  if (g[i]) {
1778  g[i]->getLabelRegions(str, local_regs, sanitize);
1779  }
1780  if (local_regs.size() == 0) {
1781  continue;
1782  }
1783 
1784  // recurse to collect all global regions comprised in each local region
1785  for (int j=0; j<local_regs.size(); j++) {
1786  ndRegions(0, 0, i, local_regs[j], regs);
1787  }
1788  }
1789 
1790  // label defined in self (nd input block)
1791  EGS_BaseGeometry::getLabelRegions(str, regs, sanitize);
1792 
1793  }
1794 
1795 
1796  void EGS_XYZGeometry::getLabelRegions(const string &str, vector<int> &regs, bool sanitize) {
1797 
1798  vector<int> local_regs;
1799 
1800  // x plane labels
1801  local_regs.clear();
1802  xp->getLabelRegions(str, local_regs, sanitize);
1803  for (int i=0; i<local_regs.size(); i++) {
1804  for (int j=0; j<ny; j++) {
1805  for (int k=0; k<nz; k++) {
1806  regs.push_back(local_regs[i] + nx*j + nxy*k);
1807  }
1808  }
1809  }
1810 
1811  // y plane labels
1812  local_regs.clear();
1813  yp->getLabelRegions(str, local_regs, sanitize);
1814  for (int j=0; j<local_regs.size(); j++) {
1815  for (int i=0; i<nx; i++) {
1816  for (int k=0; k<nz; k++) {
1817  regs.push_back(i + nx*local_regs[j] + nxy*k);
1818  }
1819  }
1820  }
1821 
1822  // z plane labels
1823  local_regs.clear();
1824  zp->getLabelRegions(str, local_regs, sanitize);
1825  for (int k=0; k<local_regs.size(); k++) {
1826  for (int i=0; i<nx; i++) {
1827  for (int j=0; j<ny; j++) {
1828  regs.push_back(i + nx*j + nxy*local_regs[k]);
1829  }
1830  }
1831  }
1832 
1833  // label defined in self (xyz input block)
1834  EGS_BaseGeometry::getLabelRegions(str, regs, sanitize);
1835 
1836  }
1837 
1838 
1839  void EGS_XYZRepeater::getLabelRegions(const string &str, vector<int> &regs, bool sanitize) {
1840 
1841  vector<int> local_regs;
1842 
1843  // label in repeated geometry
1844  local_regs.clear();
1845  g->getLabelRegions(str, local_regs, sanitize);
1846  for (int i=0; i<nx; i++) {
1847  for (int j=0; j<ny; j++) {
1848  for (int k=0; k<nz; k++) {
1849  for (int r=0; r<local_regs.size(); r++) {
1850  regs.push_back(ng*(i + nx*j + nxy*k) + local_regs[r]);
1851  }
1852  }
1853  }
1854  }
1855 
1856  // x plane labels
1857  local_regs.clear();
1858  xyz->getXLabelRegions(str, local_regs, sanitize);
1859  for (int i=0; i<local_regs.size(); i++) {
1860  for (int j=0; j<ny; j++) {
1861  for (int k=0; k<nz; k++) {
1862  for (int r=0; r<ng; r++) {
1863  regs.push_back(ng*(local_regs[i] + nx*j + nxy*k) + r);
1864  }
1865  }
1866  }
1867  }
1868 
1869  // y plane labels
1870  local_regs.clear();
1871  xyz->getYLabelRegions(str, local_regs, sanitize);
1872  for (int j=0; j<local_regs.size(); j++) {
1873  for (int i=0; i<nx; i++) {
1874  for (int k=0; k<nz; k++) {
1875  for (int r=0; r<ng; r++) {
1876  regs.push_back(ng*(i + nx*local_regs[j] + nxy*k) + r);
1877  }
1878  }
1879  }
1880  }
1881 
1882  // z plane labels
1883  local_regs.clear();
1884  xyz->getZLabelRegions(str, local_regs, sanitize);
1885  for (int k=0; k<local_regs.size(); k++) {
1886  for (int i=0; i<nx; i++) {
1887  for (int j=0; j<ny; j++) {
1888  for (int r=0; r<ng; r++) {
1889  regs.push_back(ng*(i + nx*j + nxy*local_regs[k]) + r);
1890  }
1891  }
1892  }
1893  }
1894 
1895  // label defined in self (repeater input block)
1896  EGS_BaseGeometry::getLabelRegions(str, regs, sanitize);
1897 
1898  }
1899 
1900 }
EGS_AffineTransform
A class providing affine transformations.
Definition: egs_transformations.h:432
EGS_AffineTransform::getTransformation
static EGS_AffineTransform * getTransformation(EGS_Input *inp)
Constructs an affine transformation object from the input pointed to by inp and returns a pointer to ...
Definition: egs_transformations.cpp:44
EGS_AffineTransform::transform
void transform(EGS_Vector &v) const
Transforms the vector v.
Definition: egs_transformations.h:573
EGS_BaseGeometry
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
Definition: egs_base_geometry.h:111
EGS_BaseGeometry::deref
int deref()
Decrease the reference count to this geometry.
Definition: egs_base_geometry.h:729
EGS_BaseGeometry::getBScaling
virtual EGS_Float getBScaling(int ireg) const
Get the B field scaling factor in region ireg.
Definition: egs_base_geometry.h:481
EGS_BaseGeometry::setBScaling
virtual void setBScaling(int start, int end, EGS_Float bf)
Set the B field scaling factor in regions.
Definition: egs_base_geometry.cpp:865
EGS_BaseGeometry::setMedia
void setMedia(EGS_Input *inp)
Set the media in the geometry from the input pointed to by inp.
Definition: egs_base_geometry.cpp:747
EGS_BaseGeometry::createSingleGeometry
static EGS_BaseGeometry * createSingleGeometry(EGS_Input *inp)
Create a single geometry from the input inp.
Definition: egs_base_geometry.cpp:547
EGS_BaseGeometry::nreg
int nreg
Number of local regions in this geometry.
Definition: egs_base_geometry.h:818
EGS_BaseGeometry::has_B_scaling
bool has_B_scaling
Does this geometry has B field scaling factor?
Definition: egs_base_geometry.h:866
EGS_BaseGeometry::is_convex
bool is_convex
Is this geometry convex?
Definition: egs_base_geometry.h:908
EGS_BaseGeometry::has_rho_scaling
bool has_rho_scaling
Does this geometry have relative mass density scvaling?
Definition: egs_base_geometry.h:856
EGS_BaseGeometry::hasBScaling
bool hasBScaling() const
Does this geometry object have a B field scaling feature?
Definition: egs_base_geometry.h:475
EGS_BaseGeometry::region_media
short * region_media
Array of media indeces.
Definition: egs_base_geometry.h:842
EGS_BaseGeometry::hasRhoScaling
virtual bool hasRhoScaling()
Does this geometry object have a mass density scaling feature?
Definition: egs_base_geometry.h:440
EGS_BaseGeometry::setName
void setName(EGS_Input *inp)
Set the name of the geometry from the input inp.
Definition: egs_base_geometry.cpp:613
EGS_BaseGeometry::setMedium
void setMedium(const string &Name)
Set all regions to a medium with name Name.
Definition: egs_base_geometry.cpp:484
EGS_BaseGeometry::getName
const string & getName() const
Get the name of this geometry.
Definition: egs_base_geometry.h:535
EGS_BaseGeometry::rhor
EGS_Float * rhor
Array with relative mass densities.
Definition: egs_base_geometry.h:861
EGS_BaseGeometry::getRelativeRho
virtual EGS_Float getRelativeRho(int ireg) const
Get the relative mass density in region ireg.
Definition: egs_base_geometry.h:447
EGS_BaseGeometry::bfactor
EGS_Float * bfactor
Array with B field scaling factors.
Definition: egs_base_geometry.h:871
EGS_BaseGeometry::isConvex
bool isConvex() const
Is the geometry convex?
Definition: egs_base_geometry.h:142
EGS_BaseGeometry::med
int med
Medium index.
Definition: egs_base_geometry.h:850
EGS_BaseGeometry::medium
virtual int medium(int ireg) const
Returns the medium index in region ireg.
Definition: egs_base_geometry.h:327
EGS_BaseGeometry::name
string name
Name of this geometry.
Definition: egs_base_geometry.h:834
EGS_BaseGeometry::regions
int regions() const
Returns the number of local regions in this geometry.
Definition: egs_base_geometry.h:306
EGS_BaseGeometry::setLabels
int setLabels(EGS_Input *input)
Set the labels from an input block.
Definition: egs_base_geometry.cpp:1127
EGS_BaseGeometry::setRelativeRho
virtual void setRelativeRho(int start, int end, EGS_Float rho)
Set the relative mass density in regions.
Definition: egs_base_geometry.cpp:817
EGS_BaseGeometry::printInfo
virtual void printInfo() const
Print information about this geometry.
Definition: egs_base_geometry.cpp:730
EGS_BaseGeometry::ref
int ref()
Increase the reference count to this geometry.
Definition: egs_base_geometry.h:717
EGS_BaseGeometry::getGeometry
static EGS_BaseGeometry * getGeometry(const string &Name)
Get a pointer to the geometry named Name.
Definition: egs_base_geometry.cpp:472
EGS_BaseGeometry::addMedium
static int addMedium(const string &medname)
Add a medium or get the index of an existing medium.
Definition: egs_base_geometry.cpp:492
EGS_BaseGeometry::getMediumName
static const char * getMediumName(int ind)
Get the name of medium with index ind.
Definition: egs_base_geometry.cpp:535
EGS_BaseGeometry::setBoundaryTolerance
void setBoundaryTolerance(EGS_Input *inp)
Set the value of the boundary tolerance from the input inp.
Definition: egs_base_geometry.cpp:721
EGS_BaseGeometry::getLabelRegions
virtual void getLabelRegions(const string &str, vector< int > &regs, bool sanitize=true)
Get the list of all regions labeled with str.
Definition: egs_base_geometry.cpp:1077
EGS_BaseGeometry::isWhere
virtual int isWhere(const EGS_Vector &x)=0
In which region is poisition x?
EGS_DeformedXYZ::setDeformations
int setDeformations(const char *defFile)
Definition: egs_nd_geometry.cpp:1008
EGS_Input
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
EGS_Input * takeInputItem(const string &key, bool self=true)
Get the property named key.
Definition: egs_input.cpp:229
EGS_Input::print
void print(int nind, ostream &)
Used for debugging purposes.
Definition: egs_input.cpp:1177
EGS_Input::compare
static bool compare(const string &s1, const string &s2)
Definition: egs_input.cpp:1173
EGS_Input::getInput
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:341
EGS_NDGeometry
A class modeling a N-dimensional geometry.
Definition: egs_nd_geometry.h:315
EGS_NDGeometry::type
static string type
The geometry type.
Definition: egs_nd_geometry.h:626
EGS_NDGeometry::EGS_NDGeometry
EGS_NDGeometry(int ng, EGS_BaseGeometry **G, const string &Name="", bool O=true)
Definition: egs_nd_geometry.cpp:59
EGS_NDGeometry::n
int * n
Used for calculating region indeces.
Definition: egs_nd_geometry.h:625
EGS_NDGeometry::setMedia
void setMedia(EGS_Input *inp, int nmed, const int *med_ind)
Define media.
Definition: egs_nd_geometry.cpp:115
EGS_NDGeometry::g
EGS_BaseGeometry ** g
The dimensions.
Definition: egs_nd_geometry.h:624
EGS_NDGeometry::N
int N
Number of dimensions.
Definition: egs_nd_geometry.h:623
EGS_PlanesT
A set of parallel planes.
Definition: egs_planes.h:167
EGS_Vector
A class representing 3D vectors.
Definition: egs_vector.h:57
EGS_XProjector
A projector into the x-plane.
Definition: egs_projectors.h:119
EGS_XYZGeometry
An XYZ-geometry.
Definition: egs_nd_geometry.h:794
EGS_XYZRepeater
A geometry repeated on a regular XYZ grid.
Definition: egs_nd_geometry.h:1761
EGS_YProjector
A projector into the y-plane.
Definition: egs_projectors.h:161
EGS_ZProjector
A projector into the z-plane.
Definition: egs_projectors.h:203
egs_functions.h
Global egspp functions header file.
createGeometry
EGS_GLIB_EXPORT EGS_BaseGeometry * createGeometry(EGS_Input *input)
Definition: egs_glib.cpp:84
egs_input.h
EGS_Input class header file.
egs_nd_geometry.h
N-dimensional geometries: header.
egs_transformations.h
EGS_AffineTransform and EGS_RotationMatrix class header file.
egsGetEndian
int egsGetEndian()
Get the endianess of the machine.
Definition: egs_functions.cpp:369
egsInformation
EGS_InfoFunction EGS_EXPORT egsInformation
Always use this function for reporting the progress of a simulation and any other type of information...
Definition: egs_functions.cpp:135
egsFatal
EGS_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
Definition: egs_functions.cpp:137
loopMax
const EGS_I64 loopMax
The maximum number of iterations for near-infinite loops.
Definition: egs_functions.h:96
epsilon
const EGS_Float epsilon
The epsilon constant for floating point comparisons.
Definition: egs_functions.h:62
egsSwapBytes
void egsSwapBytes(int *n)
Swap the bytes of 32 bit integers.
Definition: egs_functions.cpp:181
egsWarning
EGS_InfoFunction EGS_EXPORT egsWarning
Always use this function for reporting warnings.
Definition: egs_functions.cpp:136
Generated by   doxygen 1.9.1