diff --git a/HEN_HOUSE/egs++/Makefile b/HEN_HOUSE/egs++/Makefile
index 0e6e55e5b..1be6e3872 100644
--- a/HEN_HOUSE/egs++/Makefile
+++ b/HEN_HOUSE/egs++/Makefile
@@ -62,7 +62,7 @@ shape_libs = egs_circle egs_ellipse egs_extended_shape egs_gaussian_shape \
              egs_line_shape egs_polygon_shape egs_rectangle egs_shape_collection \
              egs_voxelized_shape
 
-aobject_libs = egs_track_scoring egs_dose_scoring egs_radiative_splitting
+aobject_libs = egs_track_scoring egs_dose_scoring egs_radiative_splitting egs_phsp_scoring
 
 all_libs = $(geometry_libs) $(source_libs) $(shape_libs) $(aobject_libs)
 lib_objects = $(addprefix $(DSO1), $(addsuffix .$(obje), $(all_libs)))
diff --git a/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/Makefile b/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/Makefile
new file mode 100644
index 000000000..04eedc2e3
--- /dev/null
+++ b/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/Makefile
@@ -0,0 +1,53 @@
+
+###############################################################################
+#
+#  EGSnrc egs++ makefile to build phase space scoring object
+#  Copyright (C) 2015 National Research Council Canada
+#
+#  This file is part of EGSnrc.
+#
+#  EGSnrc is free software: you can redistribute it and/or modify it under
+#  the terms of the GNU Affero General Public License as published by the
+#  Free Software Foundation, either version 3 of the License, or (at your
+#  option) any later version.
+#
+#  EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
+#  WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+#  FOR A PARTICULAR PURPOSE.  See the GNU Affero General Public License for
+#  more details.
+#
+#  You should have received a copy of the GNU Affero General Public License
+#  along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
+#
+###############################################################################
+#
+#  Author:          Blake Walters, 2018
+#
+#  Contributors:
+#
+###############################################################################
+
+
+include $(EGS_CONFIG)
+include $(SPEC_DIR)egspp.spec
+include $(SPEC_DIR)egspp_$(my_machine).conf
+
+DEFS = $(DEF1) -DBUILD_PHSP_SCORING_DLL
+
+library = egs_phsp_scoring
+lib_files = egs_phsp_scoring
+my_deps = $(common_ausgab_deps)
+extra_dep = $(addprefix $(DSOLIBS), $(my_deps))
+
+include $(SPEC_DIR)egspp_libs.spec
+
+INC2 = -I$(IEGS2) -I..$(DSEP).. -I$(HEN_HOUSE)iaea_phsp
+
+#this has been redefined to include iaea_phsp libraries
+link2_libs = egspp iaea_phsp
+
+$(make_depend)
+
+test:
+	@echo "common_h2: $(common_h2)"
+	@echo "extra_dep: $(extra_dep)"
diff --git a/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/egs_phsp_scoring.cpp b/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/egs_phsp_scoring.cpp
new file mode 100644
index 000000000..79059f4ee
--- /dev/null
+++ b/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/egs_phsp_scoring.cpp
@@ -0,0 +1,700 @@
+/*
+###############################################################################
+#
+#  EGSnrc egs++ phase space scoring object
+#  Copyright (C) 2015 National Research Council Canada
+#
+#  This file is part of EGSnrc.
+#
+#  EGSnrc is free software: you can redistribute it and/or modify it under
+#  the terms of the GNU Affero General Public License as published by the
+#  Free Software Foundation, either version 3 of the License, or (at your
+#  option) any later version.
+#
+#  EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
+#  WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+#  FOR A PARTICULAR PURPOSE.  See the GNU Affero General Public License for
+#  more details.
+#
+#  You should have received a copy of the GNU Affero General Public License
+#  along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
+#
+###############################################################################
+#
+#  Author:          Blake Walters, 2018
+#
+#  Contributors:
+#
+###############################################################################
+*/
+
+
+/*! \file egs_phsp_scoring.cpp
+ *  \brief A phase space scoring ausgab object: implementation
+ *  \BW
+ */
+
+#include <sstream>
+#include <fstream>
+#include <string>
+#include <cstdlib>
+
+#include "egs_phsp_scoring.h"
+#include "egs_input.h"
+#include "egs_functions.h"
+#include "iaea_phsp.h"
+
+EGS_PhspScoring::EGS_PhspScoring(const string &Name,
+                                 EGS_ObjectFactory *f) :
+    EGS_AusgabObject(Name,f), phsp_index(0), store_max(1000), phsp_file(0),
+    count(0), countg(0), emin(1.e30), emax(-1.e30), first_flush(true), is_restart(false) {
+    otype = "EGS_PhspScoring";
+}
+
+EGS_PhspScoring::~EGS_PhspScoring() {
+}
+
+void EGS_PhspScoring::setApplication(EGS_Application *App) {
+    EGS_AusgabObject::setApplication(App);
+    if (!app) {
+        return;
+    }
+
+    char buf[512];//useful character buffer
+    //set up the stack of particles to output to the phase space file
+    //1000 particles at a time
+    p_stack = new Particle[store_max];
+
+    description = "\n*******************************************\n";
+    description +=  "Phase Space Scoring Object (";
+    description += name;
+    description += ")\n";
+    description += "*******************************************\n";
+    if (score_type==0) {
+        description += "\n Will output phase space for particles crossing surfaces of geometry:\n";
+        description += phsp_geom->getName();
+        description += "\n Particles scored on: ";
+        if (scoredir == 0) {
+            description += "entering and exiting phase space geometry";
+        }
+        else if (scoredir == 1) {
+            description += "entering phase space geometry";
+        }
+        else if (scoredir == 2) {
+            description += "exiting phase space geometry";
+        }
+    }
+    else if (score_type==1) {
+        description += "\n Will output phase space for the following exit/entry region pairs:\n";
+        description += " Exit regions:\n";
+        for (int i=0; i<fromreg.size(); i++) {
+            sprintf(buf,"%d ",fromreg[i]);
+            description += buf;
+        }
+        description += "\n Entry regions:\n";
+        for (int i=0; i<fromreg.size(); i++) {
+            sprintf(buf,"%d ",toreg[i]);
+            description += buf;
+        }
+        description += "\n";
+        //now, for every region, set up a (possibly empty) array of entry regions
+        int nreg = app->getnRegions();
+        for (int ir=0; ir<nreg; ir++) {
+            vector<int> tolist;
+            for (int i=0; i<fromreg.size(); i++) {
+                if (fromreg[i]==ir) {
+                    tolist.push_back(toreg[i]);
+                }
+            }
+            from_to.push_back(tolist);
+        }
+    }
+
+    //construct name of phase space file -- opening to occur later
+    if (phspoutdir == "") {
+        phspoutdir = app->getAppDir();
+    }
+    if (oformat==0) {
+        description += "\n Data will be output in EGSnrc format.\n";
+        if (app->getNparallel()>0) {
+            sprintf(buf,"%s_w%d.egsphsp1",getObjectName().c_str(),app->getIparallel());
+        }
+        else {
+            sprintf(buf,"%s.egsphsp1",getObjectName().c_str());
+        }
+        phsp_fname=egsJoinPath(phspoutdir,buf);
+        description += "\n Phase space file name:\n";
+        description += phsp_fname;
+    }
+    else if (oformat==1) {
+        description += "\n Data will be output in IAEA format.\n";
+        if (app->getNparallel()>0) {
+            sprintf(buf,"%s_w%d.1",getObjectName().c_str(),app->getIparallel());
+        }
+        else {
+            sprintf(buf,"%s.1",getObjectName().c_str());
+        }
+        phsp_fname=egsJoinPath(phspoutdir,buf);
+
+        //need to add a null terminator
+        len = phsp_fname.length();
+        phsp_fname_char = new char[len+1];
+        phsp_fname.copy(phsp_fname_char,len,0);
+        phsp_fname_char[len]='\0'; //null terminator on string
+
+        //iaea particle charge
+        iaea_q_type[0]=2; //e-
+        iaea_q_type[1]=1; //photon
+        iaea_q_type[2]=3; //e+
+
+        iaea_id = 1;
+
+        //set up extra float and extra long array indices
+        latch_ind = 2;//type of latch as defined by iaea
+        iaea_n_extra_long=1; //only store latch
+        iaea_i_latch=0; // position of latch in array
+        if (score_mu) {
+            iaea_n_extra_float=1;
+            mu_ind = 0; //set type to generic float
+            iaea_i_mu=0; //position of mu in array
+        }
+        else {
+            iaea_n_extra_float=0; //no extra floats
+        }
+
+        description += "\n Phase space file names:\n";
+        description += "  Header file: " + phsp_fname + ".IAEAheader\n";
+        description += "  Data file: " + phsp_fname + ".IAEAphsp";
+        string xyzname[3] = {"X", "Y", "Z"};
+        for (int i=0; i<3; i++) {
+            if (xyz_is_const[i]) {
+                ostringstream xyz;
+                xyz << xyzscore[i];
+                description += "\n Data scored at constant " + xyzname[i] + " = " + xyz.str() + " cm";
+            }
+        }
+    }
+    description += "\n Particles scored: ";
+    if (ocharge == 0) {
+        description += "all";
+    }
+    else if (ocharge == 1) {
+        description += "photons";
+    }
+    else if (ocharge == 2) {
+        description += "charged";
+    }
+    if (oformat ==1 && score_mu) {
+        description += "\n mu will be scored (if available)";
+    }
+    if (oformat ==0 && score_mc) {
+        description += "\n will score multiple crossers (and descendents)";
+    }
+}
+
+//final buffer flush and then close file
+//we want to output some data about particles scored
+//may ultimately want to allow the user to define scoring zones for output
+void EGS_PhspScoring::reportResults() {
+    flushBuffer();
+    if (oformat == 1) { //iaea
+        int iaea_iostat;
+        iaea_destroy_source(&iaea_id,&iaea_iostat);
+        if (iaea_iostat<0) {
+            egsFatal("\n EGS_PhspScoring: Error closing phase space file.\n");
+        }
+    }
+    else if (oformat == 0) {
+        phsp_file.close();
+    }
+    egsInformation("\n======================================================\n");
+    egsInformation("Phase Space Scoring Object(%s)\n",name.c_str());
+    egsInformation("======================================================\n");
+    if (oformat==1) {
+        egsInformation("\n IAEA format phase space output:\n");
+        egsInformation(" Header file: %s.IAEAheader\n",phsp_fname.c_str());
+        egsInformation(" Data file: %s.IAEAphsp\n",phsp_fname.c_str());
+    }
+    else if (oformat == 0) {
+        egsInformation("\n EGSnrc format phase space output:\n");
+        egsInformation(" Data file: %s\n",phsp_fname.c_str());
+    }
+    float emintmp;
+    if (count == countg) {
+        emintmp = 0.0;
+    }
+    else {
+        emintmp = emin;
+    }
+    egsInformation("Summary of scored data:\n");
+    egsInformation("=> total no. of particles = %lld \n", count);
+    egsInformation("=> no. of photons = %lld \n", countg);
+    egsInformation("=> max. k.e. of all particles = %g MeV\n",emax);
+    egsInformation("=> min. k.e. of charged particles = %g MeV\n",emintmp);
+    egsInformation("=> no. of primary histories represented = %lld\n",last_case);
+    egsInformation("\n======================================================\n");
+}
+
+//store 1000 particles at a time in p_stack
+//if we're at 1000, actually write the particles to the file and update
+//the header info
+//also, keep track of phase space file counters, min., max. energy
+void EGS_PhspScoring::storeParticle(EGS_I64 ncase) {
+
+    //if user requested mu scoring, check if mu is available
+    if (score_mu && app->getMU() < 0) {
+        egsWarning("\nEGS_PhspScoring: User requested mu scoring, but mu is inavailable with this source.\n");
+        egsWarning("Turning off mu scoring.\n");
+        score_mu=false;
+    }
+
+    //counters, min. and max. k.e.
+    EGS_Float prm = app->getRM();
+    count++;
+    if (app->top_p.q==0) {
+        countg++;
+    }
+    if (app->top_p.E-abs(app->top_p.q)*prm > emax) {
+        emax = app->top_p.E-abs(app->top_p.q)*prm;
+    }
+    if (app->top_p.q != 0 && app->top_p.E - prm < emin) {
+        emin = app->top_p.E - prm;
+    }
+
+    //store particle data in p_stack
+    //set -ve energy marker if this is a new primary hist.
+    double E = app->top_p.E;
+    if (ncase != last_case) {
+        E = -E;
+        last_case = ncase;
+    }
+    p_stack[phsp_index].E = E;
+    p_stack[phsp_index].wt = app->top_p.wt;
+    p_stack[phsp_index].x = app->top_p.x.x;
+    p_stack[phsp_index].y = app->top_p.x.y;
+    p_stack[phsp_index].z = app->top_p.x.z;
+    p_stack[phsp_index].u = app->top_p.u.x;
+    p_stack[phsp_index].v = app->top_p.u.y;
+    p_stack[phsp_index].w = app->top_p.u.z;
+    p_stack[phsp_index].q = app->top_p.q;
+    if (score_mu) {
+        p_stack[phsp_index].mu = app->getMU();
+    }
+    p_stack[phsp_index++].latch = app->top_p.latch;
+
+    if (phsp_index > store_max - 1) {
+        //write store_max particles to the file and reset phsp_index counter
+        flushBuffer();
+    }
+}
+
+//open the phase space file for writing/appending data
+void EGS_PhspScoring::openPhspFile() const {
+//the file has already been named at this point
+    if (oformat==0) { //EGSnrc format
+        if (is_restart) {
+            phsp_file.open(phsp_fname.c_str(),ios::binary|ios::out|ios::in);
+            if (!(phsp_file))
+                egsFatal("\nEGS_PhspScoring: Failed to open phase space file %s for appending.\n",
+                         phsp_fname.c_str());
+            //check that total no. of particles in header = total no. read from .egsdat file
+            unsigned int count4;
+            phsp_file.seekg(5,ios::beg);
+            phsp_file.read((char *) &count4, sizeof(unsigned int));
+            if (count4 != countprev) {
+                egsFatal("\nEGS_PhspScoring: Particle no. mismatch between %s and .egsdat file.\n",phsp_fname.c_str());
+            }
+            //go to the end of the file
+            phsp_file.seekp(0,ios::end);
+            if (phsp_file.tellp() != 28 + countprev*sizeof(egs_phsp_write_struct)) {
+                egsFatal("\nEGS_PhspScoring: File size mismatch in %s and .egsdat file.\n",phsp_fname.c_str());
+            }
+        }
+        else {
+            phsp_file.open(phsp_fname.c_str(),ios::binary|ios::out);
+            if (!(phsp_file))
+                egsFatal("\nEGS_PhspScoring: Failed to open phase space file %s for writing.\n",
+                         phsp_fname.c_str());
+            //always MODE0 files--i.e. no ZLAST
+            phsp_file.write("MODE0",5);
+            //leave space for/skip over the rest of the header
+            phsp_file.seekp(28,ios::beg);
+        }
+    }
+    else if (oformat == 1) {  //IAEA format
+        int rwmode;
+        int iaea_iostat;
+        iaea_id = 1; //numerical index indicating this is the 1st file associated with this object scored
+        //hard coded to 1
+        if (is_restart) {
+            int rwmode = 3;
+            iaea_new_source(&iaea_id,phsp_fname_char,&rwmode,&iaea_iostat,len);
+            if (iaea_iostat < 0) {
+                egsFatal("\nEGS_PhspScoring: Failed to open phase space file %s.IAEAphsp for appending.\n",phsp_fname.c_str());
+            }
+            //check for consistency with total no. of scored particles as read from .egsdat file
+            EGS_I64 nparticle;
+            int type = -1;
+            int iaea_iostat;
+            iaea_get_max_particles(&iaea_id,&type,&nparticle);
+            if (nparticle != countprev) {
+                egsFatal("\nEGS_PhspScoring: Particle no. mismatch between %s.IAEAphsp and .egsdat file.\n",phsp_fname.c_str());
+            }
+            iaea_check_file_size_byte_order(&iaea_id,&iaea_iostat);
+            if (iaea_iostat != 0) {
+                egsFatal("\nEGS_PhspScoring: Byte order/file size mismatch in %s.IAEAphsp.\n",phsp_fname.c_str());
+            }
+        }
+        else {
+            int rwmode = 2;
+            iaea_new_source(&iaea_id,phsp_fname_char,&rwmode,&iaea_iostat,len);
+            if (iaea_iostat < 0) {
+                egsFatal("\nEGS_PhspScoring: Failed to open phase space file %s.IAEAphsp for writing.\n",phsp_fname.c_str());
+            }
+        }
+
+
+        //set up constant variables
+        for (int i=0; i<3; i++) {
+            if (xyz_is_const[i]) {
+                int index = i;
+                float constval = xyzscore[i];
+                iaea_set_constant_variable(&iaea_id,&index,&constval);
+            }
+        }
+        //set up extra floats and int indices and types
+        //need to store below in _tmp variables because this is a const function
+        int latch_ind_tmp = latch_ind;
+        int iaea_n_extra_long_tmp=iaea_n_extra_long;
+        int iaea_i_latch_tmp=iaea_i_latch;
+        int mu_ind_tmp = mu_ind;
+        int iaea_i_mu_tmp = iaea_i_mu;
+        int iaea_n_extra_float_tmp=iaea_n_extra_float;
+
+        iaea_set_extra_numbers(&iaea_id,&iaea_n_extra_float_tmp,&iaea_n_extra_long_tmp);
+        iaea_set_type_extralong_variable(&iaea_id,&iaea_i_latch_tmp,&latch_ind_tmp);
+        if (score_mu) {
+            iaea_set_type_extrafloat_variable(&iaea_id,&iaea_i_mu_tmp,&mu_ind_tmp);
+        }
+    }
+}
+
+//write phsp_index particles to phase space file
+//update header and reset phsp_index
+int EGS_PhspScoring::flushBuffer() const {
+
+    if (first_flush) {
+        openPhspFile();    //here's where we open the phase space file
+    }
+    //awkward logic because we do not know if this is a restart
+    //until after initialization
+    first_flush = false;
+
+    if (oformat == 1) { //iaea format
+        for (int j=0; j<phsp_index; j++) {
+            //fairly transparent, could probably put a lot of this in a separate method
+            //undo -ve energy marker and use n_stat to indicate new primary hist.
+            int n_stat = p_stack[j].E < 0 ? 1 : 0;
+            float E = abs(p_stack[j].E);
+            //convert charge to iaea type
+            int type = iaea_q_type[p_stack[j].q+1];
+
+            //store latch in iaea_extra_long
+            EGS_I32 *iaea_extra_long = new EGS_I32[iaea_n_extra_long];
+            iaea_extra_long[iaea_i_latch]=p_stack[j].latch;
+            float *iaea_extra_float = new float[iaea_n_extra_float];
+            if (score_mu) {
+                iaea_extra_float[iaea_i_mu] = p_stack[j].mu;
+            }
+
+            //now store double precision values in single precision reals
+            float wt = p_stack[j].wt;
+            float x = p_stack[j].x;
+            float y = p_stack[j].y;
+            float z = p_stack[j].z;
+            float u = p_stack[j].u;
+            float v = p_stack[j].v;
+            float w = p_stack[j].w;
+            //now actually write data
+
+            iaea_write_particle(&iaea_id,&n_stat,&type,&E,&wt,&x,&y,&z,&u,&v,&w,iaea_extra_float,iaea_extra_long);
+            if (n_stat < 0) {
+                egsFatal("\nEGS_PhspScoring: Failed to write particle data to phase space file.");
+            }
+        }
+        //update header with no. of primary histories
+        EGS_I64 last_case_tmp = last_case;
+        iaea_set_total_original_particles(&iaea_id,&last_case_tmp);
+        //update header file
+        int iaea_iostat;
+        iaea_update_header(&iaea_id,&iaea_iostat);
+        if (iaea_iostat<0) {
+            egsFatal("\nEGS_PhspScoring: Failed to update phase space header file.");
+        }
+    }
+    else if (oformat == 0) {  //EGSnrc format
+        if (!phsp_file) {
+            egsFatal("\nEGS_PhspScoring: phase space file is not open for writing.");
+        }
+        //don't forget that phsp_index is incremented after every particle written to p_stack
+        for (int j=0; j<phsp_index; j++) {
+            egs_phsp_write_struct ws(p_stack[j]);
+            phsp_file.write((char *) &ws, sizeof(ws));
+        }
+        //store position of end of file
+        iostream::off_type pos = (count+1)*sizeof(egs_phsp_write_struct);
+        //update header
+        phsp_file.seekp(5,ios::beg);
+        unsigned int count4 = count, countg4 = countg;
+        float pinc = last_case;
+        float emintmp;
+        if (countg4==count4) {
+            emintmp = 0.0;
+        }
+        else {
+            emintmp = emin;
+        }
+        phsp_file.write((char *) &count4, sizeof(unsigned int));
+        phsp_file.write((char *) &countg4, sizeof(unsigned int));
+        phsp_file.write((char *) &emax, sizeof(float));
+        phsp_file.write((char *) &emintmp, sizeof(float));
+        phsp_file.write((char *) &pinc, sizeof(float));
+        phsp_file.seekp(pos,ios::beg);
+    }
+
+    phsp_index=0;
+
+    return 0;
+};
+
+bool EGS_PhspScoring::storeState(ostream &data) const {
+    if (!egsStoreI64(data,count)) {
+        return false;
+    }
+    if (!egsStoreI64(data,countg)) {
+        return false;
+    }
+    //update phase space file at the end of each batch
+    flushBuffer();
+    return true;
+}
+
+bool  EGS_PhspScoring::setState(istream &data) {
+    if (!egsGetI64(data,count)) {
+        return false;
+    }
+    if (!egsGetI64(data,countg)) {
+        return false;
+    }
+    countprev = count;
+    is_restart = true;
+    return true;
+}
+
+bool  EGS_PhspScoring::addState(istream &data) {
+    EGS_I64 tmp;
+    if (!egsGetI64(data,tmp)) {  //return false;
+        egsWarning("error while reading count\n");
+        return false;
+    }
+    count += tmp;
+    if (!egsGetI64(data,tmp)) {  //return false;
+        egsWarning("error while reading countg\n");
+        return false;
+    }
+    countg += tmp;
+    return true;
+}
+
+//*********************************************************************
+// Process input for this ausgab object
+//**********************************************************************
+extern "C" {
+
+    EGS_PHSP_SCORING_EXPORT EGS_AusgabObject *createAusgabObject(EGS_Input *input,
+            EGS_ObjectFactory *f) {
+        const static char *func = "createAusgabObject(phsp_scoring)";
+        if (!input) {
+            egsWarning("%s: null input?\n",func);
+            return 0;
+        }
+        string str;
+        EGS_BaseGeometry *phspgeom;
+        int iscoremc = 0; //default to not score multiple crossers
+        vector <int> from_reg, to_reg;
+        int stype = 0; //default is to use scoring geom
+        int phspouttype;
+        int ptype;
+        int sdir;
+        int imuscore = 0;
+        float xyzconst[3];
+        bool xyzisconst[3] = {false, false, false};
+        string gname;
+        string outdir;
+        int err01 = input->getInput("phase space geometry",gname);
+        if (err01) {
+            stype = 1;
+        }
+        else {
+            phspgeom = EGS_BaseGeometry::getGeometry(gname);
+            if (!phspgeom) {
+                egsWarning("\nEGS_PhspScoring: %s does not name an existing geometry.\n"
+                           "Will assume you want to use exit/entry region pairs.\n",gname.c_str());
+                stype = 1;
+            }
+            else {
+                if (input->getInput("score particles on", str) < 0) {
+                    egsInformation("EGS_PhspScoring: No input for scoring direction.\n");
+                    egsInformation("Will score on entry and exit from phase space geometry.\n");
+                    sdir = 0;
+                }
+                else {
+                    //get scoring direction
+                    vector<string> allowed_sdir;
+                    allowed_sdir.push_back("entry and exit");
+                    allowed_sdir.push_back("entry");
+                    allowed_sdir.push_back("exit");
+                    sdir = input->getInput("score particles on",allowed_sdir,-1);
+                    if (sdir < 0) {
+                        egsFatal("\nEGS_PhspScoring: Invalid scoring direction.\n");
+                    }
+                }
+            }
+        }
+        if (stype==1) {
+            // user wants to use exit/entry region pairs
+            int err05 = input->getInput("from regions",from_reg);
+            int err06 = input->getInput("to regions",to_reg);
+            if (err05 || err06) {
+                egsFatal("\nEGS_PhspScoring: Missing/incorrect input for scoring method\n"
+                         "(scoring geometry or pairs of exit/entry regions)\n");
+            }
+            else {
+                //run some checks on exit/entry region pairs
+                vector<int>::iterator p,p1;
+                if (from_reg.size() > to_reg.size()) {
+                    p = from_reg.begin();
+                    egsWarning("\nEGS_PhspScoring: Mismatch in no. of exit/entry regions.\n"
+                               "Will only score for matched pairs.\n");
+                    p += to_reg.size();
+                    from_reg.erase(p,p+from_reg.size()-to_reg.size());
+                }
+                else if (to_reg.size() > from_reg.size()) {
+                    p = to_reg.begin();
+                    egsWarning("\nEGS_PhspScoring: Mismatch in no. of exit/entry regions.\n"
+                               "Will only score for matched pairs.\n");
+                    p += from_reg.size();
+                    to_reg.erase(p,p+to_reg.size()-from_reg.size());
+                }
+                //now go through and look for exit region = entry region
+                int i=0;
+                while (i<from_reg.size()) {
+                    if (from_reg[i]==to_reg[i]) {
+                        egsInformation("\nEGS_PhspScoring: Cannot have entry region = exit region (reg no. %d)\n",from_reg[i]);
+                        egsInformation("Will delete this pair\n");
+                        p=from_reg.begin()+i;
+                        p1 =to_reg.begin()+i;
+                        from_reg.erase(p);
+                        to_reg.erase(p1);
+                    }
+                    else {
+                        //advance counter
+                        i++;
+                    }
+                }
+            }
+        }
+        //now get common inputs for both scoring methods
+        if (input->getInput("output format", str) < 0) {
+            egsInformation("EGS_PhspScoring: No input for output format type.  Will default to EGSnrc.\n");
+            phspouttype = 0;
+        }
+        else {
+            vector<string> allowed_oformat;
+            allowed_oformat.push_back("EGSnrc");
+            allowed_oformat.push_back("IAEA");
+            phspouttype = input->getInput("output format", allowed_oformat, -1);
+            if (phspouttype < 0) {
+                egsFatal("\nEGS_PhspScoring: Invalid output format.\n");
+            }
+            //see if the user wants to specify constant X/Y/Z for IAEA format
+            if (phspouttype == 1) {
+                int err02 = input->getInput("constant X",xyzconst[0]);
+                int err03 = input->getInput("constant Y",xyzconst[1]);
+                int err04 = input->getInput("constant Z",xyzconst[2]);
+                if (!err02) {
+                    xyzisconst[0] = true;
+                }
+                if (!err03) {
+                    xyzisconst[1] = true;
+                }
+                if (!err04) {
+                    xyzisconst[2] = true;
+                }
+                //see if user wants to score mu (if available)
+                //default is not to score
+                if (!input->getInput("score mu", str)) {
+                    vector<string> allowed_muscore;
+                    allowed_muscore.push_back("no");
+                    allowed_muscore.push_back("yes");
+                    imuscore = input->getInput("score mu",allowed_muscore,-1);
+                    if (imuscore < 0) {
+                        egsWarning("\nEGS_PhspScoring: Invalid input for mu scoring.  Will not score mu.\n");
+                        imuscore = 0;
+                    }
+                }
+            }
+        }
+        if (phspouttype == 0) {
+            //see if user wants to score multiple crossers
+            if (!input->getInput("score multiple crossers", str)) {
+                vector<string> allowed_scoremc;
+                allowed_scoremc.push_back("no");
+                allowed_scoremc.push_back("yes");
+                iscoremc = input->getInput("score multiple crossers",allowed_scoremc,-1);
+                if (iscoremc < 0) {
+                    egsWarning("\nEGS_PhspScoring: Invalid input for score multiple crossers.  Will not score.\n");
+                    iscoremc = 0;
+                }
+            }
+        }
+        if (input->getInput("output directory",outdir) < 0) {
+            outdir="";
+        }
+        if (input->getInput("particle type", str) < 0) {
+            egsInformation("EGS_PhspScoring: No input for particle type.  Will score all.\n");
+            ptype = 0;
+        }
+        else {
+            //get particle type
+            vector<string> allowed_ptype;
+            allowed_ptype.push_back("all");
+            allowed_ptype.push_back("photons");
+            allowed_ptype.push_back("charged");
+            ptype = input->getInput("particle type",allowed_ptype,-1);
+            if (ptype < 0) {
+                egsFatal("\nEGS_PhspScoring: Invalid particle type.\n");
+            }
+        }
+
+        //=================================================
+
+        /* Setup phsp scoring object with input parameters */
+        EGS_PhspScoring *result = new EGS_PhspScoring("",f);
+        result->setName(input);
+        if (stype==0) {
+            result->setGeom(phspgeom);
+        }
+        else if (stype==1) {
+            result->setEntryExitReg(from_reg,to_reg);
+        }
+        result->setOType(phspouttype);
+        result->setXYZconst(xyzisconst,xyzconst);
+        result->setOutDir(outdir);
+        result->setParticleType(ptype);
+        result->setScoreDir(sdir);
+        result->setMuScore(imuscore);
+        result->setScoreMC(iscoremc);
+        return result;
+    }
+}
diff --git a/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/egs_phsp_scoring.h b/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/egs_phsp_scoring.h
new file mode 100644
index 000000000..c6966f88a
--- /dev/null
+++ b/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/egs_phsp_scoring.h
@@ -0,0 +1,551 @@
+/*
+###############################################################################
+#
+#  EGSnrc egs++ phase space scoring object headers
+#  Copyright (C) 2015 National Research Council Canada
+#
+#  This file is part of EGSnrc.
+#
+#  EGSnrc is free software: you can redistribute it and/or modify it under
+#  the terms of the GNU Affero General Public License as published by the
+#  Free Software Foundation, either version 3 of the License, or (at your
+#  option) any later version.
+#
+#  EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
+#  WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+#  FOR A PARTICULAR PURPOSE.  See the GNU Affero General Public License for
+#  more details.
+#
+#  You should have received a copy of the GNU Affero General Public License
+#  along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
+#
+###############################################################################
+#
+#  Author:          Blake Walters, 2018
+#
+#  Contributors:
+#
+###############################################################################
+*/
+
+
+/*! \file egs_phsp_scoring.h
+ *  \brief A phase space scoring ausgab object
+ *  \BW
+ */
+
+#ifndef EGS_PHSP_SCORING_
+#define EGS_PHSP_SCORING_
+
+#include "egs_ausgab_object.h"
+#include "egs_application.h"
+#include "egs_scoring.h"
+#include "egs_base_geometry.h"
+
+#ifdef WIN32
+
+    #ifdef BUILD_PHSP_SCORING_DLL
+        #define EGS_PHSP_SCORING_EXPORT __declspec(dllexport)
+    #else
+        #define EGS_PHSP_SCORING_EXPORT __declspec(dllimport)
+    #endif
+    #define EGS_PHSP_SCORING_LOCAL
+
+#else
+
+    #ifdef HAVE_VISIBILITY
+        #define EGS_PHSP_SCORING_EXPORT __attribute__ ((visibility ("default")))
+        #define EGS_PHSP_SCORING_LOCAL  __attribute__ ((visibility ("hidden")))
+    #else
+        #define EGS_PHSP_SCORING_EXPORT
+        #define EGS_PHSP_SCORING_LOCAL
+    #endif
+
+#endif
+
+/*! \brief A phase space scoring object: header
+
+\ingroup AusgabObjects
+
+This ausgab object generates phase space data for particles using one of two methods:
+1. Scores particles exiting and/or entering all surfaces of a predefined geometry.
+   The geometry can be a component of the simulation geometry or coincident with a surface
+   in the simulation geometry.
+2. Scores particles on exiting one user-specified region and entering another.
+   The user can specify multiple exit/entry region pairs.
+
+Phase space data can be scored in one of 2 possible formats:
+
+EGSnrc format: E,x,y,u,v,wt,latch
+IAEA format: iq,E,[x],[y],[z],u,v,wt,latch,[mu]
+
+Note that in IAEA format, the user has the option of specifying a fixed x, y, and/or z coordinate
+of the scoring plane/line/point, in which case the fixed coordinates shall not be scored for each
+particle but will be specified in the header (.IAEAheader) file.  Also, in this format, the
+user has the option of scoring the synchronization parameter, mu, passed from the source.
+
+Can be used in any C++ user code by entering the proper input block in
+the ausgab object definition block.
+
+\verbatim
+:start ausgab object definition:
+      :start ausgab object:
+          library                  = egs_phsp_scoring
+          name                     = some_name
+          output format            = EGSnrc or IAEA
+          constant X               = X value (cm) at which all particles are scored (IAEA format only)
+          constant Y               = Y value (cm) at which all particles are scored (IAEA format only)
+          constant Z               = Z value (cm) at which all particles are scored (IAEA format only)
+          particle type            = all, photons, or charged
+          score mu                 = yes or no [default] (IAEA format only)
+          score multiple crossers  = no (default) or yes (EGSnrc format only)
+          output directory         = name of output directory
+
+     and one of two methods of scoring particles:
+
+     Method 1: score particles on entry to/exit from a predefined geometry
+          phase space geometry     = name of previously defined geometry
+          score particles on       = entry, exit, entry and exit [default]
+
+     Method 2: score particles on exiting one region and entering another
+          from regions             = list of exit region numbers
+          to regions               = list of entry region numbers
+
+      :stop ausgab object:
+:stop ausgab object definition:
+\endverbatim
+
+Phase space data is output to the file some_name.egsphsp1 (EGSnrc format)
+or some_name.1.IAEAphsp and some_name.1.IAEAheader (IAEA format)
+Note that if the user specifies constant X/Y/Z, then particles are all assumed to be scored at the
+same X/Y/Z with this(ese) values output to .IAEAheader instead of being output for each particle
+to the .IAEAphsp file.
+
+Particles of the type indicated by the "particle type" input are scored.  If this input
+is omitted, then all particles are scored.
+
+In IAEA-format phase space files, the user can opt to score the synchronization parameter, mu.
+This option will automatically be turned off if the parameter is not available from the
+source.  Currently, this parameter can only be passed from egs_beam_source (only if the accelerator
+includes synchronized CMs), iaea_phsp_source (if scored using a BEAMnrc simulation with
+synchronized CMs or scored using this ausgab object with mu scoring turned on) and
+egs_dynamic_source (always available).
+
+The default is to not score multiple crossers (and their descendents) and, indeed, this
+is, by definition, the protocol for IAEA format phase space files.  However, if the user is scoring
+data in EGSnrc format, then they have the option to include multiple crossers and their
+descendents.  Note that the marker for a particle having been scored is to set bit 31 of the
+particle's latch high.  Thus, the marker is associated with the particle, not the scoring
+object.  This has implications if the user wishes to have multiple phase space scoring
+objects in one run: a particle scored by one phase space scoring object will not subsequently
+be scored by another, unless the user has set "score multiple crossers = yes" in the
+latter.
+
+If "output directory" is omitted or left blank then the output directory defaults to the
+application directory (i.e. $EGS_HOME/appname).
+
+When using a phase space scoring geometry, particles can be scored on entering the phase space geometry,
+exiting the geometry, or both (the default).  Be aware of how the "inside" and "outside" of the geometry
+are defined when using this option.
+
+When using pairs of exit/entry regions, the number of regions specified by the "to regions" input
+must equal that of the "from regions" input.  If the exit region number is equal to the entry
+region number, then the pair is deleted prior to the run.
+
+A note on parallel runs:
+If a phase space file is being written during a parallel run, then each job, i, outputs its phase
+space data to some_name_wi.[egsphsp1][.1.IAEAheader/phsp]. Thus, the naming scheme is the same as
+that for other output files from a parallel run.  These phase space files are not added automatically
+when the results of a parallel run are combined.  The user must either use the addphsp tool or
+program their own concatenation routine.
+
+Example:
+The following example input illustrates the two phase space scoring methods.  In both
+cases, phase space is scored in planes perpendicular to the axis of a water/air cylinder.
+\verbatim
+:start geometry definition:
+      :start geometry:
+          library = egs_planes
+          type = EGS_Zplanes
+          positions = 21.0
+          name = scoreplane
+      :stop geometry:
+
+      :start geometry:
+          library = egs_planes
+          type = EGS_Zplanes
+          positions = 0.0 10.0 21.0
+          name = cutplane
+      :stop geometry:
+
+      :start geometry:
+          library = egs_cylinders
+          type = EGS_ZCylinders
+          name = cylinders
+          radii = 1.0 5.0 15.0
+          midpoint = 0
+          :start media input:
+             media = H2O521ICRU AIR521ICRU
+             set medium = 0 0
+             set medium = 1 1
+             set medium = 2 0
+          :stop media input:
+      :stop geometry:
+
+      :start geometry:
+           name = maingeom
+           library = egs_cdgeometry
+           base geometry = cutplane
+           set geometry = 0 cylinders
+           set geometry = 1 cylinders
+      :stop geometry:
+
+      simulation geometry = maingeom
+  :stop geometry definition:
+
+  :start ausgab object definition:
+      :start ausgab object:
+        library = egs_phsp_scoring
+        name = test
+        from regions = 0 2
+        to regions = 3 5
+        output format = IAEA
+        constant Z = 10.0
+        particle type = all
+        score mu = no
+      :stop ausgab object:
+
+      :start ausgab object:
+        library = egs_phsp_scoring
+        name = test2
+        phase space geometry = scoreplane
+        output format = EGSnrc
+        particle type = all
+        score particles on = entry
+        score multiple crossers = yes
+      :stop ausgab object:
+:stop ausgab object definition:
+\endverbatim
+The first ausgab object, "test," uses scoring method 2 (exit/entry region pairs) to
+score phase space data, in IAEA format, for all particles crossing a plane perpendicular
+to the cylinder and at the mid-point of the cylinder (Z=10 cm).  Data will be scored
+for particles within the innermost cylinder (exit/entry regions 0/3) and within the outermost
+annulus (exit/entry regions 2/5).  Note that the constant Z value of the scoring plane is input,
+so data will be in 2D format.  Also note that IAEA format, by convention, does not score multiple
+crossers.  The output will be to test.1.IAEAheader (header file) and test.1.IAEAphsp (phase space data).
+
+The second ausgab object, "test2," uses scoring method 1 (predefined scoring geometry) to
+score phase space data, in EGSnrc format, for all particles crossing a single plane coincident with
+the bottom (Z=21 cm) of the cylinder.  The scoring plane is given by the geometry, "scoreplane,"
+which defines a single plane perpendicular to Z at Z=21 cm.  Recall that a single plane defines
+a single region on the +ve normal side of the plane.  In this example, a particle is considered "outside"
+scoreplane if Z < 21.0 and inside if Z > 21.0 and, thus, scoring particles on "entry" will collect
+data for all particles leaving out the bottom of the cylinder.  In this case, however, the specification
+of scoring direction is moot, since particles leaving the bottom of the cylinder are effectively leaving
+the simulation geometry.  Note that the option to score multiple crossers has been turned on.  Otherwise,
+particles previously scored in "test" would not be scored by "test2."  Phase space data will be output
+to test2.egsphsp1.
+*/
+
+class EGS_PHSP_SCORING_EXPORT EGS_PhspScoring : public EGS_AusgabObject {
+
+public:
+
+    EGS_PhspScoring(const string &Name="", EGS_ObjectFactory *f = 0);
+
+    ~EGS_PhspScoring();
+
+    int processEvent(EGS_Application::AusgabCall iarg) {
+        //only score if particle has correct charge
+        if (ocharge==0 || 1+abs(app->top_p.q)==ocharge) {
+            EGS_Vector x = app->top_p.x;
+            int ir = app->top_p.ir;
+            int latch = app->top_p.latch;
+            //only score if: 1) it has not been scored before or
+            //2) we are scoring multiple crossers (EGSnrc format only)
+            if (!(latch & bsmc()) || (oformat==0 && score_mc)) {
+                if (score_type==0) {  //using scoring geometry
+                    if (iarg == 0) {
+                        phsp_before = phsp_geom->isInside(x);
+                    }
+
+                    if (iarg == 5) {
+                        phsp_after = phsp_geom->isInside(x);
+                        if (phsp_after != phsp_before) {
+                            if (scoredir == 0 || (scoredir == 1 && phsp_after) ||
+                                    (scoredir == 2 && phsp_before)) {
+                                storeParticle(current_case);
+                            }
+                            //set bit 31 to flag this as having been scored
+                            latch = (latch | bsmc());
+                            app->setLatch(latch);
+                            return 0;
+                        }
+                    }
+                }
+                else if (score_type==1) { //pairs of exit/entry regions
+                    if (iarg == 0) {
+                        ir_before = ir;
+                    }
+
+                    if (iarg == 5) {
+                        ir_after = ir;
+                        if (from_to[ir_before].size()>0 && ir_before != ir_after) {
+                            for (int i=0; i< from_to[ir_before].size(); i++) {
+                                if (ir_after == from_to[ir_before][i]) {
+                                    storeParticle(current_case);
+                                    latch = (latch | bsmc());
+                                    app->setLatch(latch);
+                                    return 0;
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+        return 0;
+    };
+
+    int processEvent(EGS_Application::AusgabCall iarg, int ir) {
+        //same as above, we don't need the region no.
+        if (ocharge==0 || 1+abs(app->top_p.q)==ocharge) {
+            EGS_Vector x = app->top_p.x;
+            int latch = app->top_p.latch;
+            //only score if: 1) it has not been scored before or
+            //2) we are scoring multiple crossers (EGSnrc format only)
+            if (!(latch & bsmc()) || (oformat==0 && score_mc)) {
+                if (score_type==0) {  //using scoring geometry
+                    if (iarg == 0) {
+                        phsp_before = phsp_geom->isInside(x);
+                    }
+
+                    if (iarg == 5) {
+                        phsp_after = phsp_geom->isInside(x);
+                        if (phsp_after != phsp_before) {
+                            if (scoredir == 0 || (scoredir == 1 && phsp_after) ||
+                                    (scoredir == 2 && phsp_before)) {
+                                storeParticle(current_case);
+                            }
+                            //set bit 31 to flag this as having been scored
+                            latch = (latch | bsmc());
+                            app->setLatch(latch);
+                            return 0;
+                        }
+                    }
+                }
+                else if (score_type==1) { //pairs of exit/entry regions
+                    if (iarg == 0) {
+                        ir_before = ir;
+                    }
+
+                    if (iarg == 5) {
+                        ir_after = ir;
+                        if (from_to[ir_before].size()>0 && ir_before != ir_after) {
+                            for (int i=0; i< from_to[ir_before].size(); i++) {
+                                if (ir_after == from_to[ir_before][i]) {
+                                    storeParticle(current_case);
+                                    latch = (latch | bsmc());
+                                    app->setLatch(latch);
+                                    return 0;
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+        return 0;
+    };
+
+    bool needsCall(EGS_Application::AusgabCall iarg) const {
+        if (iarg == 0 || iarg == 5) {
+            return true;
+        }
+        else {
+            return false;
+        }
+    };
+
+    //below gets called from startNewShower if current_case != last_case
+    //don't update last_case yet
+    void setCurrentCase(EGS_I64 ncase) {
+        current_case=ncase;
+    };
+
+    void setGeom(EGS_BaseGeometry *phspgeom) {
+        score_type=0;
+        phsp_geom = phspgeom;
+    }
+
+    void setEntryExitReg(const vector <int> from_reg, const vector <int> to_reg) {
+        score_type=1;
+        fromreg=from_reg;
+        toreg=to_reg;
+    }
+
+    void setOType(const int phspouttype) {
+        oformat = phspouttype;
+    }
+
+    //set output directory
+    void setOutDir(const string outdir) {
+        phspoutdir =  outdir;
+    }
+
+    void setParticleType(const int ptype) {
+        ocharge = ptype;
+    }
+
+    void setScoreDir(const int sdir) {
+        scoredir = sdir;
+    }
+
+    void setMuScore(const int imuscore) {
+        if (imuscore == 1) {
+            score_mu = true;
+        }
+        else {
+            score_mu = false;
+        }
+    }
+
+    void setScoreMC(const int iscoremc) {
+        if (iscoremc==1) {
+            score_mc = true;
+        }
+        else {
+            score_mc = false;
+        }
+    }
+
+    //method below pertains only to IAEA format
+    //set array element 0/1/2 of xyz_is_constant equal to true
+    //if scoring at a constant X/Y/Z value and store the
+    //constant value in element 0/1/2 of array xyzscore
+    void setXYZconst(bool xyzisconst[3], float xyzconst[3]) {
+        for (int i=0; i<3; i++) {
+            xyz_is_const[i] = xyzisconst[i];
+            xyzscore[i]=xyzconst[i];
+        }
+    }
+
+    void storeParticle(EGS_I64 ncase);
+
+    int flushBuffer() const;
+
+    void openPhspFile() const;
+
+    void setApplication(EGS_Application *App);
+
+    void reportResults();
+
+    bool storeState(ostream &data) const;
+    bool setState(istream &data);
+    bool addState(istream &data);
+
+protected:
+
+    struct Particle {
+        int  q, latch;
+        EGS_Float E, x, y, z, u, v, w, wt, mu;
+    };
+
+    //functions, struct and variables used to write EGSnrc format phsp files
+    static unsigned int bclr() {
+        return ~((1 << 30) | (1 << 29));
+    }
+    static unsigned int bsqe() {
+        return (1 << 30);
+    }
+    static unsigned int bsqp() {
+        return (1 << 29);
+    }
+    struct egs_phsp_write_struct {
+        int   latch;
+        float E;
+        float x,y;
+        float u,v;
+        float wt;
+        egs_phsp_write_struct() {};
+        egs_phsp_write_struct(const Particle &p) {
+            latch = (p.latch & bclr());
+            if (p.q == -1) {
+                latch = (latch | bsqe());
+            }
+            else if (p.q == 1) {
+                latch = (latch | bsqp());
+            }
+            E = p.E;
+            x = p.x;
+            y = p.y;
+            u = p.u;
+            v = p.v;
+            wt = p.w >= 0 ? p.wt : -p.wt;
+        };
+    };
+    mutable fstream phsp_file; //output file -- mutable so we can write to it during storeState
+    EGS_I64 count; //total no. of particles in file
+    EGS_I64 countg; //no. of photons in file
+    float emax; //max. k.e. of particles in phsp file
+    float emin; //min. k.e. of charged particles in file
+    bool score_mc; //set to true to score multiple crossers and their descendents
+
+    //variables specific to IAEA format
+    mutable int iaea_id; //file id--mutable so we can write to it during storeState
+    int latch_ind; //IAEA id for latch variable (set to 2)
+    int mu_ind; //IAEA id for mu--no ID for now so set to generic float (0)
+    int iaea_n_extra_float, iaea_n_extra_long; //no. of extra floats, ints
+    int iaea_i_latch, iaea_i_mu; //indices of indicated variables in arrays
+    int iaea_q_type[3]; //easy conversion from q to iaea type
+    bool xyz_is_const[3]; //set to true if scoring at constant X/Y/Z
+    float xyzscore[3]; //constant X/Y/Z scoring values
+    int len; //length of name
+    char *phsp_fname_char; //need file name in char format
+    bool score_mu; //set to true if scoring mu
+    float pmu; //mu value associated with particle
+
+    //back to variables common to EGSnrc and IAEA formats
+
+    Particle *p_stack; //the stored particle stack
+
+    //below used to set bit 31 to denote the particle has been scored
+    static unsigned int bsmc() {
+        return (1 << 31);
+    }
+
+    int store_max; //max. no. of particles to store in p_stack
+    mutable int phsp_index; //index in p_stack array -- mutable so we can change it in storeState
+    mutable bool first_flush; //first time writing to file in this run -- mutable so we can change it in storeState
+
+    bool is_restart; //true if this is a restart
+    EGS_I64 countprev; //no. of particles written to file before restart
+
+    EGS_I64    last_case;   //last primary history scored
+    EGS_I64    current_case; //current primary history
+
+    int oformat;           //0 for EGSnrc format, 1 for IAEA format
+
+    int ocharge;           //particle type for output: 0--all; 1--photons; 2--charged particles
+
+    string phsp_fname; //name of phase space file
+
+    string phspoutdir; //output directory
+
+    int score_type;        //0 if scoring on exiting/entering a phase space geometry
+    //1 if using pairs of exit/entry regions
+
+    //for method 1: scoring using predefined geometry
+    EGS_BaseGeometry *phsp_geom; //geometry on entrance to/exit from which phase space data is scored
+    int scoredir;           //scoring direction: 0--on entry and exit; 1--on entry; 2--on exit
+    bool phsp_before; //true if inside scoring geometry before step
+    bool phsp_after;  //true if inside scoring geometry after step
+
+    //for method 2: scoring using exit/entry region pairs
+    vector <int> fromreg;          //array of exit regions
+    vector <int> toreg;            //array of entry regions
+    vector <vector <int> > from_to;   //from a given global exit region, an array of possible entry regions
+    int ir_before, ir_after;       //reg. no. before and after step
+};
+
+#endif
diff --git a/HEN_HOUSE/egs++/egs_advanced_application.cpp b/HEN_HOUSE/egs++/egs_advanced_application.cpp
index f73280b2b..78d729af8 100644
--- a/HEN_HOUSE/egs++/egs_advanced_application.cpp
+++ b/HEN_HOUSE/egs++/egs_advanced_application.cpp
@@ -1198,6 +1198,14 @@ void EGS_AdvancedApplication::setRadiativeSplitting(const EGS_Float &nsplit) {
     the_egsvr->nbr_split = nsplit;
 }
 
+//************************************************************
+// Utility function for ausgab phase space scoring objects
+//************************************************************
+void EGS_AdvancedApplication::setLatch(int latch) {
+    int np = the_stack->np-1;
+    the_stack->latch[np] = latch;
+}
+
 extern __extc__ void egsHowfar() {
     CHECK_GET_APPLICATION(app,"egsHowfar()");
     int np = the_stack->np-1;
diff --git a/HEN_HOUSE/egs++/egs_advanced_application.h b/HEN_HOUSE/egs++/egs_advanced_application.h
index ea95e1acf..b627ef340 100644
--- a/HEN_HOUSE/egs++/egs_advanced_application.h
+++ b/HEN_HOUSE/egs++/egs_advanced_application.h
@@ -211,6 +211,12 @@ class APP_EXPORT EGS_AdvancedApplication : public EGS_Application {
     void setEdep(EGS_Float edep);
     EGS_Float getEcut();
     EGS_Float getPcut();
+
+    //************************************************************
+    // Utility function for ausgab phase space scoring objects
+    //************************************************************
+    void setLatch(int latch);
+
     /* Needed by some sources */
     EGS_Float getRM();
     /* Turn ON/OFF radiative splitting */
diff --git a/HEN_HOUSE/egs++/egs_application.h b/HEN_HOUSE/egs++/egs_application.h
index 5e8f52cab..527513719 100644
--- a/HEN_HOUSE/egs++/egs_application.h
+++ b/HEN_HOUSE/egs++/egs_application.h
@@ -674,6 +674,18 @@ class EGS_EXPORT EGS_Application {
         geometry->getLabelRegions(str, regs);
     }
 
+    /*! \brief Returns the value of the \a mu synchronization parameter
+
+      The parameter, \a mu, is a random number on \a [0,1) associated with each
+      primary history and is retrieved from \a source.  It can be used to
+      synchronize geometric parameters throughout a simulation.  If \a mu is
+      not available in \a source (i.e., the \a getMu function has not been
+      reimplemented in \a source), then this returns -1.
+     */
+    EGS_Float getMU() {
+        return source->getMu();
+    }
+
     /*! \brief User scoring function for accumulation of results and VRT implementation
 
       This function first calls the processEvent() method of the ausgab objects
@@ -1116,6 +1128,12 @@ class EGS_EXPORT EGS_Application {
         return -1.0;
     };
     virtual void setRadiativeSplitting(const EGS_Float &nsplit) {};
+
+    //************************************************************
+    // Utility function for ausgab phase space scoring objects
+    //************************************************************
+    virtual void setLatch(int latch) {};
+
 };
 
 #define APP_MAIN(app_name) \
diff --git a/HEN_HOUSE/iaea_phsp/iaea_header.cpp b/HEN_HOUSE/iaea_phsp/iaea_header.cpp
index 1181623b4..c06c49558 100644
--- a/HEN_HOUSE/iaea_phsp/iaea_header.cpp
+++ b/HEN_HOUSE/iaea_phsp/iaea_header.cpp
@@ -687,11 +687,11 @@ int iaea_header_type::check_byte_order()
   // printf("\n \t %x %x %x %x", pf[0],pf[1],pf[2],pf[3]);
   if(pf[0] == 0 && pf[3] != 0)
   {
-    printf("\nByte order: INTEL / ALPHA,LINUX -> LITLE_ENDIAN \n");
+    //printf("\nByte order: INTEL / ALPHA,LINUX -> LITLE_ENDIAN \n");
     return(LITTLE_ENDIAN);
   }else if(pf[0] != 0 && pf[3] == 0)
   {
-    printf("\nByte order: OTHER (SGI,SUN-SOLARIS) -> BIG_ENDIAN \n ");
+    //printf("\nByte order: OTHER (SGI,SUN-SOLARIS) -> BIG_ENDIAN \n ");
     return(BIG_ENDIAN);
   }
   else