50 void EGS_PhspSource::init() {
51 otype =
"EGS_PhspSource";
79 void EGS_PhspSource::openFile(
const string &phsp_file) {
90 the_file.open(phsp_file.c_str(),ios::binary | ios::in);
92 egsWarning(
"EGS_PhspSource::openFile: failed to open binary file %s"
93 " for reading\n",phsp_file.c_str());
98 for (
int i=0; i<5; i++) {
101 egsWarning(
"EGS_PhspSource::openFile: an I/O error occured "
102 "while reading the first record of %s\n",phsp_file.c_str());
107 if (cmode ==
"MODE0") {
111 else if (cmode ==
"MODE2") {
116 egsWarning(
"EGS_PhspSource::openFile: the file %s is not a MODE0 or"
117 " MODE2 file\n",phsp_file.c_str());
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));
129 egsWarning(
"EGS_PhspSource::openFile: an I/O error occured "
130 "while reading the first record of %s\n",phsp_file.c_str());
135 if (n <= 0 || n_photon < 0 || n_photon > n || emin < 0 || emax < 0 ||
136 emax < emin || pinc < 0) {
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");
151 egsWarning(
" number of photons: %d\n",n_photon);
153 if (n_photon > n)
egsWarning(
" number of photons (%d) is "
154 "greater than number of particles (%d)\n",n_photon,n);
162 egsWarning(
" emin > emax: %g %g\n",emin,emax);
175 istream::off_type nend = n;
181 egsWarning(
"EGS_PhspSource::openFile: failed to read the last"
182 " particle in the file, this indicates some unknown error condition\n");
205 int err = input->
getInput(
"phase space file",fname);
207 egsWarning(
"EGS_PhspSource: no 'phase space file' input\n");
212 egsWarning(
"EGS_PhspSource: errors while opening the phase space file"
213 " %s\n",fname.c_str());
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]);
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);
233 if (the_filter[0] >= 0 && the_filter[0] <= 3) {
234 int nbit1 = the_latch[0].size();
236 if (!the_filter[0]) {
237 nbit2 = the_latch[1].size();
239 if (nbit1 + nbit2 > 0) {
240 int *the_bits =
new int [nbit1+nbit2];
241 for (
int j=0; j<nbit1+nbit2; j++) {
243 the_bits[j]=(int)(the_latch[0][j] ==
'1');
246 the_bits[j] = (int)(the_latch[1][j-nbit1] ==
'1');
249 setFilter(the_filter[0],nbit1,nbit2,the_bits);
254 vector<EGS_Float> wwindow;
255 err = input->
getInput(
"weight window",wwindow);
256 if (!err && wwindow.size() == 2) {
261 err = input->
getInput(
"reuse photons",ntmp);
262 if (!err && ntmp > 0) {
266 err = input->
getInput(
"recycle photons",ntmp);
267 if (!err && ntmp > 0) {
271 err = input->
getInput(
"reuse electrons",ntmp);
272 if (!err && ntmp > 0) {
276 err = input->
getInput(
"recycle electrons",ntmp);
277 if (!err && ntmp > 0) {
287 if (!recl)
egsFatal(
"EGS_PhspSource::readParticle(): the file is not "
302 EGS_Float aux = p.u*p.u+p.v*p.v;
317 if (rejectParticle()) {
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);
347 Nlast = nstart + nrun;
351 egsInformation(
"EGS_PhspSource: using phsp portion between %lld and %lld\n",
355 void EGS_PhspSource::readParticle() {
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",
369 egsFatal(
"EGS_PhspSource::readParticle(): I/O error while reading "
370 "phase space file\n");
372 data = (__egs_data32 *) &
record[0];
374 data = (__egs_data32 *) &
record[4];
376 data = (__egs_data32 *) &
record[8];
378 data = (__egs_data32 *) &
record[12];
380 data = (__egs_data32 *) &
record[16];
382 data = (__egs_data32 *) &
record[20];
384 data = (__egs_data32 *) &
record[24];
395 if (p.latch & 1073741824) {
398 else if (p.latch & 536870912) {
423 bool EGS_PhspSource::rejectParticle()
const {
424 if (particle_type < 2 && p.q != particle_type) {
427 if (particle_type == 3 && !p.q) {
430 if (p.x < Xmin || p.x > Xmax || p.y < Ymin || p.y > Ymax) {
433 if (p.wt < wmin || p.wt > wmax) {
436 if (p.latch & 2147483648UL) {
439 if (filter_type < 0) {
442 if (filter_type == 0) {
445 r1 = !(p.latch & filter1);
450 bool r2 = (p.latch & filter2);
453 if (filter_type == 1) {
454 return (p.latch & filter1);
456 int l = p.latch >> 24;
457 bool res = l & filter1;
458 return filter_type == 3 ? res : !res;
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);
471 egsWarning(
"EGS_PhspSource::setFilter: maximum number of bits is "
472 "limited to 29, you requested %d\n",ntot);
479 if (filter_type == 0) {
483 for (i=0; i<ntot; i++) {
485 unsigned int aux = 1;
496 else if (filter_type == 1) {
498 else if (filter_type == 2) {
500 else if (filter_type == 3) {
509 createSourceTemplate<EGS_PhspSource>(input,f,
"phsp source");
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.
A class representing 3D vectors.
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...
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...
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.
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
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 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.