45 #include "iaea_phsp.h"
47 EGS_PhspScoring::EGS_PhspScoring(
const string &Name,
50 count(0), countg(0), emin(1.e30), emax(-1.e30), first_flush(true), is_restart(false) {
51 otype =
"EGS_PhspScoring";
54 EGS_PhspScoring::~EGS_PhspScoring() {
66 p_stack =
new Particle[store_max];
68 description =
"\n*******************************************\n";
72 description +=
"*******************************************\n";
74 description +=
"\n Will output phase space for particles crossing surfaces of geometry:\n";
78 description +=
"entering and exiting phase space geometry";
80 else if (scoredir == 1) {
83 else if (scoredir == 2) {
87 else if (score_type==1) {
88 description +=
"\n Will output phase space for the following exit/entry region pairs:\n";
90 for (
int i=0; i<fromreg.size(); i++) {
91 sprintf(buf,
"%d ",fromreg[i]);
95 for (
int i=0; i<fromreg.size(); i++) {
96 sprintf(buf,
"%d ",toreg[i]);
101 int nreg =
app->getnRegions();
102 for (
int ir=0; ir<nreg; ir++) {
104 for (
int i=0; i<fromreg.size(); i++) {
105 if (fromreg[i]==ir) {
106 tolist.push_back(toreg[i]);
109 from_to.push_back(tolist);
114 if (phspoutdir ==
"") {
118 description +=
"\n Data will be output in EGSnrc format.\n";
129 else if (oformat==1) {
130 description +=
"\n Data will be output in IAEA format.\n";
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';
157 iaea_n_extra_float=1;
162 iaea_n_extra_float=0;
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]) {
173 description +=
"\n Data scored at constant " + xyzname[i] +
" = " + xyz.str() +
" cm";
181 else if (ocharge == 1) {
184 else if (ocharge == 2) {
187 if (oformat ==1 && score_mu) {
188 description +=
"\n mu will be scored (if available)";
190 if (oformat ==0 && score_mc) {
191 description +=
"\n will score multiple crossers (and descendents)";
198 void EGS_PhspScoring::reportResults() {
202 iaea_destroy_source(&iaea_id,&iaea_iostat);
204 egsFatal(
"\n EGS_PhspScoring: Error closing phase space file.\n");
207 else if (oformat == 0) {
210 egsInformation(
"\n======================================================\n");
212 egsInformation(
"======================================================\n");
215 egsInformation(
" Header file: %s.IAEAheader\n",phsp_fname.c_str());
218 else if (oformat == 0) {
223 if (count == countg) {
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");
242 void EGS_PhspScoring::storeParticle(EGS_I64 ncase) {
246 egsWarning(
"\nEGS_PhspScoring: User requested mu scoring, but mu is inavailable with this source.\n");
252 EGS_Float prm =
app->getRM();
267 if (ncase != last_case) {
271 p_stack[phsp_index].E = E;
281 p_stack[phsp_index].mu =
app->
getMU();
285 if (phsp_index > store_max - 1) {
292 void EGS_PhspScoring::openPhspFile()
const {
296 phsp_file.open(phsp_fname.c_str(),ios::binary|ios::out|ios::in);
298 egsFatal(
"\nEGS_PhspScoring: Failed to open phase space file %s for appending.\n",
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());
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());
314 phsp_file.open(phsp_fname.c_str(),ios::binary|ios::out);
316 egsFatal(
"\nEGS_PhspScoring: Failed to open phase space file %s for writing.\n",
319 phsp_file.write(
"MODE0",5);
321 phsp_file.seekp(28,ios::beg);
324 else if (oformat == 1) {
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());
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());
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());
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());
357 for (
int i=0; i<3; i++) {
358 if (xyz_is_const[i]) {
360 float constval = xyzscore[i];
361 iaea_set_constant_variable(&iaea_id,&index,&constval);
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;
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);
376 iaea_set_type_extrafloat_variable(&iaea_id,&iaea_i_mu_tmp,&mu_ind_tmp);
383 int EGS_PhspScoring::flushBuffer()
const {
393 for (
int j=0; j<phsp_index; j++) {
396 int n_stat = p_stack[j].E < 0 ? 1 : 0;
397 float E = abs(p_stack[j].E);
399 int type = iaea_q_type[p_stack[j].q+1];
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];
406 iaea_extra_float[iaea_i_mu] = p_stack[j].mu;
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;
419 iaea_write_particle(&iaea_id,&n_stat,&type,&E,&wt,&x,&y,&z,&u,&v,&w,iaea_extra_float,iaea_extra_long);
421 egsFatal(
"\nEGS_PhspScoring: Failed to write particle data to phase space file.");
425 EGS_I64 last_case_tmp = last_case;
426 iaea_set_total_original_particles(&iaea_id,&last_case_tmp);
429 iaea_update_header(&iaea_id,&iaea_iostat);
431 egsFatal(
"\nEGS_PhspScoring: Failed to update phase space header file.");
434 else if (oformat == 0) {
436 egsFatal(
"\nEGS_PhspScoring: phase space file is not open for writing.");
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));
444 iostream::off_type pos = (count+1)*
sizeof(egs_phsp_write_struct);
446 phsp_file.seekp(5,ios::beg);
447 unsigned int count4 = count, countg4 = countg;
448 float pinc = last_case;
450 if (countg4==count4) {
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);
469 bool EGS_PhspScoring::storeState(ostream &data)
const {
481 bool EGS_PhspScoring::setState(istream &data) {
493 bool EGS_PhspScoring::addState(istream &data) {
515 const static char *func =
"createAusgabObject(phsp_scoring)";
523 vector <int> from_reg, to_reg;
530 bool xyzisconst[3] = {
false,
false,
false};
533 int err01 = input->
getInput(
"phase space geometry",gname);
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());
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");
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);
558 egsFatal(
"\nEGS_PhspScoring: Invalid scoring direction.\n");
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");
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");
579 from_reg.erase(p,p+from_reg.size()-to_reg.size());
581 else if (to_reg.size() > from_reg.size()) {
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());
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]);
594 p=from_reg.begin()+i;
595 p1 =to_reg.begin()+i;
607 if (input->
getInput(
"output format", str) < 0) {
608 egsInformation(
"EGS_PhspScoring: No input for output format type. Will default to EGSnrc.\n");
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");
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]);
625 xyzisconst[0] =
true;
628 xyzisconst[1] =
true;
631 xyzisconst[2] =
true;
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);
641 egsWarning(
"\nEGS_PhspScoring: Invalid input for mu scoring. Will not score mu.\n");
647 if (phspouttype == 0) {
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);
655 egsWarning(
"\nEGS_PhspScoring: Invalid input for score multiple crossers. Will not score.\n");
660 if (input->
getInput(
"output directory",outdir) < 0) {
663 if (input->
getInput(
"particle type", str) < 0) {
664 egsInformation(
"EGS_PhspScoring: No input for particle type. Will score all.\n");
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);
675 egsFatal(
"\nEGS_PhspScoring: Invalid particle type.\n");
685 result->setGeom(phspgeom);
688 result->setEntryExitReg(from_reg,to_reg);
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);
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.
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.
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
int getNparallel() const
Returns the number of parallel jobs executing.
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...
void setName(EGS_Input *inp)
Set the name of the object from the information provided by inp.
int latch
latch variable (useful as a flag on many occasions)
string name
The object name.
EGS_Float getMU()
Returns the value of the mu synchronization parameter.
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.
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.