EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_phsp_scoring.cpp
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ phase space scoring object
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: Blake Walters, 2018
25 #
26 # Contributors:
27 #
28 ###############################################################################
29 */
30 
31 
37 #include <sstream>
38 #include <fstream>
39 #include <string>
40 #include <cstdlib>
41 
42 #include "egs_phsp_scoring.h"
43 #include "egs_input.h"
44 #include "egs_functions.h"
45 #include "iaea_phsp.h"
46 
47 EGS_PhspScoring::EGS_PhspScoring(const string &Name,
48  EGS_ObjectFactory *f) :
49  EGS_AusgabObject(Name,f), phsp_index(0), store_max(1000), phsp_file(0),
50  count(0), countg(0), emin(1.e30), emax(-1.e30), first_flush(true), is_restart(false) {
51  otype = "EGS_PhspScoring";
52 }
53 
54 EGS_PhspScoring::~EGS_PhspScoring() {
55 }
56 
57 void EGS_PhspScoring::setApplication(EGS_Application *App) {
59  if (!app) {
60  return;
61  }
62 
63  char buf[512];//useful character buffer
64  //set up the stack of particles to output to the phase space file
65  //1000 particles at a time
66  p_stack = new Particle[store_max];
67 
68  description = "\n*******************************************\n";
69  description += "Phase Space Scoring Object (";
70  description += name;
71  description += ")\n";
72  description += "*******************************************\n";
73  if (score_type==0) {
74  description += "\n Will output phase space for particles crossing surfaces of geometry:\n";
75  description += phsp_geom->getName();
76  description += "\n Particles scored on: ";
77  if (scoredir == 0) {
78  description += "entering and exiting phase space geometry";
79  }
80  else if (scoredir == 1) {
81  description += "entering phase space geometry";
82  }
83  else if (scoredir == 2) {
84  description += "exiting phase space geometry";
85  }
86  }
87  else if (score_type==1) {
88  description += "\n Will output phase space for the following exit/entry region pairs:\n";
89  description += " Exit regions:\n";
90  for (int i=0; i<fromreg.size(); i++) {
91  sprintf(buf,"%d ",fromreg[i]);
92  description += buf;
93  }
94  description += "\n Entry regions:\n";
95  for (int i=0; i<fromreg.size(); i++) {
96  sprintf(buf,"%d ",toreg[i]);
97  description += buf;
98  }
99  description += "\n";
100  //now, for every region, set up a (possibly empty) array of entry regions
101  int nreg = app->getnRegions();
102  for (int ir=0; ir<nreg; ir++) {
103  vector<int> tolist;
104  for (int i=0; i<fromreg.size(); i++) {
105  if (fromreg[i]==ir) {
106  tolist.push_back(toreg[i]);
107  }
108  }
109  from_to.push_back(tolist);
110  }
111  }
112 
113  //construct name of phase space file -- opening to occur later
114  if (phspoutdir == "") {
115  phspoutdir = app->getAppDir();
116  }
117  if (oformat==0) {
118  description += "\n Data will be output in EGSnrc format.\n";
119  if (app->getNparallel()>0) {
120  sprintf(buf,"%s_w%d.egsphsp1",getObjectName().c_str(),app->getIparallel());
121  }
122  else {
123  sprintf(buf,"%s.egsphsp1",getObjectName().c_str());
124  }
125  phsp_fname=egsJoinPath(phspoutdir,buf);
126  description += "\n Phase space file name:\n";
127  description += phsp_fname;
128  }
129  else if (oformat==1) {
130  description += "\n Data will be output in IAEA format.\n";
131  if (app->getNparallel()>0) {
132  sprintf(buf,"%s_w%d.1",getObjectName().c_str(),app->getIparallel());
133  }
134  else {
135  sprintf(buf,"%s.1",getObjectName().c_str());
136  }
137  phsp_fname=egsJoinPath(phspoutdir,buf);
138 
139  //need to add a null terminator
140  len = phsp_fname.length();
141  phsp_fname_char = new char[len+1];
142  phsp_fname.copy(phsp_fname_char,len,0);
143  phsp_fname_char[len]='\0'; //null terminator on string
144 
145  //iaea particle charge
146  iaea_q_type[0]=2; //e-
147  iaea_q_type[1]=1; //photon
148  iaea_q_type[2]=3; //e+
149 
150  iaea_id = 1;
151 
152  //set up extra float and extra long array indices
153  latch_ind = 2;//type of latch as defined by iaea
154  iaea_n_extra_long=1; //only store latch
155  iaea_i_latch=0; // position of latch in array
156  if (score_mu) {
157  iaea_n_extra_float=1;
158  mu_ind = 0; //set type to generic float
159  iaea_i_mu=0; //position of mu in array
160  }
161  else {
162  iaea_n_extra_float=0; //no extra floats
163  }
164 
165  description += "\n Phase space file names:\n";
166  description += " Header file: " + phsp_fname + ".IAEAheader\n";
167  description += " Data file: " + phsp_fname + ".IAEAphsp";
168  string xyzname[3] = {"X", "Y", "Z"};
169  for (int i=0; i<3; i++) {
170  if (xyz_is_const[i]) {
171  ostringstream xyz;
172  xyz << xyzscore[i];
173  description += "\n Data scored at constant " + xyzname[i] + " = " + xyz.str() + " cm";
174  }
175  }
176  }
177  description += "\n Particles scored: ";
178  if (ocharge == 0) {
179  description += "all";
180  }
181  else if (ocharge == 1) {
182  description += "photons";
183  }
184  else if (ocharge == 2) {
185  description += "charged";
186  }
187  if (oformat ==1 && score_mu) {
188  description += "\n mu will be scored (if available)";
189  }
190  if (oformat ==0 && score_mc) {
191  description += "\n will score multiple crossers (and descendents)";
192  }
193 }
194 
195 //final buffer flush and then close file
196 //we want to output some data about particles scored
197 //may ultimately want to allow the user to define scoring zones for output
198 void EGS_PhspScoring::reportResults() {
199  flushBuffer();
200  if (oformat == 1) { //iaea
201  int iaea_iostat;
202  iaea_destroy_source(&iaea_id,&iaea_iostat);
203  if (iaea_iostat<0) {
204  egsFatal("\n EGS_PhspScoring: Error closing phase space file.\n");
205  }
206  }
207  else if (oformat == 0) {
208  phsp_file.close();
209  }
210  egsInformation("\n======================================================\n");
211  egsInformation("Phase Space Scoring Object(%s)\n",name.c_str());
212  egsInformation("======================================================\n");
213  if (oformat==1) {
214  egsInformation("\n IAEA format phase space output:\n");
215  egsInformation(" Header file: %s.IAEAheader\n",phsp_fname.c_str());
216  egsInformation(" Data file: %s.IAEAphsp\n",phsp_fname.c_str());
217  }
218  else if (oformat == 0) {
219  egsInformation("\n EGSnrc format phase space output:\n");
220  egsInformation(" Data file: %s\n",phsp_fname.c_str());
221  }
222  float emintmp;
223  if (count == countg) {
224  emintmp = 0.0;
225  }
226  else {
227  emintmp = emin;
228  }
229  egsInformation("Summary of scored data:\n");
230  egsInformation("=> total no. of particles = %lld \n", count);
231  egsInformation("=> no. of photons = %lld \n", countg);
232  egsInformation("=> max. k.e. of all particles = %g MeV\n",emax);
233  egsInformation("=> min. k.e. of charged particles = %g MeV\n",emintmp);
234  egsInformation("=> no. of primary histories represented = %lld\n",last_case);
235  egsInformation("\n======================================================\n");
236 }
237 
238 //store 1000 particles at a time in p_stack
239 //if we're at 1000, actually write the particles to the file and update
240 //the header info
241 //also, keep track of phase space file counters, min., max. energy
242 void EGS_PhspScoring::storeParticle(EGS_I64 ncase) {
243 
244  //if user requested mu scoring, check if mu is available
245  if (score_mu && app->getMU() < 0) {
246  egsWarning("\nEGS_PhspScoring: User requested mu scoring, but mu is inavailable with this source.\n");
247  egsWarning("Turning off mu scoring.\n");
248  score_mu=false;
249  }
250 
251  //counters, min. and max. k.e.
252  EGS_Float prm = app->getRM();
253  count++;
254  if (app->top_p.q==0) {
255  countg++;
256  }
257  if (app->top_p.E-abs(app->top_p.q)*prm > emax) {
258  emax = app->top_p.E-abs(app->top_p.q)*prm;
259  }
260  if (app->top_p.q != 0 && app->top_p.E - prm < emin) {
261  emin = app->top_p.E - prm;
262  }
263 
264  //store particle data in p_stack
265  //set -ve energy marker if this is a new primary hist.
266  double E = app->top_p.E;
267  if (ncase != last_case) {
268  E = -E;
269  last_case = ncase;
270  }
271  p_stack[phsp_index].E = E;
272  p_stack[phsp_index].wt = app->top_p.wt;
273  p_stack[phsp_index].x = app->top_p.x.x;
274  p_stack[phsp_index].y = app->top_p.x.y;
275  p_stack[phsp_index].z = app->top_p.x.z;
276  p_stack[phsp_index].u = app->top_p.u.x;
277  p_stack[phsp_index].v = app->top_p.u.y;
278  p_stack[phsp_index].w = app->top_p.u.z;
279  p_stack[phsp_index].q = app->top_p.q;
280  if (score_mu) {
281  p_stack[phsp_index].mu = app->getMU();
282  }
283  p_stack[phsp_index++].latch = app->top_p.latch;
284 
285  if (phsp_index > store_max - 1) {
286  //write store_max particles to the file and reset phsp_index counter
287  flushBuffer();
288  }
289 }
290 
291 //open the phase space file for writing/appending data
292 void EGS_PhspScoring::openPhspFile() const {
293 //the file has already been named at this point
294  if (oformat==0) { //EGSnrc format
295  if (is_restart) {
296  phsp_file.open(phsp_fname.c_str(),ios::binary|ios::out|ios::in);
297  if (!(phsp_file))
298  egsFatal("\nEGS_PhspScoring: Failed to open phase space file %s for appending.\n",
299  phsp_fname.c_str());
300  //check that total no. of particles in header = total no. read from .egsdat file
301  unsigned int count4;
302  phsp_file.seekg(5,ios::beg);
303  phsp_file.read((char *) &count4, sizeof(unsigned int));
304  if (count4 != countprev) {
305  egsFatal("\nEGS_PhspScoring: Particle no. mismatch between %s and .egsdat file.\n",phsp_fname.c_str());
306  }
307  //go to the end of the file
308  phsp_file.seekp(0,ios::end);
309  if (phsp_file.tellp() != 28 + countprev*sizeof(egs_phsp_write_struct)) {
310  egsFatal("\nEGS_PhspScoring: File size mismatch in %s and .egsdat file.\n",phsp_fname.c_str());
311  }
312  }
313  else {
314  phsp_file.open(phsp_fname.c_str(),ios::binary|ios::out);
315  if (!(phsp_file))
316  egsFatal("\nEGS_PhspScoring: Failed to open phase space file %s for writing.\n",
317  phsp_fname.c_str());
318  //always MODE0 files--i.e. no ZLAST
319  phsp_file.write("MODE0",5);
320  //leave space for/skip over the rest of the header
321  phsp_file.seekp(28,ios::beg);
322  }
323  }
324  else if (oformat == 1) { //IAEA format
325  int iaea_iostat;
326  iaea_id = 1; //numerical index indicating this is the 1st file associated with this object scored
327  //hard coded to 1
328  if (is_restart) {
329  int rwmode = 3;
330  iaea_new_source(&iaea_id,phsp_fname_char,&rwmode,&iaea_iostat,len);
331  if (iaea_iostat < 0) {
332  egsFatal("\nEGS_PhspScoring: Failed to open phase space file %s.IAEAphsp for appending.\n",phsp_fname.c_str());
333  }
334  //check for consistency with total no. of scored particles as read from .egsdat file
335  EGS_I64 nparticle;
336  int type = -1;
337  int iaea_iostat;
338  iaea_get_max_particles(&iaea_id,&type,&nparticle);
339  if (nparticle != countprev) {
340  egsFatal("\nEGS_PhspScoring: Particle no. mismatch between %s.IAEAphsp and .egsdat file.\n",phsp_fname.c_str());
341  }
342  iaea_check_file_size_byte_order(&iaea_id,&iaea_iostat);
343  if (iaea_iostat != 0) {
344  egsFatal("\nEGS_PhspScoring: Byte order/file size mismatch in %s.IAEAphsp.\n",phsp_fname.c_str());
345  }
346  }
347  else {
348  int rwmode = 2;
349  iaea_new_source(&iaea_id,phsp_fname_char,&rwmode,&iaea_iostat,len);
350  if (iaea_iostat < 0) {
351  egsFatal("\nEGS_PhspScoring: Failed to open phase space file %s.IAEAphsp for writing.\n",phsp_fname.c_str());
352  }
353  }
354 
355 
356  //set up constant variables
357  for (int i=0; i<3; i++) {
358  if (xyz_is_const[i]) {
359  int index = i;
360  float constval = xyzscore[i];
361  iaea_set_constant_variable(&iaea_id,&index,&constval);
362  }
363  }
364  //set up extra floats and int indices and types
365  //need to store below in _tmp variables because this is a const function
366  int latch_ind_tmp = latch_ind;
367  int iaea_n_extra_long_tmp=iaea_n_extra_long;
368  int iaea_i_latch_tmp=iaea_i_latch;
369  int mu_ind_tmp = mu_ind;
370  int iaea_i_mu_tmp = iaea_i_mu;
371  int iaea_n_extra_float_tmp=iaea_n_extra_float;
372 
373  iaea_set_extra_numbers(&iaea_id,&iaea_n_extra_float_tmp,&iaea_n_extra_long_tmp);
374  iaea_set_type_extralong_variable(&iaea_id,&iaea_i_latch_tmp,&latch_ind_tmp);
375  if (score_mu) {
376  iaea_set_type_extrafloat_variable(&iaea_id,&iaea_i_mu_tmp,&mu_ind_tmp);
377  }
378  }
379 }
380 
381 //write phsp_index particles to phase space file
382 //update header and reset phsp_index
383 int EGS_PhspScoring::flushBuffer() const {
384 
385  if (first_flush) {
386  openPhspFile(); //here's where we open the phase space file
387  }
388  //awkward logic because we do not know if this is a restart
389  //until after initialization
390  first_flush = false;
391 
392  if (oformat == 1) { //iaea format
393  for (int j=0; j<phsp_index; j++) {
394  //fairly transparent, could probably put a lot of this in a separate method
395  //undo -ve energy marker and use n_stat to indicate new primary hist.
396  int n_stat = p_stack[j].E < 0 ? 1 : 0;
397  float E = abs(p_stack[j].E);
398  //convert charge to iaea type
399  int type = iaea_q_type[p_stack[j].q+1];
400 
401  //store latch in iaea_extra_long
402  EGS_I32 *iaea_extra_long = new EGS_I32[iaea_n_extra_long];
403  iaea_extra_long[iaea_i_latch]=p_stack[j].latch;
404  float *iaea_extra_float = new float[iaea_n_extra_float];
405  if (score_mu) {
406  iaea_extra_float[iaea_i_mu] = p_stack[j].mu;
407  }
408 
409  //now store double precision values in single precision reals
410  float wt = p_stack[j].wt;
411  float x = p_stack[j].x;
412  float y = p_stack[j].y;
413  float z = p_stack[j].z;
414  float u = p_stack[j].u;
415  float v = p_stack[j].v;
416  float w = p_stack[j].w;
417  //now actually write data
418 
419  iaea_write_particle(&iaea_id,&n_stat,&type,&E,&wt,&x,&y,&z,&u,&v,&w,iaea_extra_float,iaea_extra_long);
420  if (n_stat < 0) {
421  egsFatal("\nEGS_PhspScoring: Failed to write particle data to phase space file.");
422  }
423  }
424  //update header with no. of primary histories
425  EGS_I64 last_case_tmp = last_case;
426  iaea_set_total_original_particles(&iaea_id,&last_case_tmp);
427  //update header file
428  int iaea_iostat;
429  iaea_update_header(&iaea_id,&iaea_iostat);
430  if (iaea_iostat<0) {
431  egsFatal("\nEGS_PhspScoring: Failed to update phase space header file.");
432  }
433  }
434  else if (oformat == 0) { //EGSnrc format
435  if (!phsp_file) {
436  egsFatal("\nEGS_PhspScoring: phase space file is not open for writing.");
437  }
438  //don't forget that phsp_index is incremented after every particle written to p_stack
439  for (int j=0; j<phsp_index; j++) {
440  egs_phsp_write_struct ws(p_stack[j]);
441  phsp_file.write((char *) &ws, sizeof(ws));
442  }
443  //store position of end of file
444  iostream::off_type pos = (count+1)*sizeof(egs_phsp_write_struct);
445  //update header
446  phsp_file.seekp(5,ios::beg);
447  unsigned int count4 = count, countg4 = countg;
448  float pinc = last_case;
449  float emintmp;
450  if (countg4==count4) {
451  emintmp = 0.0;
452  }
453  else {
454  emintmp = emin;
455  }
456  phsp_file.write((char *) &count4, sizeof(unsigned int));
457  phsp_file.write((char *) &countg4, sizeof(unsigned int));
458  phsp_file.write((char *) &emax, sizeof(float));
459  phsp_file.write((char *) &emintmp, sizeof(float));
460  phsp_file.write((char *) &pinc, sizeof(float));
461  phsp_file.seekp(pos,ios::beg);
462  }
463 
464  phsp_index=0;
465 
466  return 0;
467 };
468 
469 bool EGS_PhspScoring::storeState(ostream &data) const {
470  if (!egsStoreI64(data,count)) {
471  return false;
472  }
473  if (!egsStoreI64(data,countg)) {
474  return false;
475  }
476  //update phase space file at the end of each batch
477  flushBuffer();
478  return true;
479 }
480 
481 bool EGS_PhspScoring::setState(istream &data) {
482  if (!egsGetI64(data,count)) {
483  return false;
484  }
485  if (!egsGetI64(data,countg)) {
486  return false;
487  }
488  countprev = count;
489  is_restart = true;
490  return true;
491 }
492 
493 bool EGS_PhspScoring::addState(istream &data) {
494  EGS_I64 tmp;
495  if (!egsGetI64(data,tmp)) { //return false;
496  egsWarning("error while reading count\n");
497  return false;
498  }
499  count += tmp;
500  if (!egsGetI64(data,tmp)) { //return false;
501  egsWarning("error while reading countg\n");
502  return false;
503  }
504  countg += tmp;
505  return true;
506 }
507 
508 //*********************************************************************
509 // Process input for this ausgab object
510 //**********************************************************************
511 extern "C" {
512 
513  EGS_PHSP_SCORING_EXPORT EGS_AusgabObject *createAusgabObject(EGS_Input *input,
514  EGS_ObjectFactory *f) {
515  const static char *func = "createAusgabObject(phsp_scoring)";
516  if (!input) {
517  egsWarning("%s: null input?\n",func);
518  return 0;
519  }
520  string str;
521  EGS_BaseGeometry *phspgeom;
522  int iscoremc = 0; //default to not score multiple crossers
523  vector <int> from_reg, to_reg;
524  int stype = 0; //default is to use scoring geom
525  int phspouttype;
526  int ptype;
527  int sdir=0;
528  int imuscore = 0;
529  float xyzconst[3];
530  bool xyzisconst[3] = {false, false, false};
531  string gname;
532  string outdir;
533  int err01 = input->getInput("phase space geometry",gname);
534  if (err01) {
535  stype = 1;
536  }
537  else {
538  phspgeom = EGS_BaseGeometry::getGeometry(gname);
539  if (!phspgeom) {
540  egsWarning("\nEGS_PhspScoring: %s does not name an existing geometry.\n"
541  "Will assume you want to use exit/entry region pairs.\n",gname.c_str());
542  stype = 1;
543  }
544  else {
545  if (input->getInput("score particles on", str) < 0) {
546  egsInformation("EGS_PhspScoring: No input for scoring direction.\n");
547  egsInformation("Will score on entry and exit from phase space geometry.\n");
548  sdir = 0;
549  }
550  else {
551  //get scoring direction
552  vector<string> allowed_sdir;
553  allowed_sdir.push_back("entry and exit");
554  allowed_sdir.push_back("entry");
555  allowed_sdir.push_back("exit");
556  sdir = input->getInput("score particles on",allowed_sdir,-1);
557  if (sdir < 0) {
558  egsFatal("\nEGS_PhspScoring: Invalid scoring direction.\n");
559  }
560  }
561  }
562  }
563  if (stype==1) {
564  // user wants to use exit/entry region pairs
565  int err05 = input->getInput("from regions",from_reg);
566  int err06 = input->getInput("to regions",to_reg);
567  if (err05 || err06) {
568  egsFatal("\nEGS_PhspScoring: Missing/incorrect input for scoring method\n"
569  "(scoring geometry or pairs of exit/entry regions)\n");
570  }
571  else {
572  //run some checks on exit/entry region pairs
573  vector<int>::iterator p,p1;
574  if (from_reg.size() > to_reg.size()) {
575  p = from_reg.begin();
576  egsWarning("\nEGS_PhspScoring: Mismatch in no. of exit/entry regions.\n"
577  "Will only score for matched pairs.\n");
578  p += to_reg.size();
579  from_reg.erase(p,p+from_reg.size()-to_reg.size());
580  }
581  else if (to_reg.size() > from_reg.size()) {
582  p = to_reg.begin();
583  egsWarning("\nEGS_PhspScoring: Mismatch in no. of exit/entry regions.\n"
584  "Will only score for matched pairs.\n");
585  p += from_reg.size();
586  to_reg.erase(p,p+to_reg.size()-from_reg.size());
587  }
588  //now go through and look for exit region = entry region
589  int i=0;
590  while (i<from_reg.size()) {
591  if (from_reg[i]==to_reg[i]) {
592  egsInformation("\nEGS_PhspScoring: Cannot have entry region = exit region (reg no. %d)\n",from_reg[i]);
593  egsInformation("Will delete this pair\n");
594  p=from_reg.begin()+i;
595  p1 =to_reg.begin()+i;
596  from_reg.erase(p);
597  to_reg.erase(p1);
598  }
599  else {
600  //advance counter
601  i++;
602  }
603  }
604  }
605  }
606  //now get common inputs for both scoring methods
607  if (input->getInput("output format", str) < 0) {
608  egsInformation("EGS_PhspScoring: No input for output format type. Will default to EGSnrc.\n");
609  phspouttype = 0;
610  }
611  else {
612  vector<string> allowed_oformat;
613  allowed_oformat.push_back("EGSnrc");
614  allowed_oformat.push_back("IAEA");
615  phspouttype = input->getInput("output format", allowed_oformat, -1);
616  if (phspouttype < 0) {
617  egsFatal("\nEGS_PhspScoring: Invalid output format.\n");
618  }
619  //see if the user wants to specify constant X/Y/Z for IAEA format
620  if (phspouttype == 1) {
621  int err02 = input->getInput("constant X",xyzconst[0]);
622  int err03 = input->getInput("constant Y",xyzconst[1]);
623  int err04 = input->getInput("constant Z",xyzconst[2]);
624  if (!err02) {
625  xyzisconst[0] = true;
626  }
627  if (!err03) {
628  xyzisconst[1] = true;
629  }
630  if (!err04) {
631  xyzisconst[2] = true;
632  }
633  //see if user wants to score mu (if available)
634  //default is not to score
635  if (!input->getInput("score mu", str)) {
636  vector<string> allowed_muscore;
637  allowed_muscore.push_back("no");
638  allowed_muscore.push_back("yes");
639  imuscore = input->getInput("score mu",allowed_muscore,-1);
640  if (imuscore < 0) {
641  egsWarning("\nEGS_PhspScoring: Invalid input for mu scoring. Will not score mu.\n");
642  imuscore = 0;
643  }
644  }
645  }
646  }
647  if (phspouttype == 0) {
648  //see if user wants to score multiple crossers
649  if (!input->getInput("score multiple crossers", str)) {
650  vector<string> allowed_scoremc;
651  allowed_scoremc.push_back("no");
652  allowed_scoremc.push_back("yes");
653  iscoremc = input->getInput("score multiple crossers",allowed_scoremc,-1);
654  if (iscoremc < 0) {
655  egsWarning("\nEGS_PhspScoring: Invalid input for score multiple crossers. Will not score.\n");
656  iscoremc = 0;
657  }
658  }
659  }
660  if (input->getInput("output directory",outdir) < 0) {
661  outdir="";
662  }
663  if (input->getInput("particle type", str) < 0) {
664  egsInformation("EGS_PhspScoring: No input for particle type. Will score all.\n");
665  ptype = 0;
666  }
667  else {
668  //get particle type
669  vector<string> allowed_ptype;
670  allowed_ptype.push_back("all");
671  allowed_ptype.push_back("photons");
672  allowed_ptype.push_back("charged");
673  ptype = input->getInput("particle type",allowed_ptype,-1);
674  if (ptype < 0) {
675  egsFatal("\nEGS_PhspScoring: Invalid particle type.\n");
676  }
677  }
678 
679  //=================================================
680 
681  /* Setup phsp scoring object with input parameters */
682  EGS_PhspScoring *result = new EGS_PhspScoring("",f);
683  result->setName(input);
684  if (stype==0) {
685  result->setGeom(phspgeom);
686  }
687  else if (stype==1) {
688  result->setEntryExitReg(from_reg,to_reg);
689  }
690  result->setOType(phspouttype);
691  result->setXYZconst(xyzisconst,xyzconst);
692  result->setOutDir(outdir);
693  result->setParticleType(ptype);
694  result->setScoreDir(sdir);
695  result->setMuScore(imuscore);
696  result->setScoreMC(iscoremc);
697  return result;
698  }
699 }
A phase space scoring object: header.
bool EGS_EXPORT egsGetI64(istream &data, EGS_I64 &n)
Reads a 64 bit integer from the stream data and assigns it to n. Returns true on success, false on failure.
EGS_Float E
particle energy in MeV
string description
A short ausgab object description.
const string & getObjectName() const
Get the object name.
int q
particle charge
EGS_Input class header file.
EGS_Float y
y-component
Definition: egs_vector.h:61
bool EGS_EXPORT egsStoreI64(ostream &data, EGS_I64 n)
Writes the 64 bit integer n to the output stream data and returns true on success, false on failure.
string egsJoinPath(const string &first, const string &second)
Join two path variables (or a path and a file name) using the platform specific directory separator a...
int getIparallel() const
Returns the job number in a parallel run.
Global egspp functions header file.
EGS_Vector x
position
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
int getNparallel() const
Returns the number of parallel jobs executing.
EGS_Float z
z-component
Definition: egs_vector.h:62
static EGS_BaseGeometry * getGeometry(const string &Name)
Get a pointer to the geometry named Name.
EGS_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
A phase space scoring ausgab object.
EGS_InfoFunction EGS_EXPORT egsInformation
Always use this function for reporting the progress of a simulation and any other type of information...
An object factory.
void setName(EGS_Input *inp)
Set the name of the object from the information provided by inp.
EGS_Float x
x-component
Definition: egs_vector.h:60
int latch
latch variable (useful as a flag on many occasions)
EGS_Vector u
direction
string name
The object name.
EGS_Float getMU()
Returns the value of the mu synchronization parameter.
A class for storing information in a tree-like structure of key-value pairs. This class is used throu...
Definition: egs_input.h:182
const string & getName() const
Get the name of this geometry.
const string & getAppDir() const
Returns the absolute path to the user code directory.
EGS_Particle top_p
The top particle on the stack (i.e., the particle being transported)
string otype
The object type.
EGS_Float wt
statistical weight
virtual void setApplication(EGS_Application *App)
Set the application this object belongs to.
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
Base class for advanced EGSnrc C++ applications.
EGS_Application * app
The application this object belongs to.
EGS_InfoFunction EGS_EXPORT egsWarning
Always use this function for reporting warnings.