EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_advanced_application.cpp
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ advanced application
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: Ernesto Mainegra-Hing
27 # Frederic Tessier
28 # Blake Walters
29 # Randle Taylor
30 # Reid Townson
31 # Max Orok
32 # Martin Martinov
33 #
34 ###############################################################################
35 */
36 
37 
47 #include "egs_functions.h"
48 #include "egs_interface2.h"
49 #include "egs_run_control.h"
50 #include "egs_base_source.h"
51 #include "egs_input.h"
52 #include "egs_interpolator.h"
53 #include "egs_rndm.h"
54 #include "egs_math.h"
55 #include "egs_ausgab_object.h"
56 
57 #include <string>
58 #include <cstdio>
59 #include <cstdarg>
60 #include <cstdlib>
61 #include <cstring>
62 
63 using namespace std;
64 
65 #ifdef DEBUG_APPLICATION
66 #define CHECK_GET_APPLICATION(a,b) \
67  EGS_Application *a = EGS_Application::activeApplication(); \
68  if( !a ) egsFatal("%s: no active application!\n",b);
69 #else
70 #define CHECK_GET_APPLICATION(a,b) \
71  EGS_Application *a = EGS_Application::activeApplication();
72 #endif
73 
75 
76 /*
77 extern __extc__ struct EGS_Stack F77_OBJ(stack,STACK);
78 extern __extc__ struct EGS_Bounds F77_OBJ(bounds,BOUNDS);
79 extern __extc__ struct EGS_Thresh F77_OBJ(thresh,THRESH);
80 extern __extc__ struct EGS_Epcont F77_OBJ(epcont,EPCONT);
81 extern __extc__ struct EGS_EtControl F77_OBJ_(et_control,ET_CONTROL);
82 extern __extc__ struct EGS_Media F77_OBJ(media,MEDIA);
83 extern __extc__ struct EGS_Useful F77_OBJ(useful,USEFUL);
84 extern __extc__ struct EGS_XOptions F77_OBJ_(xsection_options,XSECTION_OPTIONS);
85 extern __extc__ struct EGS_IO F77_OBJ_(egs_io,EGS_IO);
86 extern __extc__ struct EGS_VarianceReduction F77_OBJ_(egs_vr,EGS_VR);
87 
88 struct EGS_Stack *the_stack = & F77_OBJ(stack,STACK);
89 struct EGS_Bounds *the_bounds = & F77_OBJ(bounds,BOUNDS);
90 struct EGS_Thresh *the_thresh = & F77_OBJ(thresh,THRESH);
91 struct EGS_Epcont *the_epcont = & F77_OBJ(epcont,EPCONT);
92 struct EGS_EtControl *the_etcontrol = & F77_OBJ_(et_control,ET_CONTROL);
93 struct EGS_Media *the_media = & F77_OBJ(media,MEDIA);
94 struct EGS_Useful *the_useful = & F77_OBJ(useful,USEFUL);
95 struct EGS_XOptions *the_xoptions = & F77_OBJ_(xsection_options,XSECTION_OPTIONS);
96 struct EGS_IO *the_egsio = & F77_OBJ_(egs_io,EGS_IO);
97 struct EGS_VarianceReduction *the_egsvr = & F77_OBJ_(egs_vr,EGS_VR);
98 */
99 
100 #define egsGetRNGPointers F77_OBJ_(egs_get_rng_pointers,EGS_GET_RNG_POINTERS)
101 extern __extc__ void egsGetRNGPointers(EGS_I32 *, EGS_I32 *);
102 #define egsGetRNGArray F77_OBJ_(egs_get_rng_array,EGS_GET_RNG_ARRAY)
103 extern __extc__ void egsGetRNGArray(EGS_Float *);
104 #define egsSetRNGState F77_OBJ_(egs_set_rng_state,EGS_SET_RNG_STATE)
105 extern __extc__ void egsSetRNGState(const EGS_I32 *, const EGS_Float *);
106 #define egsGetSteps F77_OBJ_(egs_get_steps,EGS_GET_STEPS)
107 extern __extc__ void egsGetSteps(double *, double *);
108 #define egsSetSteps F77_OBJ_(egs_set_steps,EGS_SET_STEPS)
109 extern __extc__ void egsSetSteps(const double *, const double *);
110 #define egsOpenUnits F77_OBJ_(egs_open_units,EGS_OPEN_UNITS)
111 extern __extc__ void egsOpenUnits(const EGS_I32 *);
112 #define egsGetElectronData F77_OBJ_(egs_get_electron_data,EGS_GET_ELECTRON_DATA)
113 extern __extc__ void egsGetElectronData(void (*func)(EGS_I32 *,EGS_Float *,
114  EGS_Float *,EGS_Float *,EGS_Float *),const EGS_I32 *,const EGS_I32 *);
115 #define egsGetPhotonData F77_OBJ_(egs_get_photon_data,EGS_GET_PHOTON_DATA)
116 extern __extc__ void egsGetPhotonData(void (*func)(EGS_I32 *,EGS_Float *,
117  EGS_Float *,EGS_Float *,EGS_Float *),const EGS_I32 *,const EGS_I32 *);
118 
119 static EGS_Float *__help1, *__help2, *__help3, *__help4;
120 static EGS_I32 *__ihelp;
121 static void __help_get_data(EGS_I32 *nbin,EGS_Float *emin, EGS_Float *emax,
122  EGS_Float *a, EGS_Float *b) {
123  __ihelp = nbin;
124  __help1 = emin;
125  __help2 = emax;
126  __help3 = a;
127  __help4 = b;
128 }
129 
130 static void __init_interpolator(int q, int imed, int type,
131  EGS_Interpolator &i) {
132  if (q) {
133  egsGetElectronData(__help_get_data,&imed,&type);
134  }
135  else {
136  egsGetPhotonData(__help_get_data,&imed,&type);
137  }
138  i.initialize(*__ihelp,*__help1,*__help2,__help3,__help4);
139 }
140 
142  \
143  EGS_Application::setAusgabCall(call,on_or_off);
144  the_epcont->iausfl[call] = (int) on_or_off;
145 }
146 
147 static EGS_LOCAL void __egs_iovar(int nio, int ns,
148  const char *source, char *var) {
149  if (ns > nio-1) {
150  egsFatal("Mortran variable is not long enough to hold %s\n",source);
151  }
152  int j;
153  for (j=0; j<ns; j++) {
154  var[j] = source[j];
155  }
156  var[ns] = 0;
157  for (j=ns+1; j<nio; j++) {
158  var[j] = ' ';
159  }
160 }
161 
162 static EGS_LOCAL char __write_buf[8192];
163 extern __extc__ void F77_OBJ_(egs_write_string,EGS_WRITE_STRING)(int *,
164  const char *,int);
165 void EGS_LOCAL __write_to_fortran_file(const char *msg, ...) {
166  va_list ap;
167  va_start(ap, msg);
168  vsprintf(__write_buf,msg,ap);
169  va_end(ap);
170  int ounit=the_egsio->i_log;
171  F77_OBJ_(egs_write_string,EGS_WRITE_STRING)(&ounit,__write_buf,
172  strlen(__write_buf));
173 }
174 void EGS_LOCAL __write_to_fortran_file_and_exit(const char *msg, ...) {
175  va_list ap;
176  va_start(ap, msg);
177  vsprintf(__write_buf,msg,ap);
178  va_end(ap);
179  int ounit=the_egsio->i_log;
180  F77_OBJ_(egs_write_string,EGS_WRITE_STRING)(&ounit,__write_buf,
181  strlen(__write_buf));
182  exit(1);
183 }
184 
185 extern __extc__ void F77_OBJ_(egs_set_defaults,EGS_SET_DEFAULTS)();
186 extern __extc__ void F77_OBJ_(egs_init1,EGS_INIT1)();
188  F77_OBJ_(egs_set_defaults,EGS_SET_DEFAULTS)();
189  __egs_iovar(64,app_name.size(),app_name.c_str(),the_egsio->user_code);
190  __egs_iovar(128,hen_house.size(),hen_house.c_str(),the_egsio->hen_house);
191  __egs_iovar(128,egs_home.size(),egs_home.c_str(),the_egsio->egs_home);
192  __egs_iovar(256,input_file.size(),input_file.c_str(),the_egsio->input_file);
193  __egs_iovar(256,final_output_file.size(),final_output_file.c_str(),
194  the_egsio->output_file);
195  __egs_iovar(256,pegs_file.size(),pegs_file.c_str(),the_egsio->pegs_file);
196  the_egsio->i_parallel = i_parallel;
197  the_egsio->n_parallel = n_parallel;
198  the_egsio->is_pegsless = is_pegsless;
199  if (run) {
200  the_egsio->n_chunk = run->getNchunk();
201  }
202  the_egsio->is_batch = (int) batch_run;
203  F77_OBJ_(egs_init1,EGS_INIT1)();
204  /*
205  if( batch_run ) {
206  egsSetInfoFunction(Information,__write_to_fortran_file);
207  egsSetInfoFunction(Warning,__write_to_fortran_file);
208  egsSetInfoFunction(Fatal,__write_to_fortran_file_and_exit);
209  }
210  */
211  io_flag = 1;
212  return 0;
213 }
214 
215 /* The following 2 functions were used in an attempt to redirect
216  all fortran I/O to use egsInformationn a degsWarning.
217  Unfortunately, there were segmentation violations in the
218  fortran I/O functions that I don't understand.
219 
220 extern __extc__ void F77_OBJ_(egs_info_output,EGS_INFO_OUTPUT)
221  (char *msg, int len) {
222  int j; for(j=len-1; j>=0; j--) if( !isspace(msg[j]) ) break;
223  if( j<len-1 ) ++j;
224  char c = msg[j]; msg[j] = 0;
225  egsInformation("%s\n",msg);
226  msg[j] = c;
227 }
228 
229 extern __extc__ void F77_OBJ_(egs_warning_output,EGS_WARNING_OUTPUT)
230  (char *msg, int len) {
231  int j; for(j=len-1; j>=0; j--) if( !isspace(msg[j]) ) break;
232  if( j<len-1 ) ++j;
233  char c = msg[j]; msg[j] = 0;
234  egsWarning("%s\n",msg);
235  msg[j] = c;
236 }
237 */
238 
239 #ifndef SKIP_DOXYGEN
247 class EGS_LOCAL EGS_TransportProperty {
248 
249  string name;
250  vector<string> options;
251  vector<string> *str_v;
252  vector<EGS_Float> *f_v;
253  string charinp;
254  int type;
255  EGS_I32 *int_input;
256  EGS_Float *float_input;
257  int len,nitem;
258  char *char_input;
259 
260 public:
261 
262  /* Integer input type = 0 */
263  EGS_TransportProperty(const char *Name, EGS_I32 *val) :
264  name(Name), type(0), int_input(val) {};
265  /* Real input type = 1 */
266  EGS_TransportProperty(const char *Name, EGS_Float *val) :
267  name(Name), type(1), float_input(val) {};
268  /* Character string input type = 3 */
269  EGS_TransportProperty(const char *Name, int L, char *val) :
270  name(Name), type(3), len(L), char_input(val) {};
271  /* Character string array input type = 4 */
272  EGS_TransportProperty(const char *Name,int L,int N,vector<string> *str):
273  name(Name), str_v(str), type(4), len(L), nitem(N) {};
274  /* Real array input type = 5 */
275  EGS_TransportProperty(const char *Name,int N,vector<EGS_Float> *f):
276  name(Name), f_v(f), type(5), nitem(N) {};
277  /* Integer input with allowed values => type = 2 */
278  void addOption(const char *opt) {
279  options.push_back(opt);
280  type = 2;
281  };
282  void setOption(const int &index, const char *opt) {
283  options[index] = opt;
284  };
285  void getInput(EGS_Input *input) {
286  if (type == 1) {
287  EGS_Float aux;
288  int err = input->getInput(name,aux);
289  if (!err) {
290  *float_input = aux;
291  }
292  }
293  else if (type == 0) {
294  EGS_I32 aux;
295  int err = input->getInput(name,aux);
296  if (!err) {
297  *int_input = aux;
298  }
299  }
300  else if (type == 3) {
301  int err = input->getInput(name,charinp);
302  if (!err) {
303  int l = charinp.size();
304  if (l > len) {
305  l = len;
306  }
307  int j;
308  for (j=0; j<l; j++) {
309  char_input[j] = charinp[j];
310  }
311  for (j=l; j<len; j++) {
312  char_input[j] = ' ';
313  }
314  }
315  else {
316  charinp = string(char_input).substr(0,len);
317  }
318  }
319  else if (type == 4) {
320  int err = input->getInput(name,*str_v);
321  if (!err) {
322  vector<string> str = *str_v;
323  str_v->clear();
324  int number = str.size();
325  if (number>nitem) {
326  str.erase(str.begin()+nitem,str.end());
327  number=nitem;
328  }
329  for (int i=0; i<number; i++) {
330  if (str[i].size()>len) {
331  str[i] = str[i].substr(0,len);
332  }
333  }
334  *str_v = str;
335  }
336  }
337  else if (type == 5) {
338  int err = input->getInput(name,*f_v);
339  if (!err) {
340  vector<EGS_Float> fv = *f_v;
341  f_v->clear();
342  int number = fv.size();
343  if (number>nitem) {
344  fv.erase(fv.begin()+nitem,fv.end());
345  number=nitem;
346  }
347  *f_v = fv;
348  }
349  }
350  else {
351  bool is_ok;
352  EGS_I32 iaux = input->getInput(name,options,0,&is_ok);
353  if (is_ok) {
354  *int_input = iaux;
355  }
356  egsWarning("Property %s: %d\n",name.c_str(),*int_input);
357  }
358  };
359  void info(int nc) {
360  egsInformation("%s",name.c_str());
361  for (int j=name.size(); j<nc; j++) {
362  egsInformation(" ");
363  }
364  if (type == 0) {
365  egsInformation("%d\n",*int_input);
366  }
367  else if (type == 1) {
368  egsInformation("%g\n",*float_input);
369  }
370  else if (type == 3) {
371  egsInformation("%s\n",(string(char_input).substr(0,len)).c_str());
372  }
373  else if (type == 5) {
374  for (vector<EGS_Float>::iterator it = f_v->begin() ; it != f_v->end(); ++it) {
375  egsInformation("%g ",*it);
376  }
377  egsInformation("\n");
378  }
379  else {
380  egsInformation("%s\n",options[*int_input].c_str());
381  }
382  };
383  int size() {
384  if (type == 1) {
385  return sizeof(EGS_Float);
386  }
387  else if (type == 3) {
388  return len;
389  }
390  else if (type == 4) {
391  return nitem;
392  }
393  else if (type == 5) {
394  return f_v->size();
395  }
396  else {
397  return sizeof(EGS_I32);
398  }
399  };
400 
401 };
402 #endif
403 
405  EGS_Application(argc,argv), nmed(0), n_rng_buffer(0), final_job(false), io_flag(0) { }
406 
408  if (n_rng_buffer > 0) {
409  delete [] rng_buffer;
410  }
411  if (nmed > 0) {
412  delete [] i_ededx;
413  delete [] i_pdedx;
414  delete [] i_esig;
415  delete [] i_psig;
416  delete [] i_ebr1;
417  delete [] i_pbr1;
418  delete [] i_pbr2;
419  delete [] i_gmfp;
420  delete [] i_gbr1;
421  delete [] i_gbr2;
422  delete [] i_cohe;
423  delete [] i_photonuc;
424  }
425 }
426 
429  if (final_job) {
430  helpInit(0,false);
431  }
432 }
433 
435  egsWarning("In initCrossSections(): spin effects = %d\n",
437  EGS_Input *transportp = input->getInputItem("MC transport parameter");
438  if (!transportp) {
439  transportp = input->getInputItem("transport parameter");
440  }
441  return helpInit(transportp,true);
442 }
443 
444 extern "C" void F77_OBJ_(set_elastic_parameter,SET_ELASTIC_PARAMETER)();
445 
446 int EGS_AdvancedApplication::helpInit(EGS_Input *transportp, bool do_hatch) {
447  if (!geometry) {
448  egsWarning("initCrossSections(): no geometry?\n");
449  return 1;
450  }
452  int nmed_new = EGS_BaseGeometry::nMedia(), j;
453  if (nmed_new < 1) {
454  egsWarning("initCrossSections(): no media in this geometry?\n");
455  return 2;
456  }
457  int nmed_old = nmed;
458  nmed = nmed_new;
459  int *ind = new int [nmed];
460  for (j=0; j<nmed; j++) {
461  const char *medname = geometry->getMediumName(j);
462  int len = strlen(medname);
463  ind[j] = egsAddMedium(medname, len);
464  }
465 
466  vector<string> ff_media, ff_names;
467  string eii_xsect;
468  vector<EGS_Float> efield_v, bfield_v;
469  //
470  // The Intel compiler uses -1 for .true. => all logical
471  // variables don't work!
472  // This is a fix for this
473  //
474  if (the_xoptions->spin_effects != 0) {
476  }
477  //if( the_xoptions->exact_bca != 0 ) the_xoptions->exact_bca=1;
478 
479  //
480  // See if there are transport parameters set in the input file.
481  //
482  EGS_TransportProperty efield("Electric Field",3,&efield_v);
483  EGS_TransportProperty bfield("Magnetic Field",3,&bfield_v);
484  EGS_TransportProperty estepem("EM ESTEPE",&the_emf->EMLMTIN);
485  EGS_TransportProperty ecut("Global Ecut",&the_bounds->ecut);
486  EGS_TransportProperty pcut("Global Pcut",&the_bounds->pcut);
487  EGS_TransportProperty smax("Global Smax",&the_etcontrol->smaxir);
488  EGS_TransportProperty estep("ESTEPE",&the_etcontrol->estepe);
489  EGS_TransportProperty ximax("Ximax",&the_etcontrol->ximax);
490  EGS_TransportProperty skind("Skin depth for BCA",
492  EGS_TransportProperty pxsec("Photon cross sections",16,the_media->pxsec);
493  EGS_TransportProperty pxsec_out("Photon cross-sections output",&the_egsio->xsec_out);
494  pxsec_out.addOption("Off");
495  pxsec_out.addOption("On");
496  EGS_TransportProperty cxsec("Compton cross sections",16,the_media->compxsec);
497  EGS_TransportProperty bc("Bound Compton scattering",&the_xoptions->ibcmp);
498  bc.addOption("Off");
499  bc.addOption("On");
500  bc.addOption("Simple");
501  bc.addOption("norej");
502  EGS_TransportProperty radc("Radiative Compton corrections",
504  radc.addOption("Off");
505  radc.addOption("On");
506  /* Photonuclear input parameters */
507  EGS_TransportProperty photonucxsec("Photonuclear cross sections",16,the_media->photonucxsec);
508  EGS_TransportProperty photonuc("Photonuclear attenuation",&the_xoptions->iphotonuc);
509  photonuc.addOption("Off");
510  photonuc.addOption("On");
511 
512  EGS_TransportProperty rayl("Rayleigh scattering",&the_xoptions->iraylr);
513  rayl.addOption("Off");
514  rayl.addOption("On");
515  rayl.addOption("custom");
516  EGS_TransportProperty ff_med("ff media names",24,MXMED,&ff_media);
517  EGS_TransportProperty ff_files("ff file names",128,MXMED, &ff_names);
518 
519  EGS_TransportProperty relax("Atomic relaxations",&the_xoptions->iedgfl);
520  relax.addOption("Off");
521  relax.addOption("On");
522  relax.addOption("eadl");
523  relax.addOption("simple");
524 
525  EGS_TransportProperty iphter("Photoelectron angular sampling",
526  &the_xoptions->iphter);
527  iphter.addOption("Off");
528  iphter.addOption("On");
529  EGS_TransportProperty spin("Spin effects",&the_xoptions->spin_effects);
530  spin.addOption("Off");
531  spin.addOption("On");
532  EGS_TransportProperty trip("Triplet production",&the_xoptions->itriplet);
533  trip.addOption("Off");
534  trip.addOption("On");
535  EGS_TransportProperty eii("Electron Impact Ionization",16,the_media->eiixsec);
536  EGS_TransportProperty brem("Brems cross sections",&the_xoptions->ibr_nist);
537  brem.addOption("BH");
538  brem.addOption("NIST");
539  brem.addOption("NRC");
540  EGS_TransportProperty bang("Brems angular sampling",&the_xoptions->ibrdst);
541  bang.addOption("Simple");
542  bang.addOption("KM");
543  EGS_TransportProperty pair("Pair cross sections",&the_xoptions->pair_nrc);
544  pair.addOption("BH");
545  pair.addOption("NRC");
546  EGS_TransportProperty pang("Pair angular sampling",&the_xoptions->iprdst);
547  pang.addOption("Off");
548  pang.addOption("Simple");
549  pang.addOption("KM");
550  pang.addOption("Uniform");
551  pang.addOption("Blend");
552  EGS_TransportProperty tran("Electron-step algorithm",
554  tran.addOption("EGSnrc");
555  tran.addOption("PRESTA-I");
556  tran.addOption("PRESTA-II");
557  tran.addOption("default");
558  EGS_TransportProperty bca("Boundary crossing algorithm",
560  bca.addOption("Exact");
561  bca.addOption("PRESTA-I");
562 
563  if (transportp) {
564  efield.getInput(transportp);
565  bfield.getInput(transportp);
566  estepem.getInput(transportp);
567  ecut.getInput(transportp);
568  pcut.getInput(transportp);
569  smax.getInput(transportp);
571  estep.getInput(transportp);
572  ximax.getInput(transportp);
573  skind.getInput(transportp);
574  pxsec.getInput(transportp);
575  pxsec_out.getInput(transportp);
576  cxsec.getInput(transportp);
577  photonucxsec.getInput(transportp);
578  photonuc.getInput(transportp);
579  bc.getInput(transportp);
580  radc.getInput(transportp);
581  rayl.getInput(transportp);
582  if (the_xoptions->iraylr>1) {
583  ff_med.getInput(transportp);
584  ff_files.getInput(transportp);
585  setRayleighData(ff_media,ff_names);
586  }
587  relax.getInput(transportp);
588  if (the_xoptions->iedgfl == 0 || the_xoptions->iedgfl == 3) {
589  the_xoptions->eadl_relax = 0;
590  }
591  else {
592  the_xoptions->eadl_relax = 1;
593  // Now, 'On' means EADL relaxation
594  if (the_xoptions->iedgfl == 1) {
595  relax.setOption(1,"eadl");
596  }
597  }
598  iphter.getInput(transportp);
599  spin.getInput(transportp);
600  trip.getInput(transportp);
601  eii.getInput(transportp);
602  setEIIData(eii.size());
603  brem.getInput(transportp);
604  bang.getInput(transportp);
605  pair.getInput(transportp);
606  pang.getInput(transportp);
607  tran.getInput(transportp);
610  }
611  bca.getInput(transportp);
612 
613  if (egsEquivStr(string("mcdf-xcom "),
614  string(the_media->pxsec).substr(0,pxsec.size()))) {
615  the_xoptions->mcdf_pe_xsections = 1;
616  string("xcom ").copy(the_media->pxsec,16,0);
617  }
618  else if (egsEquivStr(string("mcdf-epdl "),
619  string(the_media->pxsec).substr(0,pxsec.size()))) {
620  the_xoptions->mcdf_pe_xsections = 1;
621  string("epdl ").copy(the_media->pxsec,16,0);
622  }
623  else {
624  the_xoptions->mcdf_pe_xsections = 0;
625  }
626 
627  // iedgfl == 3 implies eadl_relax == 0
628  if (the_xoptions->iedgfl == 3 &&
629  the_xoptions->mcdf_pe_xsections == 1) {
630  egsWarning("\n**** Warning:"
631  "\n Simplified atomic relaxation not allowed"
632  "\n with shellwise PE cross sections. Resetting"
633  "\n to detailed EADL atomic relaxation!!!\n\n");
634  the_xoptions->eadl_relax = 1;
635  the_xoptions->iedgfl = 2;
636  relax.setOption(1,"eadl");
637  }
638  delete transportp;
639  }
640 
641  if (do_hatch) {
642 
643  egsHatch();
644  F77_OBJ_(set_elastic_parameter,SET_ELASTIC_PARAMETER)();
645 
648 
649  if (nmed_old > 0) {
650  delete [] i_ededx;
651  delete [] i_pdedx;
652  delete [] i_esig;
653  delete [] i_psig;
654  delete [] i_ebr1;
655  delete [] i_pbr1;
656  delete [] i_pbr2;
657  delete [] i_gmfp;
658  delete [] i_gbr1;
659  delete [] i_gbr2;
660  delete [] i_cohe;
661  delete [] i_photonuc;
662  }
663  if (nmed > 0) {
664  i_ededx = new EGS_Interpolator [nmed];
665  i_pdedx = new EGS_Interpolator [nmed];
666  i_esig = new EGS_Interpolator [nmed];
667  i_psig = new EGS_Interpolator [nmed];
668  i_ebr1 = new EGS_Interpolator [nmed];
669  i_pbr1 = new EGS_Interpolator [nmed];
670  i_pbr2 = new EGS_Interpolator [nmed];
671  i_gmfp = new EGS_Interpolator [nmed];
672  i_gbr1 = new EGS_Interpolator [nmed];
673  i_gbr2 = new EGS_Interpolator [nmed];
674  i_cohe = new EGS_Interpolator [nmed];
676  for (int imed=0; imed<nmed; imed++) {
677  int type;
678  int imed1 = imed+1;
679  type = 3;
680  __init_interpolator(1,imed1,type,i_ededx[imed]);
681  type = 4;
682  __init_interpolator(1,imed1,type,i_pdedx[imed]);
683  type = 1;
684  __init_interpolator(1,imed1,type,i_esig[imed]);
685  type = 2;
686  __init_interpolator(1,imed1,type,i_psig[imed]);
687  type = 5;
688  __init_interpolator(1,imed1,type,i_ebr1[imed]);
689  type = 6;
690  __init_interpolator(1,imed1,type,i_pbr1[imed]);
691  type = 7;
692  __init_interpolator(1,imed1,type,i_pbr2[imed]);
693  type = 1;
694  __init_interpolator(0,imed1,type,i_gmfp[imed]);
695  type = 2;
696  __init_interpolator(0,imed1,type,i_gbr1[imed]);
697  type = 3;
698  __init_interpolator(0,imed1,type,i_gbr2[imed]);
699  type = 4;
700  __init_interpolator(0,imed1,type,i_cohe[imed]);
701  type = 5;
702  __init_interpolator(0,imed1,type,i_photonuc[imed]);
703  }
704  }
705  }
706 
707  egsInformation("\n\nThe following media are defined:\n"
708  "================================\n\n");
709  EGS_Float Emax = source ? source->getEmax() : 0;
710  bool data_ok = true;
711  for (j=0; j<nmed; j++) {
712  int imed = ind[j]-1;
713  egsInformation("%3d %-24s AE=%7.4f AP=%7.4f %d\n",j,
714  geometry->getMediumName(j),the_thresh->ae[imed],
715  the_thresh->ap[imed],imed);
716  if (Emax > the_thresh->ue[imed]-the_useful->rm ||
717  Emax > the_thresh->up[imed]) {
718  egsInformation(" upper limits (UE=%g UP=%g) not enough for "
719  "maximum source energy (%g)\n",the_thresh->ue[imed],
720  the_thresh->up[imed],Emax);
721  data_ok = false;
722  }
723  }
724  if (!data_ok) {
725  delete [] ind;
726  return 3;
727  }
728 
729  egsInformation("\n\nTransport parameter and cross section options:\n"
730  "==============================================\n");
731  int nc = 50;
732  // Recover initial input string to reflect what's actually being used
733  // the_media->pxsec was already used in egsHatch() to read photon data
734  if (the_xoptions->mcdf_pe_xsections &&
735  egsEquivStr(string("xcom "),
736  string(the_media->pxsec).substr(0,pxsec.size()))) {
737  string("mcdf-xcom ").copy(the_media->pxsec,16,0);
738  }
739  else if (the_xoptions->mcdf_pe_xsections &&
740  egsEquivStr(string("epdl "),
741  string(the_media->pxsec).substr(0,pxsec.size()))) {
742  string("mcdf-epdl ").copy(the_media->pxsec,16,0);
743  }
744 
745  if (!isspace(the_media->pxsec[0])) {
746  pxsec.info(nc);
747  }
748  if (!isspace(the_media->compxsec[0])) {
749  cxsec.info(nc);
750  }
751  pcut.info(nc);
752  pair.info(nc);
753  pang.info(nc);
754  trip.info(nc);
755  bc.info(nc);
756  radc.info(nc);
757  rayl.info(nc);
758  relax.info(nc);
759  iphter.info(nc);
760  photonuc.info(nc);
761  photonucxsec.info(nc);
762  egsInformation("\n");
763  ecut.info(nc);
764  brem.info(nc);
765  bang.info(nc);
766  spin.info(nc);
767  eii.info(nc);
768  smax.info(nc);
769  estep.info(nc);
770  ximax.info(nc);
771  bca.info(nc);
772  skind.info(nc);
773  tran.info(nc);
774  if (efield.size()==3) {
775  efield.info(nc);
776  the_emf->ExIN=efield_v[0];
777  the_emf->EyIN=efield_v[1];
778  the_emf->EzIN=efield_v[2];
779  if (efield_v[0]*efield_v[0]+efield_v[1]*efield_v[1]+efield_v[2]*efield_v[2] > 0) {
780  the_emf->emfield_on = true;
781  }
782 
783  }
784  if (bfield.size()==3) {
785  bfield.info(nc);
786  the_emf->BxIN=bfield_v[0];
787  the_emf->ByIN=bfield_v[1];
788  the_emf->BzIN=bfield_v[2];
789  the_emf->Bx=bfield_v[0];
790  the_emf->By=bfield_v[1];
791  the_emf->Bz=bfield_v[2];
792  }
793  if (efield.size()==3 || bfield.size()==3) {
794  estepem.info(nc);
795  if (bfield_v[0]*bfield_v[0]+bfield_v[1]*bfield_v[1]+bfield_v[2]*bfield_v[2] > 0) {
796  the_emf->emfield_on = true;
797  }
798  }
799  egsInformation("==============================================\n\n");
800 
801  delete [] ind;
802  return 0;
803 }
804 
805 /**********************************************************************
806  Set media and corresponding ff file names for custom Rayleigh data
807 
808  Since the_xoptions->iraylr is reset to 1 here, output will only
809  indicate coherent scatter is on, but not whether custom input
810  requested. Perhaps should make EGS check for iraylm(medium)>0
811  instead of iraylm(medium)=1. The latter has been kept over
812  the years for historic reasons.
813  EMH April 6th 2011
814 **********************************************************************/
816  const vector<string> &str_medium,
817  const vector<string> &str_file) {
818  the_xoptions->iraylr=1;
819  int i,k;
820  for (i=0; i<str_medium.size(); i++) {
821  for (k=0; k<str_medium[i].size(); k++) {
822  the_rayleigh->ff_media[i][k]= str_medium[i][k];
823  }
824  }
825  for (i=0; i<str_file.size(); i++) {
826  for (k=0; k<str_file[i].size(); k++) {
827  the_rayleigh->ff_file[i][k]=str_file[i][k];
828  }
829  }
830 }
831 
832 /**********************************************************************
833  Set EII flag and xsection file name.
834 
835  The EII input is of a mixed type, i.e., one should be able to turn
836  this option on/off, but there is also the possibility of using an
837  arbitrary EII xsection compilation, including the available EII
838  xsections with EGSnrc.
839  EMH April 6th 2011
840 **********************************************************************/
842  string str_eii,off,on;
843  the_xoptions->eii_flag = 1;
844  string ik="ik";
845  for (int i=0; i<len; i++) {
846  str_eii += toupper(the_media->eiixsec[i]);
847  off+=' ';
848  on+=' ';
849  }
850  off[0]='O';
851  off[1]='F';
852  off[2]='F';
853  on[0]='O';
854  on[1]='N';
855  if (str_eii == off) {
856  the_xoptions->eii_flag = 0;
857  }
858  else if (str_eii == on) {
859  ik.copy(the_media->eiixsec,16,0);
860  }
861 }
862 
863 #ifdef GDEBUG
864  #define MAX_STEP 100000
865  EGS_Vector steps_x[MAX_STEP], steps_u[MAX_STEP];
866  int steps_ireg[MAX_STEP], steps_inew[MAX_STEP];
867  EGS_Float steps_ustepi[MAX_STEP], steps_ustepf[MAX_STEP];
868  int steps_n = 0;
869 #endif
870 
872  if (!p.wt) {
873  return 0;
874  }
875  the_stack->E[0] = (p.q) ? p.E + the_useful->rm : p.E;
876  the_stack->x[0] = p.x.x;
877  the_stack->y[0] = p.x.y;
878  the_stack->z[0] = p.x.z;
879  the_stack->u[0] = p.u.x;
880  the_stack->v[0] = p.u.y;
881  the_stack->w[0] = p.u.z;
882  the_stack->dnear[0] = 0;
883  the_stack->wt[0] = p.wt;
884  the_stack->ir[0] = p.ir + 2;
885  the_stack->iq[0] = p.q;
886  the_stack->latch[0] = p.latch;
887  the_stack->np = 1;
888 #ifdef GDEBUG
889  steps_n = 0;
890 #endif
891  egsShower();
892  return 0;
893 }
894 
896  egsFinish();
897  //egsSetDefaultIOFunctions();
898  io_flag = 0;
899  // so that if extra output is being produced somewhere, it does not
900  // go to fort.6 as it would without this as all fortran units are
901  // closed after egsFinish().
902 }
903 
906  egsInformation("finishSimulation(%s) %d\n",app_name.c_str(),err);
907  if (err <= 0) {
908  return err;
909  }
910  //egsInformation("\n\n********** I'm last job! **********\n\n");
911  //return 0;
912  // if err is positive, this is the last job in a parallel run
913  // => we have to combine the results
914  // EGS_Application::finishSimulation() has already called
915  // our finishRun() which has called egs_finish => all data has been
916  // moved to the user code directory and the working directory has been
917  // reset to be the user code directory. We must now reset the
918  // output_file name and re-open units.
920  the_egsio->i_parallel = 0;
921  i_parallel=0;
922  int flag = 0;
923  egsOpenUnits(&flag);
924  // The following is necessary because finishRun() was called from
925  // EGS_Application::finishSimulation() and this resets I/O
926  // to stdout and stderr. But if we are here, we are combining
927  // the results of a parallel run and we want the output to go
928  // into the combined run log file.
929  //egsSetInfoFunction(Information,__write_to_fortran_file);
930  //egsSetInfoFunction(Warning,__write_to_fortran_file);
931  //egsSetInfoFunction(Fatal,__write_to_fortran_file_and_exit);
932  io_flag = 1;
933  final_job = true;
936  err = combineResults();
937  if (err) {
938  return err;
939  }
941  for (int j=0; j<a_objects_list.size(); ++j) {
942  a_objects_list[j]->reportResults();
943  }
944  outputResults();
945  finishRun();
946  return 0;
947 }
948 
950  int err = EGS_Application::outputData();
951  if (err) {
952  return err;
953  }
954  EGS_I32 np, ip;
955  egsGetRNGPointers(&np,&ip);
956  if (np < 1) {
957  return 11;
958  }
959  if (np > 10000000) {
960  egsWarning("EGS_AdvancedApplication::outputData(): egsGetRNGPointers"
961  " returns a huge number? (%d)\n",np);
962  return 12;
963  }
964  EGS_Float *array = new EGS_Float [np];
965  egsGetRNGArray(array);
966  (*data_out) << " " << np << " " << ip << endl;
967  for (int j=0; j<np; j++) {
968  (*data_out) << array[j] << " ";
969  }
970  (*data_out) << endl;
971  double ch_steps, all_steps;
972  egsGetSteps(&ch_steps,&all_steps);
973  (*data_out) << ch_steps << " " << all_steps << endl;
974  delete [] array;
975  return data_out->good() ? 0 : 13;
976 }
977 
979  int err = EGS_Application::readData();
980  if (err) {
981  return err;
982  }
983  int np, ip;
984  (*data_in) >> np >> ip;
985  if (np < 1) {
986  return 11;
987  }
988  if (np > 10000000) {
989  egsWarning("EGS_AdvancedApplication::readData(): got huge size "
990  "for the mortran random array? (%d)\n",np);
991  return 12;
992  }
993  EGS_Float *array = new EGS_Float [np];
994  for (int j=0; j<np; j++) {
995  (*data_in) >> array[j];
996  }
997  if (!data_in->good()) {
998  return 13;
999  }
1000  egsSetRNGState(&ip,array);
1001  delete [] array;
1002  double ch_steps, all_steps;
1003  (*data_in) >> ch_steps >> all_steps;
1004  egsSetSteps(&ch_steps,&all_steps);
1005  return data_in->good() ? 0 : 13;
1006 }
1007 
1009  int err = EGS_Application::addState(data);
1010  if (err) {
1011  return err;
1012  }
1013  int np, ip;
1014  data >> np >> ip;
1015  if (np < 1) {
1016  egsWarning("Got np=%d\n",np);
1017  return 11;
1018  }
1019  if (np > 10000000) {
1020  egsWarning("EGS_AdvancedApplication::addState(): got huge size "
1021  "for the mortran random array? (%d)\n",np);
1022  return 12;
1023  }
1024  EGS_Float *array = new EGS_Float [np];
1025  for (int j=0; j<np; j++) {
1026  data >> array[j];
1027  }
1028  if (!data.good()) {
1029  return 13;
1030  }
1031  egsSetRNGState(&ip,array);
1032  delete [] array;
1033  double ch_steps, all_steps;
1034  data >> ch_steps >> all_steps;
1035  double ch_steps_old, all_steps_old;
1036  egsGetSteps(&ch_steps_old,&all_steps_old);
1037  ch_steps += ch_steps_old;
1038  all_steps += all_steps_old;
1039  egsSetSteps(&ch_steps,&all_steps);
1040  return data.good() ? 0 : 13;
1041 }
1042 
1045  double ch_steps = 0, all_steps = 0;
1046  egsSetSteps(&ch_steps,&all_steps);
1047 }
1048 
1050  EGS_I64 nused = EGS_Application::randomNumbersUsed();
1051  EGS_I32 np, ip;
1052  egsGetRNGPointers(&np,&ip);
1053  if (ip <= np) {
1054  nused -= (np+1-ip);
1055  }
1056  return nused;
1057 }
1058 
1059 void EGS_AdvancedApplication::appInformation(const char *msg) {
1060  if (!msg) {
1061  return;
1062  }
1063  if (!io_flag) {
1065  }
1066  else {
1067  int ounit = the_egsio->i_log;
1068  F77_OBJ_(egs_write_string,EGS_WRITE_STRING)(&ounit,msg,strlen(msg));
1069  checkDeviceFull(stdout);
1070  }
1071 }
1072 
1073 void EGS_AdvancedApplication::appWarning(const char *msg) {
1074  if (!msg) {
1075  return;
1076  }
1077  if (!io_flag) {
1079  }
1080  else {
1081  int ounit = the_egsio->i_log;
1082  F77_OBJ_(egs_write_string,EGS_WRITE_STRING)(&ounit,msg,strlen(msg));
1083  checkDeviceFull(stdout);
1084  }
1085 }
1086 
1087 void EGS_AdvancedApplication::appFatal(const char *msg) {
1088  if (!io_flag) {
1090  }
1091  else {
1092  if (msg) {
1093  int ounit = the_egsio->i_log;
1094  F77_OBJ_(egs_write_string,EGS_WRITE_STRING)(&ounit,msg,strlen(msg));
1095  }
1096  exit(1);
1097  }
1098 }
1099 
1101  double &all_steps) const {
1102  egsGetSteps(&ch_steps,&all_steps);
1103 }
1104 
1106  if (geometry->hasRhoScaling()) {
1107  int ireg = the_stack->ir[the_stack->np-1] - 2;
1108  EGS_Float rho = geometry->getRelativeRho(ireg);
1109  the_useful->rhor = rho;
1110  the_useful->rhor_new = rho;
1111  }
1112  else {
1113  the_useful->rhor = 1;
1114  the_useful->rhor_new = 1;
1115  }
1116  if (geometry->hasBScaling()) {
1117  int ireg = the_stack->ir[the_stack->np-1] - 2;
1118  EGS_Float bf = geometry->getBScaling(ireg);
1119  the_emf->Bx = bf*the_emf->BxIN;
1120  the_emf->By = bf*the_emf->ByIN;
1121  the_emf->Bz = bf*the_emf->BzIN;
1122  the_emf->Bx_new = the_emf->Bx;
1123  the_emf->By_new = the_emf->By;
1124  the_emf->Bz_new = the_emf->Bz;
1125  }
1126  else {
1127  the_emf->Bx = the_emf->BxIN;
1128  the_emf->By = the_emf->ByIN;
1129  the_emf->Bz = the_emf->BzIN;
1130  the_emf->Bx_new = the_emf->Bx;
1131  the_emf->By_new = the_emf->By;
1132  the_emf->Bz_new = the_emf->Bz;
1133  }
1134 }
1135 
1136 void EGS_AdvancedApplication::enterNewRegion() {
1137  if (geometry->hasRhoScaling()) {
1138  int ireg = the_epcont->irnew-2;
1139  if (ireg >= 0) {
1140  EGS_Float rho = geometry->getRelativeRho(ireg);
1141  the_useful->rhor_new = rho;
1142  }
1143  }
1144  else {
1145  the_useful->rhor_new = 1;
1146  }
1147  if (geometry->hasBScaling()) {
1148  int ireg = the_epcont->irnew-2;
1149  if (ireg >= 0) {
1150  EGS_Float bf = geometry->getBScaling(ireg);
1151  the_emf->Bx_new = bf*the_emf->BxIN;
1152  the_emf->By_new = bf*the_emf->ByIN;
1153  the_emf->Bz_new = bf*the_emf->BzIN;
1154  }
1155  }
1156  else {
1157  the_emf->Bx_new = the_emf->BxIN;
1158  the_emf->By_new = the_emf->ByIN;
1159  the_emf->Bz_new = the_emf->BzIN;
1160  }
1161 }
1162 
1164  rndm->saveState();
1165  int np,ip;
1166  egsGetRNGPointers(&np,&ip);
1167  if (n_rng_buffer > 0) {
1168  if (np != n_rng_buffer) egsFatal("EGS_AdvancedApplication::saveRNGState:\n"
1169  " new buffer size (%d) is not the same as previous size(%d) ?\n",
1170  np,n_rng_buffer);
1171  }
1172  else {
1173  rng_buffer = new EGS_Float [np];
1174  n_rng_buffer = np;
1175  }
1176  i_rng_buffer = ip;
1177  egsGetRNGArray(rng_buffer);
1178 }
1179 
1181  rndm->resetState();
1182  if (n_rng_buffer > 0) {
1183  egsSetRNGState(&i_rng_buffer,rng_buffer);
1184  }
1185 }
1186 
1187 //************************************************************
1188 // Utility functions for use with ausgab dose scoring object
1189 //************************************************************
1190 // Returns density for medium ind
1191 EGS_Float EGS_AdvancedApplication::getMediumRho(int ind) {
1192  // handle the negative medium index for vacuum
1193  if (ind < 0) {
1194  return 0.0;
1195  }
1196  return the_media->rho[ind];
1197 }
1198 // Returns edep
1199 EGS_Float EGS_AdvancedApplication::getEdep() {
1200  return the_epcont->edep;
1201 }
1202 // Sets edep
1203 void EGS_AdvancedApplication::setEdep(EGS_Float edep) {
1204  the_epcont->edep = edep;
1205 }
1206 // Get the ecut
1207 EGS_Float EGS_AdvancedApplication::getEcut() {
1208  return the_bounds->ecut;
1209 }
1210 // Get the pcut
1211 EGS_Float EGS_AdvancedApplication::getPcut() {
1212  return the_bounds->pcut;
1213 }
1214 // Returns rest mass
1215 EGS_Float EGS_AdvancedApplication::getRM() {
1216  return the_useful->rm;
1217 }
1218 
1219 //************************************************************
1220 // Utility functions for fluence scoring objects
1221 //************************************************************
1222 
1223 // Turns ON/OFF EGSnrc internal radiative splitting (UBS)
1224 void EGS_AdvancedApplication::setRadiativeSplitting(const EGS_Float &nsplit) {
1225  the_egsvr->nbr_split = nsplit;
1226 }
1227 
1228 // Turns ON/OFF EGSnrc internal Russian Roulette + UBS
1229 void EGS_AdvancedApplication::setRussianRoulette(const EGS_Float &iSwitchRR) {
1230  if (iSwitchRR > 1.0) {
1231  the_egsvr->i_play_RR = 1;
1232  the_egsvr->prob_RR = 1.0/iSwitchRR;
1233  the_egsvr->nbr_split = iSwitchRR;
1234  }
1235  else {
1236  the_egsvr->i_play_RR = 0;
1237  the_egsvr->nbr_split = 1;
1238  }
1239 }
1240 
1241 // Splits top particle into nsplit particles uniformly in 4Pi
1242 void EGS_AdvancedApplication::splitTopParticleIsotropically(const EGS_Float &fsplit) {
1243  // Reset particle pointer
1245  /* Initialize local variables */
1246  int np = the_stack->np-1,
1247  the_latch = the_stack->latch[np],
1248  ir = the_stack->ir[np],
1249  iq = the_stack->iq[np];
1250  the_stack->wt[np] /= fsplit;
1251  double E = the_stack->E[np];
1252  EGS_Float x = the_stack->x[np], y = the_stack->y[np], z = the_stack->z[np],
1253  wthin = the_stack->wt[np], dnear = the_stack->dnear[np];
1254  EGS_Float u,v,w;
1255  /* If fsplit is a non-integer, sample between int(fsplit) and int(split)+1 */
1256  int nsplit = int(fsplit);
1257  EGS_Float dsplit = fsplit - nsplit;
1258  if (dsplit > 0) { // non-integer splitting number
1259  if (rndm->getUniform() < dsplit) {
1260  ++nsplit;
1261  }
1262  }
1263  for (int i=0; i < nsplit; i++) {
1264  np++;
1265  if (np >= MXSTACK) {
1266  egsFatal("\n\n******************************************\n"
1267  "ERROR: In EGS_AdvancedApplication::splitTopParticleIsotropically() :\n"
1268  "max. stack depth MXSTACK=%d < np=%d\n"
1269  "Stack overflow due to splitting!\n"
1270  "******************************************\n"
1271  ,MXSTACK,np);
1272  }
1273  the_stack->x[np] = x;
1274  the_stack->y[np] = y;
1275  the_stack->z[np] = z;
1276  the_stack->iq[np]= iq;
1277  the_stack->dnear[np] = dnear;
1278  the_stack->latch[np] = the_latch;
1279  the_stack->ir[np]= ir;
1280  the_stack->E[np] = E;
1281  the_stack->wt[np]=wthin;
1282  // Particles isotropically distributed in space
1283  w = 2*rndm->getUniform()-1;
1284  EGS_Float sinz = 1-w*w;
1285  if (sinz > epsilon) {
1286  sinz = sqrt(sinz);
1287  EGS_Float cphi, sphi;
1288  EGS_Float phi = 2*M_PI*rndm->getUniform();
1289  cphi = cos(phi);
1290  sphi = sin(phi);
1291  u = sinz*cphi;
1292  v = sinz*sphi;
1293  }
1294  else {
1295  u = 0;
1296  v = 0;
1297  }
1298 
1299  the_stack->u[np] = u;
1300  the_stack->v[np] = v;
1301  the_stack->w[np] = w;
1302 
1303  }
1304 
1305  the_stack->np = np+1;
1306 }
1307 
1308 //************************************************************
1309 // Utility functions for fluence scoring objects
1310 //************************************************************
1311 EGS_Float EGS_AdvancedApplication::getTVSTEP() {
1312  return the_epcont->tvstep;
1313 }
1314 
1315 EGS_Interpolator *EGS_AdvancedApplication::getDEDX(const int &imed, const int &iq) {
1316  if (iq == -1) {
1317  return &i_ededx[imed];
1318  }
1319  else if (iq == 1) {
1320  return &i_pdedx[imed];
1321  }
1322  else {
1323  return 0;
1324  }
1325 }
1326 
1327 void EGS_AdvancedApplication::setLatch(const int &ip, const int &latch) {
1328  the_stack->latch[ip] = latch;
1329 }
1330 
1331 void EGS_AdvancedApplication::incLatch(const int &ip, const int &increment) {
1332  the_stack->latch[ip] += increment;
1333 }
1334 
1335 int EGS_AdvancedApplication::getNp() {
1336  return the_stack->np-1;;
1337 }
1338 
1339 int EGS_AdvancedApplication::getNpOld() {
1340  return the_stack->npold-1;
1341 }
1342 
1343 //************************************************************
1344 // Utility function for ausgab phase space scoring objects
1345 //************************************************************
1346 void EGS_AdvancedApplication::setLatch(int latch) {
1347  int np = the_stack->np-1;
1348  the_stack->latch[np] = latch;
1349 }
1350 
1351 extern __extc__ void egsHowfar() {
1352  CHECK_GET_APPLICATION(app,"egsHowfar()");
1353  int np = the_stack->np-1;
1354  int ireg = the_stack->ir[np]-2;
1355  if (ireg < 0 || the_stack->wt[np] == 0) {
1356  the_epcont->idisc = 1;
1357  return;
1358  }
1359  int newmed;
1360  int inew = app->howfar(ireg,
1361  EGS_Vector(the_stack->x[np],the_stack->y[np],the_stack->z[np]),
1362  EGS_Vector(the_stack->u[np],the_stack->v[np],the_stack->w[np]),
1363  the_epcont->ustep,&newmed);
1364 #ifdef GDEBUG
1365  EGS_Float tsave = the_epcont->ustep;
1366  if (steps_n < MAX_STEP) {
1367  steps_x[steps_n] =
1368  EGS_Vector(the_stack->x[np],the_stack->y[np],the_stack->z[np]);
1369  steps_u[steps_n] =
1370  EGS_Vector(the_stack->u[np],the_stack->v[np],the_stack->w[np]);
1371  steps_ireg[steps_n] = ireg;
1372  steps_inew[steps_n] = inew;
1373  steps_ustepi[steps_n] = tsave;
1374  steps_ustepf[steps_n++] = the_epcont->ustep;
1375  }
1376  if (the_epcont->ustep < -1e-4) {
1377  egsWarning("Negative step: %g\n",the_epcont->ustep);
1378  for (int j=0; j<steps_n; j++)
1379  egsWarning("%d x=(%g,%g,%g) u=(%g,%g,%g) ireg=%d inew=%d ustepi=%g"
1380  " ustepf=%g\n",j,steps_x[j].x,steps_x[j].y,steps_x[j].z,
1381  steps_u[j].x,steps_u[j].y,steps_u[j].z,
1382  steps_ireg[j],steps_inew[j],steps_ustepi[j],
1383  steps_ustepf[j]);
1384  exit(1);
1385  }
1386 #endif
1387  if (inew != ireg) {
1388  the_epcont->irnew = inew+2;
1389  the_useful->medium_new = newmed+1;
1390  app->enterNewRegion();
1391  if (inew < 0) {
1392  the_epcont->idisc = -1; // i.e. discard after the step.
1393  the_useful->medium_new = 0;
1394  }
1395  }
1396 }
1397 
1398 extern __extc__ void egsHownear(EGS_Float *tperp) {
1399  CHECK_GET_APPLICATION(app,"egsHownear()");
1400  int np = the_stack->np-1;
1401  int ireg = the_stack->ir[np]-2;
1402  *tperp = app->hownear(ireg,EGS_Vector(the_stack->x[np],the_stack->y[np],
1403  the_stack->z[np]));
1404 }
1405 
1406 extern __extc__ void egsFillRandomArray(const EGS_I32 *n, EGS_Float *rarray) {
1407  CHECK_GET_APPLICATION(app,"egsFillRandomArray()");
1408  app->fillRandomArray(*n,rarray);
1409 }
1410 
1411 //extern __extc__ int egsAusgab(const EGS_I32 *iarg) {
1412 // CHECK_GET_APPLICATION(app,"egsAusgab()");
1413 // return app->ausgab(*iarg);
1414 //}
1415 extern __extc__ void egsAusgab(EGS_I32 *iarg) {
1416  CHECK_GET_APPLICATION(app,"egsAusgab()");
1417  //*iarg = app->ausgab(*iarg);
1418  int np = the_stack->np-1;
1419  app->top_p.E = the_stack->E[np];
1420  app->top_p.wt = the_stack->wt[np];
1421  app->top_p.x = EGS_Vector(the_stack->x[np],the_stack->y[np],the_stack->z[np]);
1422  app->top_p.u = EGS_Vector(the_stack->u[np],the_stack->v[np],the_stack->w[np]);
1423  app->top_p.q = the_stack->iq[np];
1424  app->top_p.ir = the_stack->ir[np]-2;
1425  app->top_p.latch = the_stack->latch[np];
1426  app->Np = the_stack->np-1;
1427  *iarg = app->userScoring(*iarg);
1428 }
1429 
1430 extern __extc__ void egsStartParticle() {
1431  CHECK_GET_APPLICATION(app,"egsStartParticle()");
1432  int np = the_stack->np - 1;
1433  int ir = the_stack->ir[np]-2;
1434  if (ir < 0) {
1435  the_epcont->idisc = 1;
1436  return;
1437  }
1438  the_epcont->idisc = 0;
1439  the_useful->medium = app->getMedium(ir)+1;
1440  //egsInformation("start particle: ir=%d medium=%d\n",ir,the_useful->medium);
1441  app->startNewParticle();
1442 }
void resetCounter()
Reset the application to a 'pristine' state.
int readData()
Read intermediate results.
int initCrossSections()
Initialize the run-time cross section data.
int io_flag
determines how to write info
EGS_I64 randomNumbersUsed() const
Get the number of random numbers used.
virtual void finishRun()
Finish a simulation.
virtual ~EGS_AdvancedApplication()
Destructor.
EGS_Interpolator * i_gbr2
photon branching 2 interpolator
EGS_Interpolator * i_esig
electron cross section interpolator
EGS_Float * rng_buffer
RNG buffer.
EGS_Interpolator * i_pbr1
positron branching 1 interpolator
int i_rng_buffer
Pointer to the RNG buffer.
void setEIIData(EGS_I32 len)
Set EII flag and xsection file name.
int outputData()
Output intermediate results.
static string base_revision
Holds the CVS revision number of the egs_advanced_application.cpp file.
EGS_Interpolator * i_psig
positron cross section interpolator
EGS_Interpolator * i_ededx
electron stopping power interpolator
bool final_job
Is this the final job of a parallel run ?
EGS_Interpolator * i_gmfp
photon mean-free-path interpolator
EGS_Interpolator * i_cohe
photon Rayleigh interpolator
void startNewParticle()
Start the transport of a new particle.
EGS_Interpolator * i_ebr1
electron branching interpolator
int finishSimulation()
Finish the simulation.
void setRayleighData(const vector< string > &str_medium, const vector< string > &str_file)
Custom Rayleigh data setup.
void getElectronSteps(double &ch_steps, double &all_steps) const
Get the number of condensed history and all electron steps.
void setAusgabCall(AusgabCall call, bool on_or_off)
Turns on or off a call to the user scoring function ausgab()
EGS_Interpolator * i_photonuc
photonuclear interpolator
EGS_Interpolator * i_pdedx
positron stopping power interpolator
int initEGSnrcBackEnd()
Initialize the EGSnrc mortran back-end.
EGS_Interpolator * i_pbr2
positron branching 2 interpolator
int n_rng_buffer
Size of the RNG buffer.
int helpInit(EGS_Input *, bool do_hatch)
Helper function used in initCrossSections() and describeSimulation().
EGS_AdvancedApplication(int argc, char **argv)
Constructor.
int shower()
Simulate a single shower.
EGS_Interpolator * i_gbr1
photon branching 1 interpolator
void describeSimulation()
Describe the simulation.
int addState(istream &)
Add data from a parallel job.
Base class for advanced EGSnrc C++ applications.
virtual EGS_I64 randomNumbersUsed() const
Returns the number of random numbers used.
EGS_RandomGenerator * rndm
the random number generator
virtual void resetCounter()
Reset the application to a 'pristine' state.
EGS_SimpleContainer< EGS_AusgabObject * > a_objects_list
The ausgab objects.
string output_file
The output file name (no extension)
virtual void setAusgabCall(AusgabCall call, bool on_or_off)
Turns on or off a call to the user scoring function ausgab.
virtual void describeSimulation()
Describe the simulation.
virtual int addState(istream &data)
Add data from a parallel job.
virtual int finishSimulation()
Analyze and output the results.
virtual void appFatal(const char *)
Write a warning message and exit.
virtual int outputData()
Output intermediate results.
EGS_Input * input
the input to this simulation.
string app_name
The application name.
virtual void outputResults()
Output the simulation results.
string final_output_file
The final output file name.
EGS_BaseGeometry * geometry
the geometry of this simulation
AusgabCall
Possible calls to the user scoring function ausgab().
int i_parallel
Job index in parallel runs.
istream * data_in
data input stream
int app_index
the index of this application.
virtual void appInformation(const char *)
Write an information message.
virtual void describeUserCode() const
Describe the user code.
ostream * data_out
data output stream
virtual void appWarning(const char *)
Write a warning message.
EGS_RunControl * run
the run control object.
virtual int combineResults()
Combine results from parallel runs.
void checkDeviceFull(FILE *)
Check if a device holding a given stream is full.
EGS_BaseSource * source
the particle source
virtual int readData()
Read intermediate results.
static int nMedia()
Get the number of media registered so far by all geometries.
virtual EGS_Float getBScaling(int ireg) const
Get the B field scaling factor in region ireg.
bool hasBScaling() const
Does this geometry object have a B field scaling feature?
static void setActiveGeometryList(int list)
Set the currently active geometry list.
virtual EGS_Float getRelativeRho(int ireg) const
Get the relative mass density in region ireg.
bool hasRhoScaling() const
Does this geometry object have a mass density scaling feature?
static const char * getMediumName(int ind)
Get the name of medium with index ind.
virtual EGS_Float getEmax() const =0
Return the maximum energy of this source.
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 * getInputItem(const string &key) const
Same as the previous function but now ownership remains with the EGS_Input object.
Definition: egs_input.cpp:245
int getInput(const string &key, vector< string > &values) const
Assign values to an array of strings from an input identified by key.
Definition: egs_input.cpp:338
A class for fast run-time interpolations.
void initialize(int nbin, EGS_Float Xmin, EGS_Float Xmax, const EGS_Float *values)
Initialize the interpolator.
virtual void saveState()=0
Save the RNG state.
EGS_Float getUniform()
Returns a random number uniformly distributed between zero (inclusive) and 1 (exclusive).
Definition: egs_rndm.h:103
virtual void resetState()=0
Reset the state to a previously saved state.
virtual int finishSimulation()
Finish the simulation.
A class representing 3D vectors.
Definition: egs_vector.h:57
EGS_Float y
y-component
Definition: egs_vector.h:62
EGS_Float z
z-component
Definition: egs_vector.h:63
EGS_Float x
x-component
Definition: egs_vector.h:61
__extc__ void egsHownear(EGS_Float *tperp)
Calculate minimum perpendicular distance to any surrounding boundary and return it in tperp.
__extc__ void egsAusgab(EGS_I32 *iarg)
User scoring function.
__extc__ void egsFillRandomArray(const EGS_I32 *n, EGS_Float *rarray)
Get n random numbers stored in the array pointed to by rarray.
EGS_AdvancedApplication class header file.
EGS_AusgabObject interface class header file.
EGS_BaseSource class header file.
Global egspp functions header file.
EGS_Input class header file.
This file defines the C/C++ interface to the EGSnrc mortran back-end.
__extc__ struct EGS_VarianceReduction * the_egsvr
The address of the mortran egs_vr common block as a pointer to a C-structure of type EGS_VarianceRedu...
__extc__ struct EGS_Useful * the_useful
The address of the mortran USEFUL common block as a pointer to a C-structure of type EGS_Useful.
__extc__ struct EGS_Bounds * the_bounds
The address of the morrtan BOUNDS common block as a pointer to a C-structure of type EGS_Bounds.
__extc__ struct EGS_Epcont * the_epcont
The address of the mortran EPCONT common block as a pointer to a C-structure of type EGS_Epcont.
#define egsHatch
Initialize cross section data from a PEGS4 data file.
__extc__ struct EGS_emfInputs * the_emf
The address of the mortran EMF-INPUTS common block as a pointer to a C-structure of type EGS_emfInput...
#define egsHowfar
Calculate distance to the next boundary along the current direction of motion.
#define egsShower
Transport the particles currently on the particle stack.
__extc__ EGS_I32 egsAddMedium(const char *medname, EGS_I32 length)
Add a medium to the list of media.
__extc__ struct EGS_EtControl * the_etcontrol
The address of the mortran Et_Control common block as a pointer to a C-structure of type EGS_EtContro...
__extc__ struct EGS_XOptions * the_xoptions
The address of the mortran cross section options common block as a pointer to a C-structure of type E...
#define egsStartParticle
Start the transport of a new particle.
#define egsFinish
Finish a simulation.
__extc__ struct EGS_Rayleigh * the_rayleigh
The address of the mortran rayleigh_inputs common block as a pointer to a C-structure of type EGS_Ray...
__extc__ struct EGS_Stack * the_stack
The address of the mortran STACK common block as a pointer to a C-structure of type EGS_Stack.
__extc__ struct EGS_Media * the_media
The address of the mortran MEDIA common block as a pointer to a C-structure of type EGS_Media
__extc__ struct EGS_IO * the_egsio
The address of the mortran egs_io common block as a pointer to a C-structure of type EGS_IO.
__extc__ struct EGS_Thresh * the_thresh
The address of the mortran THRESH common block as a pointer to a C-structure of type EGS_Thresh.
EGS_Interpolator class header file.
Attempts to fix broken math header files.
EGS_RandomGenerator class header file.
EGS_RunControl and EGS_JCFControl class header file.
EGS_InfoFunction EGS_EXPORT egsInformation
Always use this function for reporting the progress of a simulation and any other type of information...
EGS_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
const EGS_Float epsilon
The epsilon constant for floating point comparisons.
Definition: egs_functions.h:62
EGS_InfoFunction EGS_EXPORT egsWarning
Always use this function for reporting warnings.
EGS_Float pcut_new
EGS_Float pcut
EGS_Float ecut
EGS_Float ecut_new
EGS_Float tvstep
EGS_Float ustep
EGS_I32 idisc
EGS_I32 irnew
EGS_I32 iausfl[MXAUS]
EGS_Float smaxir
EGS_Float estepe
EGS_Float ximax
EGS_Float smax_new
EGS_Float skindepth_for_bca
EGS_I32 bca_algorithm
EGS_I32 transport_algorithm
char pxsec[16]
char eiixsec[16]
char photonucxsec[16]
char compxsec[16]
EGS_Float rho[MXMED]
EGS_Vector x
position
EGS_Float E
particle energy in MeV
EGS_Vector u
direction
int ir
particle region index
int latch
latch variable (useful as a flag on many occasions)
EGS_Float wt
statistical weight
int q
particle charge
char ff_file[MXMED][128]
char ff_media[MXMED][24]
EGS_Float z[MXSTACK]
double E[MXSTACK]
EGS_Float x[MXSTACK]
EGS_Float v[MXSTACK]
EGS_Float w[MXSTACK]
EGS_I32 ir[MXSTACK]
EGS_I32 iq[MXSTACK]
EGS_Float dnear[MXSTACK]
EGS_Float y[MXSTACK]
EGS_Float u[MXSTACK]
EGS_Float wt[MXSTACK]
EGS_I32 npold
EGS_I32 np
EGS_I32 latch[MXSTACK]
EGS_Float ap[MXMED]
EGS_Float up[MXMED]
EGS_Float ae[MXMED]
EGS_Float ue[MXMED]
EGS_Float rhor
Mass density ratio.
EGS_I32 medium
Current medium.
EGS_Float rm
Electron rest energy.
EGS_Float rhor_new
Mass density ratio in the new region.
EGS_I32 medium_new
Medium in the new region.
EGS_I32 ibr_nist
EGS_I32 pair_nrc
EGS_I32 spin_effects
EGS_I32 iphotonuc
EGS_I32 radc_flag
EGS_I32 itriplet
EGS_I32 eii_flag
EGS_Float ExIN