EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
iaea_phsp_source.cpp
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ IAEA format phase-space source
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, 2013
25 #
26 # Contributors: Reid Townson
27 # Hubert Ho
28 #
29 ###############################################################################
30 */
31 
32 
38 #include "iaea_phsp.h"
39 #include "iaea_phsp_source.h"
40 #include "egs_input.h"
41 #include "egs_functions.h"
42 
43 IAEA_PhspSource::IAEA_PhspSource(const string &phsp_file,
44  const string &Name, EGS_ObjectFactory *f) : EGS_BaseSource(Name,f) {
45  init();
46  openFile(phsp_file);
47 }
48 
49 void IAEA_PhspSource::init() {
50  otype = "IAEA_PhspSource";
51  Xmin = -veryFar;
52  Xmax = veryFar;
53  Ymin = -veryFar;
54  Ymax = veryFar;
55  is_valid = false;
56  mode2 = false;
57  swap_bytes = false;
58  the_file_name = "no file";
59  filter_type = -1;
60  particle_type = 2;
61  description = "Invalid IAEA phase space source";
62  Nread = 0;
63  count = 0;
64  Nrestart = 0;
65  Npos = 0;
66  Nlast = 0;
67  wmin = -veryFar;
68  wmax = veryFar;
69  Nrecycle_g = 0;
70  Nrecycle_e = 0;
71  Nrecycle = 0;
72  Nuse = -1;
73  first = true;
74 }
75 
76 void IAEA_PhspSource::openFile(const string &phsp_file) {
77  the_file_name = "no file";
78  is_valid = false;
79  int rwmode=1; //read only
80  iaea_fileid=0; //this is just an array index for the iaea routines, not an actual
81  //unit no., as in Fortran, so it can really be any integer value
82  //not above MAX_NUM_SOURCES (30)
83 
84  //below is required because iaea opening routine gets file name as char array
85  int len = phsp_file.length();
86  char phsp_name_tmp[len+1];
87  phsp_file.copy(phsp_name_tmp,len,0);
88  //now put a null character on the end
89  phsp_name_tmp[len]='\0';
90 
91  //use iaea function for opening phase space file
92  iaea_new_source(&iaea_fileid,phsp_name_tmp,&rwmode,&iaea_iostat,len);
93  if (iaea_iostat==105) {
94  egsWarning("IAEA_PhspSource::openFile: no header file name supplied.\n");
95  return;
96  }
97  else if (iaea_iostat==-96) {
98  egsWarning("IAEA_PhspSource::openFile: failed to open header file %s.IAEAheader\n"
99  " for reading\n",phsp_file.c_str());
100  return;
101  }
102  else if (iaea_iostat==-94) {
103  egsWarning("IAEA_PhspSource::openFile: failed to open phase space file %s.IAEAphsp\n"
104  " for reading\n",phsp_file.c_str());
105  return;
106  }
107  else if (iaea_iostat==-1) {
108  egsWarning("IAEA_PhspSource::openFile: failed to initialize phase space source %s\n",phsp_file.c_str());
109  return;
110  }
111  else if (iaea_iostat==-91) {
112  egsWarning("IAEA_PhspSource::openFile: failed to get record contents from header file %s.IAEAheader\n",phsp_file.c_str());
113  return;
114  }
115  else if (iaea_iostat <0) {
116  egsWarning("IAEA_PhspSource::openFile: I/O error ocurred on opening phase space file %s\n",phsp_file.c_str());
117  return;
118  }
119 
120  //now check file size and byte order
121 
122  swap_bytes = false;
123  iaea_check_file_size_byte_order(&iaea_fileid,&iaea_iostat);
124  if (iaea_iostat==-1) {
125  egsWarning("IAEA_PhspSource::openFile: error reading file size/byte order from header %s.IAEAheader\n",phsp_file.c_str());
126  return;
127  }
128  else if (iaea_iostat==-4) {
129  egsWarning("IAEA_PhspSource::openFile: byte mismatch between phase space file and machine. Will swap bytes.\n");
130  swap_bytes = true;
131  }
132  else if (iaea_iostat==-3 || iaea_iostat==-5) {
133  egsWarning("IAEA_PhspSource::openFile: mismatch between file size in header and actual file size of %s\n",phsp_file.c_str());
134  return;
135  }
136 
137  //now get some info from the header
138  EGS_I64 n, n_photon, pinc;
139  float emax;
140  int iaea_type=-1; //for getting no. of particles
141  iaea_get_max_particles(&iaea_fileid,&iaea_type,&n);
142  if (n<0) {
143  egsWarning("IAEA_PhspSource::openFile: failed to get total no. of particles from %s.IAEAheader\n",phsp_file.c_str());
144  return;
145  }
146  iaea_type=1; //for getting no. of photons
147  iaea_get_max_particles(&iaea_fileid,&iaea_type,&n_photon);
148  if (n_photon<0) {
149  egsWarning("IAEA_PhspSource::openFile: failed to get no. of photons from %s.IAEAheader\n",phsp_file.c_str());
150  return;
151  }
152  //now get max. energy
153  iaea_get_maximum_energy(&iaea_fileid,&emax);
154  if (emax<0) {
155  egsWarning("IAEA_PhspSource::openFile: failed to get max. energy from %s.IAEAheader\n",phsp_file.c_str());
156  return;
157  }
158  //no. of particles incident from original source
159  iaea_get_total_original_particles(&iaea_fileid,&pinc);
160  if (pinc<0) {
161  egsWarning("IAEA_PhspSource::openFile: failed to get no. of original particles from %s.IAEAheader\n",phsp_file.c_str());
162  return;
163  }
164 
165  //get number of extra float and extra long variables so we can set array dimensions
166  iaea_get_extra_numbers(&iaea_fileid,&n_extra_floats,&n_extra_longs);
167  if (n_extra_floats==-1 || n_extra_longs==-1) {
168  egsWarning("IAEA_PhspSource::openFile: failed to get no. of extra floats and longs stored in %s.IAEAphsp\n",phsp_file.c_str());
169  return;
170  }
171  //determine if Zlast is stored in the file and, if so, its array index, i_zlast
172  int extrafloat_types[MAXEXTRAS], extralong_types[MAXEXTRAS];
173  iaea_get_type_extra_variables(&iaea_fileid,&iaea_iostat,extralong_types,extrafloat_types);
174  if (iaea_iostat==-1) {
175  egsWarning("IAEA_PhspSource::openFile: failed to get Mode of data %s.IAEAheader\n",phsp_file.c_str());
176  return;
177  }
178  mode2=false;
179  i_zlast = -1;
180  for (int i=0; i< n_extra_floats; i++) {
181  if (extrafloat_types[i]==3) {
182  mode2=true;
183  i_zlast=i;
184  break;
185  }
186  }
187  //now see if mu index is stored in the file
188  //assume it is the first variable stored in extra_float
189  //after zlast
190  i_mu=-1;
191  mu_stored=false;
192  if (n_extra_floats > i_zlast+1) {
193  if (extrafloat_types[i_zlast+1] == 0) {
194  i_mu=i_zlast+1;
195  mu_stored=true;
196  egsInformation("IAEA_PhspSource::openFile: Mu index included in data in %s.IAEAphsp\n",phsp_file.c_str());
197  }
198  }
199 
200  //determine array index of latch
201  latch_stored=false;
202  for (int i=0; i< n_extra_longs; i++) {
203  if (extralong_types[i]==2) {
204  latch_stored=true;
205  i_latch=i;
206  break;
207  }
208  }
209  if (!latch_stored) {
210  egsWarning("IAEA_PhspSource::openFile: LATCH is not stored in data in %s.IAEAphsp\n",phsp_file.c_str());
211  }
212 
213  Npos = 0;
214  Nlast = n;
215  Nfirst = 1;
216  // at this point the position should be at the first particle in the file
217  Emax = emax;
218  Pinc = pinc;
219  Nparticle = n;
220  Nphoton = n_photon;
221  is_valid = true;
222  the_file_name = phsp_file;
223 }
224 
226  EGS_BaseSource(input,f) {
227  init();
228  string fname;
229  int err = input->getInput("iaea phase space file",fname);
230  if (err) {
231  egsWarning("IAEA_PhspSource: no 'iaea phase space file' input\n");
232  return;
233  }
234  openFile(fname);
235  if (!isValid()) {
236  egsWarning("IAEA_PhspSource: errors while opening the phase space file"
237  " %s\n",fname.c_str());
238  return;
239  }
240  vector<EGS_Float> cutout;
241  err = input->getInput("cutout",cutout);
242  if (!err && cutout.size() == 4) {
243  setCutout(cutout[0],cutout[1],cutout[2],cutout[3]);
244  }
245  vector<string> ptype;
246  ptype.push_back("electrons");
247  ptype.push_back("photons");
248  ptype.push_back("positrons");
249  ptype.push_back("all");
250  ptype.push_back("charged");
251  particle_type = input->getInput("particle type",ptype,3)-1;
252  vector<int> the_filter;
253  vector<string> the_latch;
254  err = input->getInput("filter type",the_filter);
255  int err1 = input->getInput("latch setting",the_latch);
256  if (!err && !err1) {
257  if (the_filter[0] >= 0 && the_filter[0] <= 3) {
258  int nbit1 = the_latch[0].size();
259  int nbit2 = 0;
260  if (!the_filter[0]) {
261  nbit2 = the_latch[1].size();
262  }
263  if (nbit1 + nbit2 > 0) {
264  int *the_bits = new int [nbit1+nbit2];
265  for (int j=0; j<nbit1+nbit2; j++) {
266  if (j < nbit1) {
267  the_bits[j]=(int)(the_latch[0][j] == '1');
268  }
269  else {
270  the_bits[j] = (int)(the_latch[1][j-nbit1] == '1');
271  }
272  }
273  setFilter(the_filter[0],nbit1,nbit2,the_bits);
274  delete [] the_bits;
275  }
276  }
277  }
278  vector<EGS_Float> wwindow;
279  err = input->getInput("weight window",wwindow);
280  if (!err && wwindow.size() == 2) {
281  wmin = wwindow[0];
282  wmax = wwindow[1];
283  }
284  int ntmp;
285  err = input->getInput("reuse photons",ntmp);
286  if (!err && ntmp > 0) {
287  Nrecycle_g = ntmp;
288  }
289  else {
290  err = input->getInput("recycle photons",ntmp);
291  if (!err && ntmp > 0) {
292  Nrecycle_g = ntmp;
293  }
294  }
295  err = input->getInput("reuse electrons",ntmp);
296  if (!err && ntmp > 0) {
297  Nrecycle_e = ntmp;
298  }
299  else {
300  err = input->getInput("recycle electrons",ntmp);
301  if (!err && ntmp > 0) {
302  Nrecycle_e = ntmp;
303  }
304  }
305  description = "IAEA phase space source from ";
307 }
308 
309 EGS_I64 IAEA_PhspSource::getNextParticle(EGS_RandomGenerator *, int &q,
310  int &latch, EGS_Float &E, EGS_Float &wt, EGS_Vector &x, EGS_Vector &u) {
311  /*
312  if( Nuse >= Nrecycle ) {
313  do { readParticle(); } while ( rejectParticle() );
314  }
315  */
316  int nstat,extrainttemp[MAXEXTRAS];
317  float extrafloattemp[MAXEXTRAS];
318  if (Nuse > Nrecycle || Nuse < 0) { //get a new particle
319  if ((++Npos) > Nlast) {
320  egsWarning("IAEA_PhspSource::getNextParticle(): reached the end of the "
321  "phase space file chunk (%lld)\n will start from the beginning "
322  "of the chunk (%lld) but this "
323  "implies that uncertainty estimates will be inaccurate\n",
324  Nlast,Nfirst);
325  iaea_set_record(&iaea_fileid,&Nfirst,&iaea_iostat);
326  if (iaea_iostat<0) {
327  egsFatal("IAEA_PhspSource::getNextParticle(): error restarting phase space chunk\n");
328  }
329  Nrestart++;
330  Npos = Nfirst;
331  }
332  iaea_get_particle(&iaea_fileid,&nstat,&p.q,&p.E,&p.wt,&p.x,&p.y,&p.z,&p.u,&p.v,&p.w,extrafloattemp,extrainttemp);
333  ++Nread;
334  p.latch=0; //important if we are using latch to do vr
335  /*
336  if (latch_stored) {
337  p.latch = extrainttemp[i_latch];
338  }
339  */
340  if (mode2) {
341  p.zlast = extrafloattemp[i_zlast];
342  }
343  if (mu_stored) {
344  p.mu = extrafloattemp[i_mu];
345  }
346  if (swap_bytes) {
347  egsSwapBytes(&p.q);
348  egsSwapBytes(&nstat);
349  egsSwapBytes(&p.zlast);
350  egsSwapBytes(&p.latch);
351  egsSwapBytes(&p.E);
352  egsSwapBytes(&p.wt);
353  egsSwapBytes(&p.x);
354  egsSwapBytes(&p.y);
355  egsSwapBytes(&p.z);
356  egsSwapBytes(&p.u);
357  egsSwapBytes(&p.v);
358  egsSwapBytes(&p.w);
359  egsSwapBytes(&p.mu);
360  }
361  //do check here because we need the swapped version of nstat
362  if (nstat<0) {
363  egsFatal("IAEA_PhspSource::getNextParticle(): error reading particle number %i\n",Npos);
364  }
365  //convert charge from iaea type
366  if (p.q==1) {
367  p.q=0;
368  }
369  else if (p.q==2) {
370  p.q=-1;
371  }
372  else if (p.q==3) {
373  p.q=1;
374  }
375  else {
376  egsFatal("IAEA_PhspSource::getNextParticle: unknown charge on particle\n");
377  }
378  //note: we don't have to convert to p.E to K.E. because that's what IAEA format stores
379  if (first || nstat>0) {
380  count++; //increment primary history counter
381  }
382  first= false;
383  //store Nrecycle
384  if (p.q) {
386  }
387  else {
389  }
390  //reset Nuse
391  Nuse = 0;
392  p.wt /= (Nrecycle+1);
393  }
394 
395  //energy, wt, position and direction cosines
396  x.x = p.x;
397  x.y = p.y;
398  x.z = p.z;
399  u.x = p.u;
400  u.y = p.v, u.z = p.w;
401  wt=p.wt;
402  if (rejectParticle()) {
403  wt = 0;
404  }
405  E=p.E;
406  q=p.q;
407  //latch = 0;
408  latch= p.latch;
409  //rejectParticle uses BeamParticle, p
410  ++Nuse;
411  return count;
412 }
413 
414 #ifndef SKIP_DOXYGEN
415 
419 union __egs_data32 {
420  int i;
421  float f;
422  char c[4];
423 };
424 #endif
425 
426 void IAEA_PhspSource::setSimulationChunk(EGS_I64 nstart, EGS_I64 nrun) {
427  if (nstart < 0 || nrun < 1 || nstart + nrun > Nparticle) {
428  egsWarning("IAEA_PhspSource::setSimulationChunk(): illegal attempt "
429  "to set the simulation chunk between %lld and %lld ignored\n",
430  nstart+1,nstart+nrun);
431  return;
432  }
433  Nfirst = nstart+1;
434  Nlast = nstart + nrun;
435  Npos = nstart;
436  iaea_set_record(&iaea_fileid,&Nfirst,&iaea_iostat);
437  if (iaea_iostat<0) {
438  egsWarning("IAEA_PhspSource::setSimulationChunk(): error setting phase space chunk\n");
439  return;
440  }
441  egsInformation("IAEA_PhspSource: using phsp portion between %lld and %lld\n",
442  Nfirst,Nlast);
443 }
444 
445 bool IAEA_PhspSource::rejectParticle() const {
446  if (particle_type < 2 && p.q != particle_type) {
447  return true;
448  }
449  if (particle_type == 3 && !p.q) {
450  return true;
451  }
452  if (p.x < Xmin || p.x > Xmax || p.y < Ymin || p.y > Ymax) {
453  return true;
454  }
455  if (p.wt < wmin || p.wt > wmax) {
456  return true;
457  }
458  if (p.latch & 2147483648UL) {
459  return true;
460  }
461  if (filter_type < 0) {
462  return false;
463  }
464  if (filter_type == 1) {
465  bool r1;
466  if (filter1) {
467  r1 = !(p.latch & filter1);
468  }
469  else {
470  r1 = false;
471  }
472  bool r2 = (p.latch & filter2);
473  return (r1 || r2);
474  }
475  if (filter_type == 1) {
476  return (p.latch & filter1);
477  }
478  int l = p.latch >> 24;
479  bool res = l & filter1;
480  return filter_type == 3 ? res : !res;
481 }
482 
483 void IAEA_PhspSource::setFilter(int type, int nbit1, int nbit2, const int *bits) {
484  if (type < 0 || type > 3) {
485  egsWarning("IAEA_PhspSource::setFilter: invalid filter type %d\n",type);
486  return;
487  }
488  int ntot = nbit1;
489  if (type == 0) {
490  ntot += nbit2;
491  }
492  if (ntot > 29) {
493  egsWarning("IAEA_PhspSource::setFilter: maximum number of bits is "
494  "limited to 29, you requested %d\n",ntot);
495  return;
496  }
497  if (!ntot) {
498  return;
499  }
500  filter_type = type;
501  if (filter_type == 0) {
502  int i;
503  filter1 = 0;
504  filter2 = 0;
505  for (i=0; i<ntot; i++) {
506  if (bits[i]) {
507  unsigned int aux = 1;
508  aux = (aux << i);
509  if (i < nbit1) {
510  filter1 += aux;
511  }
512  else {
513  filter2 += aux;
514  }
515  }
516  }
517  }
518  else if (filter_type == 1) { // TODO
519  }
520  else if (filter_type == 2) { // TODO
521  }
522  else if (filter_type == 3) { // TODO
523  }
524 }
525 
526 extern "C" {
527 
528  IAEA_PHSP_SOURCE_EXPORT EGS_BaseSource *createSource(EGS_Input *input,
529  EGS_ObjectFactory *f) {
530  return
531  createSourceTemplate<IAEA_PhspSource>(input,f,"iaea phsp source");
532  }
533 
534 }
EGS_I64 Nread
Number of particles read so far.
int n_extra_floats
no. of extra floats stored in phsp file
EGS_I64 Nfirst
first record this source can use
int i_zlast
index of zlast in extra_floats array
int i_mu
index of mu index in extra_floats array
bool mode2
true, if a MODE2 file (i.e. storing Zlast)
int Nrecycle_e
Number of times to recycle a charged particle.
EGS_Input class header file.
EGS_Float y
y-component
Definition: egs_vector.h:61
A class representing 3D vectors.
Definition: egs_vector.h:56
int Nrestart
Number of times the file was restarted.
Global egspp functions header file.
EGS_I64 Npos
Next record to be read.
An IAEA format phase-space file source.
const EGS_Float veryFar
A very large float.
void egsSwapBytes(int *n)
Swap the bytes of 32 bit integers.
int iaea_iostat
iostat on read/write of iaea phsp file
int Nuse
Number of times current particle was used so far.
Base random number generator class. All random number generators should be derived from this class...
Definition: egs_rndm.h:67
EGS_Float z
z-component
Definition: egs_vector.h:62
EGS_I64 Nparticle
Number of particles in the file.
EGS_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
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.
EGS_I64 Nphoton
Number of photons in the file.
EGS_Float Pinc
Number of incident particles that created the file.
EGS_Float Emax
Maximum k.e. (obtained from the phsp file)
int iaea_fileid
phsp file unit no.
EGS_Float x
x-component
Definition: egs_vector.h:60
int i_latch
index of latch in extra_floats array
EGS_I64 Nlast
Last record this source can use.
bool mu_stored
true if mu index stored
string the_file_name
The phase-space file name.
A class for storing information in a tree-like structure of key-value pairs. This class is used throu...
Definition: egs_input.h:182
int Nrecycle_g
Number of times to recycle a photon.
string otype
The object type.
bool latch_stored
true if LATCH is stored in data
int Nrecycle
Number of times to recycle current particle.
IAEA_PhspSource(const string &phsp_file, const string &Name="", EGS_ObjectFactory *f=0)
Constructor.
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 source class. All particle sources must be derived from this class.
string description
A short source description.
EGS_InfoFunction EGS_EXPORT egsWarning
Always use this function for reporting warnings.
int n_extra_longs
no. of extra longs stored in phsp file