38 #include "iaea_phsp.h"
49 void IAEA_PhspSource::init() {
50 otype =
"IAEA_PhspSource";
76 void IAEA_PhspSource::openFile(
const string &phsp_file) {
85 int len = phsp_file.length();
86 char phsp_name_tmp[len+1];
87 phsp_file.copy(phsp_name_tmp,len,0);
89 phsp_name_tmp[len]=
'\0';
94 egsWarning(
"IAEA_PhspSource::openFile: no header file name supplied.\n");
98 egsWarning(
"IAEA_PhspSource::openFile: failed to open header file %s.IAEAheader\n"
99 " for reading\n",phsp_file.c_str());
103 egsWarning(
"IAEA_PhspSource::openFile: failed to open phase space file %s.IAEAphsp\n"
104 " for reading\n",phsp_file.c_str());
108 egsWarning(
"IAEA_PhspSource::openFile: failed to initialize phase space source %s\n",phsp_file.c_str());
112 egsWarning(
"IAEA_PhspSource::openFile: failed to get record contents from header file %s.IAEAheader\n",phsp_file.c_str());
116 egsWarning(
"IAEA_PhspSource::openFile: I/O error ocurred on opening phase space file %s\n",phsp_file.c_str());
125 egsWarning(
"IAEA_PhspSource::openFile: error reading file size/byte order from header %s.IAEAheader\n",phsp_file.c_str());
129 egsWarning(
"IAEA_PhspSource::openFile: byte mismatch between phase space file and machine. Will swap bytes.\n");
133 egsWarning(
"IAEA_PhspSource::openFile: mismatch between file size in header and actual file size of %s\n",phsp_file.c_str());
138 EGS_I64 n, n_photon, pinc;
141 iaea_get_max_particles(&
iaea_fileid,&iaea_type,&n);
143 egsWarning(
"IAEA_PhspSource::openFile: failed to get total no. of particles from %s.IAEAheader\n",phsp_file.c_str());
147 iaea_get_max_particles(&
iaea_fileid,&iaea_type,&n_photon);
149 egsWarning(
"IAEA_PhspSource::openFile: failed to get no. of photons from %s.IAEAheader\n",phsp_file.c_str());
155 egsWarning(
"IAEA_PhspSource::openFile: failed to get max. energy from %s.IAEAheader\n",phsp_file.c_str());
159 iaea_get_total_original_particles(&
iaea_fileid,&pinc);
161 egsWarning(
"IAEA_PhspSource::openFile: failed to get no. of original particles from %s.IAEAheader\n",phsp_file.c_str());
168 egsWarning(
"IAEA_PhspSource::openFile: failed to get no. of extra floats and longs stored in %s.IAEAphsp\n",phsp_file.c_str());
172 int extrafloat_types[MAXEXTRAS], extralong_types[MAXEXTRAS];
175 egsWarning(
"IAEA_PhspSource::openFile: failed to get Mode of data %s.IAEAheader\n",phsp_file.c_str());
181 if (extrafloat_types[i]==3) {
192 if (n_extra_floats >
i_zlast+1) {
193 if (extrafloat_types[
i_zlast+1] == 0) {
196 egsInformation(
"IAEA_PhspSource::openFile: Mu index included in data in %s.IAEAphsp\n",phsp_file.c_str());
203 if (extralong_types[i]==2) {
210 egsWarning(
"IAEA_PhspSource::openFile: LATCH is not stored in data in %s.IAEAphsp\n",phsp_file.c_str());
229 int err = input->
getInput(
"iaea phase space file",fname);
231 egsWarning(
"IAEA_PhspSource: no 'iaea phase space file' input\n");
236 egsWarning(
"IAEA_PhspSource: errors while opening the phase space file"
237 " %s\n",fname.c_str());
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]);
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);
257 if (the_filter[0] >= 0 && the_filter[0] <= 3) {
258 int nbit1 = the_latch[0].size();
260 if (!the_filter[0]) {
261 nbit2 = the_latch[1].size();
263 if (nbit1 + nbit2 > 0) {
264 int *the_bits =
new int [nbit1+nbit2];
265 for (
int j=0; j<nbit1+nbit2; j++) {
267 the_bits[j]=(int)(the_latch[0][j] ==
'1');
270 the_bits[j] = (int)(the_latch[1][j-nbit1] ==
'1');
273 setFilter(the_filter[0],nbit1,nbit2,the_bits);
278 vector<EGS_Float> wwindow;
279 err = input->
getInput(
"weight window",wwindow);
280 if (!err && wwindow.size() == 2) {
285 err = input->
getInput(
"reuse photons",ntmp);
286 if (!err && ntmp > 0) {
290 err = input->
getInput(
"recycle photons",ntmp);
291 if (!err && ntmp > 0) {
295 err = input->
getInput(
"reuse electrons",ntmp);
296 if (!err && ntmp > 0) {
300 err = input->
getInput(
"recycle electrons",ntmp);
301 if (!err && ntmp > 0) {
316 int nstat,extrainttemp[MAXEXTRAS];
317 float extrafloattemp[MAXEXTRAS];
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",
327 egsFatal(
"IAEA_PhspSource::getNextParticle(): error restarting phase space chunk\n");
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);
341 p.zlast = extrafloattemp[
i_zlast];
344 p.mu = extrafloattemp[
i_mu];
363 egsFatal(
"IAEA_PhspSource::getNextParticle(): error reading particle number %i\n",
Npos);
376 egsFatal(
"IAEA_PhspSource::getNextParticle: unknown charge on particle\n");
379 if (first || nstat>0) {
400 u.
y = p.v, u.
z = p.w;
402 if (rejectParticle()) {
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);
434 Nlast = nstart + nrun;
438 egsWarning(
"IAEA_PhspSource::setSimulationChunk(): error setting phase space chunk\n");
441 egsInformation(
"IAEA_PhspSource: using phsp portion between %lld and %lld\n",
445 bool IAEA_PhspSource::rejectParticle()
const {
446 if (particle_type < 2 && p.q != particle_type) {
449 if (particle_type == 3 && !p.q) {
452 if (p.x < Xmin || p.x > Xmax || p.y < Ymin || p.y > Ymax) {
455 if (p.wt < wmin || p.wt > wmax) {
458 if (p.latch & 2147483648UL) {
461 if (filter_type < 0) {
464 if (filter_type == 1) {
467 r1 = !(p.latch & filter1);
472 bool r2 = (p.latch & filter2);
475 if (filter_type == 1) {
476 return (p.latch & filter1);
478 int l = p.latch >> 24;
479 bool res = l & filter1;
480 return filter_type == 3 ? res : !res;
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);
493 egsWarning(
"IAEA_PhspSource::setFilter: maximum number of bits is "
494 "limited to 29, you requested %d\n",ntot);
501 if (filter_type == 0) {
505 for (i=0; i<ntot; i++) {
507 unsigned int aux = 1;
518 else if (filter_type == 1) {
520 else if (filter_type == 2) {
522 else if (filter_type == 3) {
531 createSourceTemplate<IAEA_PhspSource>(input,f,
"iaea phsp source");
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.
A class representing 3D vectors.
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...
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...
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.
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.
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.
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