EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_phsp_source.cpp
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ 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: Iwan Kawrakow, 2005
25 #
26 # Contributors: Ernesto Mainegra-Hing
27 # Frederic Tessier
28 # Reid Townson
29 #
30 ###############################################################################
31 */
32 
33 
39 #include "egs_phsp_source.h"
40 #include "egs_input.h"
41 #include "egs_functions.h"
42 #include "egs_application.h"
43 
44 EGS_PhspSource::EGS_PhspSource(const string &phsp_file,
45  const string &Name, EGS_ObjectFactory *f) : EGS_BaseSource(Name,f) {
46  init();
47  openFile(phsp_file);
48 }
49 
50 void EGS_PhspSource::init() {
51  otype = "EGS_PhspSource";
52  recl = 0;
53  Xmin = -veryFar;
54  Xmax = veryFar;
55  Ymin = -veryFar;
56  Ymax = veryFar;
57  is_valid = false;
58  record = 0;
59  mode2 = false;
60  swap_bytes = false;
61  the_file_name = "no file";
62  filter_type = -1;
63  particle_type = 2;
64  description = "Invalid phase space source";
65  Nread = 0;
66  count = 0;
67  Nrestart = 0;
68  Npos = 0;
69  Nlast = 0;
70  wmin = -veryFar;
71  wmax = veryFar;
72  Nrecycle_g = 0;
73  Nrecycle_e = 0;
74  Nrecycle = 0;
75  Nuse = -1;
76  first = true;
77 }
78 
79 void EGS_PhspSource::openFile(const string &phsp_file) {
80  if (the_file.is_open()) {
81  the_file.close();
82  }
83  the_file_name = "no file";
84  is_valid = false;
85  recl = 0;
86  if (record) {
87  delete [] record;
88  record = 0;
89  }
90  the_file.open(phsp_file.c_str(),ios::binary | ios::in);
91  if (!the_file.is_open()) {
92  egsWarning("EGS_PhspSource::openFile: failed to open binary file %s"
93  " for reading\n",phsp_file.c_str());
94  return;
95  }
96  string cmode;
97  char auxc;
98  for (int i=0; i<5; i++) {
99  the_file.get(auxc);
100  if (the_file.eof() || !the_file.good()) {
101  egsWarning("EGS_PhspSource::openFile: an I/O error occured "
102  "while reading the first record of %s\n",phsp_file.c_str());
103  return;
104  }
105  cmode += auxc;
106  }
107  if (cmode == "MODE0") {
108  mode2 = false;
109  recl = 28;
110  }
111  else if (cmode == "MODE2") {
112  mode2 = true;
113  recl = 32;
114  }
115  else {
116  egsWarning("EGS_PhspSource::openFile: the file %s is not a MODE0 or"
117  " MODE2 file\n",phsp_file.c_str());
118  return;
119  }
120  record = new char [recl];
121  int n, n_photon;
122  float emax, emin, pinc;
123  the_file.read((char *) &n, sizeof(int));
124  the_file.read((char *) &n_photon, sizeof(int));
125  the_file.read((char *) &emax, sizeof(float));
126  the_file.read((char *) &emin, sizeof(float));
127  the_file.read((char *) &pinc, sizeof(float));
128  if (the_file.eof() || !the_file.good()) {
129  egsWarning("EGS_PhspSource::openFile: an I/O error occured "
130  "while reading the first record of %s\n",phsp_file.c_str());
131  recl = 0;
132  return;
133  }
134  swap_bytes = false;
135  if (n <= 0 || n_photon < 0 || n_photon > n || emin < 0 || emax < 0 ||
136  emax < emin || pinc < 0) {
137  swap_bytes = true;
138  egsSwapBytes(&n);
139  egsSwapBytes(&n_photon);
140  egsSwapBytes(&emin);
141  egsSwapBytes(&emax);
142  egsSwapBytes(&pinc);
143  if (n <= 0 || n_photon < 0 || n_photon > n || emin < 0 || emax < 0 ||
144  emax < emin || pinc < 0) {
145  egsWarning("EGS_PhspSource::openFile: phase space file header"
146  " contains meaningless values with and without byte swaping:\n");
147  if (n <= 0) {
148  egsWarning(" number of particles: %d\n",n);
149  }
150  if (n_photon < 0) {
151  egsWarning(" number of photons: %d\n",n_photon);
152  }
153  if (n_photon > n) egsWarning(" number of photons (%d) is "
154  "greater than number of particles (%d)\n",n_photon,n);
155  if (emin < 0) {
156  egsWarning(" minimum energy: %g\n",emin);
157  }
158  if (emax < 0) {
159  egsWarning(" maximum energy: %g\n",emax);
160  }
161  if (emin > emax) {
162  egsWarning(" emin > emax: %g %g\n",emin,emax);
163  }
164  if (pinc < 0) {
165  egsWarning(" incident particles: %g\n",pinc);
166  }
167  recl = 0;
168  return;
169  }
170  }
171  // at this points we have passed a set of checks and think that
172  // we have some meaningful information about number of particles, etc.
173  // to be completely sure, it is a good idea to read the last particle
174  // in the file and check for errors.
175  istream::off_type nend = n;
176  nend = nend*recl;
177  //the_file.seekg(recl*n,ios::beg);
178  the_file.seekg(nend,ios::beg);
179  the_file.read(record,recl*sizeof(char));
180  if (the_file.bad() || the_file.fail()) {
181  egsWarning("EGS_PhspSource::openFile: failed to read the last"
182  " particle in the file, this indicates some unknown error condition\n");
183  recl = 0;
184  return;
185  }
186  the_file.clear();
187  the_file.seekg(recl,ios::beg);
188  Npos = 0;
189  Nlast = n;
190  Nfirst = 1;
191  // at this point the position should be at the first particle in the file
192  Emax = emax;
193  Emin = emin;
194  Pinc = pinc;
195  Nparticle = n;
196  Nphoton = n_photon;
197  is_valid = true;
198  the_file_name = phsp_file;
199 }
200 
202  EGS_BaseSource(input,f) {
203  init();
204  string fname;
205  int err = input->getInput("phase space file",fname);
206  if (err) {
207  egsWarning("EGS_PhspSource: no 'phase space file' input\n");
208  return;
209  }
210  openFile(fname);
211  if (!isValid()) {
212  egsWarning("EGS_PhspSource: errors while opening the phase space file"
213  " %s\n",fname.c_str());
214  return;
215  }
216  vector<EGS_Float> cutout;
217  err = input->getInput("cutout",cutout);
218  if (!err && cutout.size() == 4) {
219  setCutout(cutout[0],cutout[1],cutout[2],cutout[3]);
220  }
221  vector<string> ptype;
222  ptype.push_back("electrons");
223  ptype.push_back("photons");
224  ptype.push_back("positrons");
225  ptype.push_back("all");
226  ptype.push_back("charged");
227  particle_type = input->getInput("particle type",ptype,3)-1;
228  vector<int> the_filter;
229  vector<string> the_latch;
230  err = input->getInput("filter type",the_filter);
231  int err1 = input->getInput("latch setting",the_latch);
232  if (!err && !err1) {
233  if (the_filter[0] >= 0 && the_filter[0] <= 3) {
234  int nbit1 = the_latch[0].size();
235  int nbit2 = 0;
236  if (!the_filter[0]) {
237  nbit2 = the_latch[1].size();
238  }
239  if (nbit1 + nbit2 > 0) {
240  int *the_bits = new int [nbit1+nbit2];
241  for (int j=0; j<nbit1+nbit2; j++) {
242  if (j < nbit1) {
243  the_bits[j]=(int)(the_latch[0][j] == '1');
244  }
245  else {
246  the_bits[j] = (int)(the_latch[1][j-nbit1] == '1');
247  }
248  }
249  setFilter(the_filter[0],nbit1,nbit2,the_bits);
250  delete [] the_bits;
251  }
252  }
253  }
254  vector<EGS_Float> wwindow;
255  err = input->getInput("weight window",wwindow);
256  if (!err && wwindow.size() == 2) {
257  wmin = wwindow[0];
258  wmax = wwindow[1];
259  }
260  int ntmp;
261  err = input->getInput("reuse photons",ntmp);
262  if (!err && ntmp > 0) {
263  Nrecycle_g = ntmp;
264  }
265  else {
266  err = input->getInput("recycle photons",ntmp);
267  if (!err && ntmp > 0) {
268  Nrecycle_g = ntmp;
269  }
270  }
271  err = input->getInput("reuse electrons",ntmp);
272  if (!err && ntmp > 0) {
273  Nrecycle_e = ntmp;
274  }
275  else {
276  err = input->getInput("recycle electrons",ntmp);
277  if (!err && ntmp > 0) {
278  Nrecycle_e = ntmp;
279  }
280  }
281  description = "Phase space source from ";
283 }
284 
285 EGS_I64 EGS_PhspSource::getNextParticle(EGS_RandomGenerator *, int &q,
286  int &latch, EGS_Float &E, EGS_Float &wt, EGS_Vector &x, EGS_Vector &u) {
287  if (!recl) egsFatal("EGS_PhspSource::readParticle(): the file is not "
288  "open yet\n");
289  /*
290  if( Nuse >= Nrecycle ) {
291  do { readParticle(); } while ( rejectParticle() );
292  }
293  */
294  if (Nuse > Nrecycle || Nuse < 0) {
295  readParticle();
296  }
297  x.x = p.x;
298  x.y = p.y;
299  x.z = 0;
300  u.x = p.u;
301  u.y = p.v;
302  EGS_Float aux = p.u*p.u+p.v*p.v;
303  if (aux < 1) {
304  aux = sqrt(1-aux);
305  }
306  else {
307  aux = 0;
308  }
309  if (p.wt > 0) {
310  u.z = aux;
311  wt = p.wt;
312  }
313  else {
314  u.z = -aux;
315  wt = -p.wt;
316  }
317  if (rejectParticle()) {
318  wt = 0;
319  }
320  E = p.E;
321  q = p.q;
322  latch = 0; //latch = p.latch;
323  ++Nuse;
324  return count;
325 }
326 
327 #ifndef SKIP_DOXYGEN
328 
332 union __egs_data32 {
333  int i;
334  float f;
335  char c[4];
336 };
337 #endif
338 
339 void EGS_PhspSource::setSimulationChunk(EGS_I64 nstart, EGS_I64 nrun) {
340  if (nstart < 0 || nrun < 1 || nstart + nrun > Nparticle) {
341  egsWarning("EGS_PhspSource::setSimulationChunk(): illegal attempt "
342  "to set the simulation chunk between %lld and %lld ignored\n",
343  nstart+1,nstart+nrun);
344  return;
345  }
346  Nfirst = nstart+1;
347  Nlast = nstart + nrun;
348  Npos = nstart;
349  istream::off_type pos = Nfirst*recl;
350  the_file.seekg(pos,ios::beg);
351  egsInformation("EGS_PhspSource: using phsp portion between %lld and %lld\n",
352  Nfirst,Nlast);
353 }
354 
355 void EGS_PhspSource::readParticle() {
356  if ((++Npos) > Nlast) {
357  egsWarning("EGS_PhspSource::readParticle(): reached the end of the "
358  "phase space file chunk (%lld)\n will start from the beginning "
359  "of the chunk (%lld) but this "
360  "implies that uncertainty estimates will be inaccurate\n",
361  Nlast,Nfirst);
362  the_file.seekg(recl,ios::beg);
363  Nrestart++;
364  Npos = Nfirst;
365  }
366  the_file.read(record,recl*sizeof(char));
367  ++Nread;
368  if (the_file.eof() || !the_file.good())
369  egsFatal("EGS_PhspSource::readParticle(): I/O error while reading "
370  "phase space file\n");
371  __egs_data32 *data;
372  data = (__egs_data32 *) &record[0];
373  p.latch = data->i;
374  data = (__egs_data32 *) &record[4];
375  p.E = data->f;
376  data = (__egs_data32 *) &record[8];
377  p.x = data->f;
378  data = (__egs_data32 *) &record[12];
379  p.y = data->f;
380  data = (__egs_data32 *) &record[16];
381  p.u = data->f;
382  data = (__egs_data32 *) &record[20];
383  p.v = data->f;
384  data = (__egs_data32 *) &record[24];
385  p.wt = data->f;
386  if (swap_bytes) {
387  egsSwapBytes(&p.latch);
388  egsSwapBytes(&p.E);
389  egsSwapBytes(&p.wt);
390  egsSwapBytes(&p.x);
391  egsSwapBytes(&p.y);
392  egsSwapBytes(&p.u);
393  egsSwapBytes(&p.v);
394  }
395  if (p.latch & 1073741824) {
396  p.q = -1;
397  }
398  else if (p.latch & 536870912) {
399  p.q = 1;
400  }
401  else {
402  p.q = 0;
403  }
404  if (p.E < 0) {
405  count++;
406  p.E = -p.E;
407  }
408  else if (first) {
409  ++count;
410  }
411  first = false;
412  if (p.q) {
413  p.E -= EGS_Application::activeApplication()->getRM();
415  }
416  else {
418  }
419  Nuse = 0;
420  p.wt /= (Nrecycle+1);
421 }
422 
423 bool EGS_PhspSource::rejectParticle() const {
424  if (particle_type < 2 && p.q != particle_type) {
425  return true;
426  }
427  if (particle_type == 3 && !p.q) {
428  return true;
429  }
430  if (p.x < Xmin || p.x > Xmax || p.y < Ymin || p.y > Ymax) {
431  return true;
432  }
433  if (p.wt < wmin || p.wt > wmax) {
434  return true;
435  }
436  if (p.latch & 2147483648UL) {
437  return true;
438  }
439  if (filter_type < 0) {
440  return false;
441  }
442  if (filter_type == 0) {
443  bool r1;
444  if (filter1) {
445  r1 = !(p.latch & filter1);
446  }
447  else {
448  r1 = false;
449  }
450  bool r2 = (p.latch & filter2);
451  return (r1 || r2);
452  }
453  if (filter_type == 1) {
454  return (p.latch & filter1);
455  }
456  int l = p.latch >> 24;
457  bool res = l & filter1;
458  return filter_type == 3 ? res : !res;
459 }
460 
461 void EGS_PhspSource::setFilter(int type, int nbit1, int nbit2, const int *bits) {
462  if (type < 0 || type > 3) {
463  egsWarning("EGS_PhspSource::setFilter: invalid filter type %d\n",type);
464  return;
465  }
466  int ntot = nbit1;
467  if (type == 0) {
468  ntot += nbit2;
469  }
470  if (ntot > 29) {
471  egsWarning("EGS_PhspSource::setFilter: maximum number of bits is "
472  "limited to 29, you requested %d\n",ntot);
473  return;
474  }
475  if (!ntot) {
476  return;
477  }
478  filter_type = type;
479  if (filter_type == 0) {
480  int i;
481  filter1 = 0;
482  filter2 = 0;
483  for (i=0; i<ntot; i++) {
484  if (bits[i]) {
485  unsigned int aux = 1;
486  aux = (aux << i);
487  if (i < nbit1) {
488  filter1 += aux;
489  }
490  else {
491  filter2 += aux;
492  }
493  }
494  }
495  }
496  else if (filter_type == 1) { // TODO
497  }
498  else if (filter_type == 2) { // TODO
499  }
500  else if (filter_type == 3) { // TODO
501  }
502 }
503 
504 extern "C" {
505 
506  EGS_PHSP_SOURCE_EXPORT EGS_BaseSource *createSource(EGS_Input *input,
507  EGS_ObjectFactory *f) {
508  return
509  createSourceTemplate<EGS_PhspSource>(input,f,"phsp source");
510  }
511 
512 }
ifstream the_file
Phase space data stream.
EGS_Float Emin
Minimum energy (obtained from the phsp file)
EGS_Float Pinc
Number of incident particles that created the file.
int Nrestart
Number of times the file was restarted.
char * record
Memory to read a particle into.
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
Global egspp functions header file.
EGS_I64 Nlast
Last record this source can use.
bool mode2
true, if a MODE2 file
const EGS_Float veryFar
A very large float.
void egsSwapBytes(int *n)
Swap the bytes of 32 bit integers.
EGS_I64 Nread
Number of particles read 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_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
EGS_I64 Nparticle
Number of particles in the file.
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.
static EGS_Application * activeApplication()
Get the active application.
EGS_I64 Npos
Next record to be read.
string the_file_name
The phase-space file name.
int Nrecycle
Number of times to recycle current particle.
EGS_Float x
x-component
Definition: egs_vector.h:60
int Nrecycle_g
Number of times to recycle a photon.
EGS_PhspSource(const string &phsp_file, const string &Name="", EGS_ObjectFactory *f=0)
Constructor.
EGS_I64 Nfirst
first record this source can use
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_e
Number of times to recycle a charged particle.
string otype
The object type.
int recl
The particle record length.
EGS_I64 Nphoton
Number of photons in the file.
EGS_Float Emax
Maximum energy (obtained from the phsp file)
EGS_Application class header file.
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
int Nuse
Number of times current particle was used so far.
Base source class. All particle sources must be derived from this class.
A phase-space file source.
string description
A short source description.
EGS_InfoFunction EGS_EXPORT egsWarning
Always use this function for reporting warnings.