diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5139f77ed40bb2078ad1761840c581f2d4f087e3..d60d4935d7008900e07f8788686a5e0a7cb7fea6 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -5,7 +5,7 @@
 ################################################################################
 cmake_minimum_required (VERSION 2.8)
 project (pythia6m)
-set (PYTHIA6M_MAJOR_VERSION 1)
+set (PYTHIA6M_MAJOR_VERSION 2)
 set (PYTHIA6M_MINOR_VERSION 0)
 set (PYTHIA6M_PATCH_VERSION 0)
 set (PYTHIA6M_VERSION 
@@ -27,7 +27,8 @@ add_subdirectory(pythia6m/fmotion)
 add_subdirectory(pythia6m/gmc)
 add_subdirectory(pythia6m/pythia6)
 add_subdirectory(pythia6m/radgen)
-add_subdirectory(pythia6m/util)
+add_subdirectory(pythia6m/interface)
+add_subdirectory(pythia6m/core)
 
 add_subdirectory(program)
 
diff --git a/cmake/pythia6m-config.cmake.in b/cmake/pythia6m-config.cmake.in
index b5900e933673da81b5594a510d1fbcc2a2f5733e..e988dd1346365e0a7ea0374975e9268c26ed71dc 100644
--- a/cmake/pythia6m-config.cmake.in
+++ b/cmake/pythia6m-config.cmake.in
@@ -23,7 +23,8 @@ set(PYTHIA6M_LIBRARIES
   pythia6m_gmc
   pythia6m_pythia
   pythia6m_radgen
-  pythia6m_util)
+  pythia6m_interface
+  pythia6m_core)
 set(PYTHIA6M_LIBRARIES ${PYTHIA6M_LIBRARIES} ${PYTHIA6M_LIBRARIES})
 
 # Add the external libraries
diff --git a/program/CMakeLists.txt b/program/CMakeLists.txt
index cc43539a57b82554d820935a30ad2cd91981f690..eff5035b65610b51fe5317240111dd3395e60924 100644
--- a/program/CMakeLists.txt
+++ b/program/CMakeLists.txt
@@ -38,12 +38,14 @@ target_link_libraries(${EXE}
   pythia6m_gmc
   pythia6m_pythia
   pythia6m_radgen
-  pythia6m_util
+  pythia6m_interface
+  pythia6m_core
   pythia6m_fmotion
   pythia6m_gmc
   pythia6m_pythia
   pythia6m_radgen
-  pythia6m_util
+  pythia6m_interface
+  pythia6m_core
   ${ROOT_LIBRARIES} ${Boost_LIBRARIES} ${NANOCERNLIB_LIBRARIES})
 set_target_properties(${EXE} PROPERTIES VERSION ${PYTHIA6M_VERSION})
 
diff --git a/program/example-config.json b/program/example-config.json
index e5768c9f4707d8bdfbd0ebed0fabba6fc0f2668f..047aa9c560705823d5d65a19bcb4a8b5ddc2f6ee 100644
--- a/program/example-config.json
+++ b/program/example-config.json
@@ -1,30 +1,29 @@
 {
   "mc" : {
-    "type": "MC Data",
-    "target": {
-      "A": 1,
-      "Z": 1, 
-      "nucleon": "p"
-    },
-    "beam": {
-      "energy": 11
-    },
-    "Q2": {
-      "min": 0.9,
-      "max": 60
-    },
-    "x": {
-      "min": 0.002,
-      "max": 0.99
-    },
-    "y": {
-      "min": 0.05,
-      "max": 0.95 
+    "type": "pythia6m",
+    "write_lund": true,
+    "generator": {
+      "type": "pythia6",
+      "target": {
+        "A": 1,
+        "Z": 1, 
+        "nucleon": "p"
+      },
+      "beam": {
+        "energy": 11
+      },
+      "Q2": {
+        "min": 0.9,
+        "max": 60
+      },
+      "x": {
+        "min": 0.002,
+        "max": 0.99
+      },
+      "y": {
+        "min": 0.05,
+        "max": 0.95 
+      }
     }
-  },
-  "output": {
-    "type": "MC Output",
-    "root": true,
-    "text": true
   }
 }
diff --git a/program/pythia6m.cc b/program/pythia6m.cc
index e0da8111d63e2f9634df11bb090cfceb322dc9c8..efd49ea80284eed1d35a662c437f05fbf9f49ccd 100644
--- a/program/pythia6m.cc
+++ b/program/pythia6m.cc
@@ -1,421 +1,60 @@
-#include <pythia6m/util/configuration.hh>
-#include <pythia6m/util/exception.hh>
-#include <pythia6m/util/framework.hh>
-#include <pythia6m/util/logger.hh>
-#include <pythia6m/util/assert.hh>
-
-#include <pythia6m/gmc/gen_set.h>
-#include <pythia6m/gmc/gmc.h>
-#include <pythia6m/gmc/mc_set.h>
-#include <pythia6m/gmc/mcevnt.h>
-#include <pythia6m/gmc/random.h>
-
-#include <pythia6m/fmotion/fmotion.h>
-
-#include <pythia6m/pythia6/pythia6.h>
-
-#include <cstring>
+#include <Math/Vector4D.h>
+#include <TFile.h>
 #include <fstream>
 #include <memory>
-
-#include <TFile.h>
-#include <TLorentzVector.h>
-#include <TTree.h>
+#include <pythia6m/core/configuration.hh>
+#include <pythia6m/core/exception.hh>
+#include <pythia6m/core/framework.hh>
+#include <pythia6m/core/logger.hh>
+#include <pythia6m/core/progress_meter.hh>
+#include <pythia6m/interface/event.hh>
+#include <pythia6m/interface/event_out.hh>
+#include <pythia6m/interface/pythia6m_generator.hh>
 
 using namespace pythia6m;
 
-GeneratorSettings glb_gen;
-MCSettings glb_mc;
-EventData glb_event;
-
-// global for output configuration for ...sigh... consistency
-struct output_cfg {
-  TFile* root;
-  std::unique_ptr<std::ofstream> text;
-};
-output_cfg glb_out;
-
-// global event buffer...
-struct part_buffer {
-  int32_t n_part;
-  std::vector<int32_t> index;
-  std::vector<int32_t> status;
-  std::vector<int32_t> type;
-  std::vector<int32_t> parent;
-  std::vector<int32_t> daughter; // this is the first daughter, or zero
-  std::vector<int32_t> charge;
-  std::vector<double> mom;
-  std::vector<double> px;
-  std::vector<double> py;
-  std::vector<double> pz;
-  std::vector<double> E;
-  std::vector<double> mass;
-  std::vector<double> vx;
-  std::vector<double> vy;
-  std::vector<double> vz;
-  std::vector<bool> init;
-  std::vector<bool> scat;
-  std::vector<bool> virt;
-  std::vector<bool> lund;
-
-  void clear() {
-    n_part = 0;
-    index.clear();
-    status.clear();
-    type.clear();
-    parent.clear();
-    daughter.clear();
-    charge.clear();
-    mom.clear();
-    px.clear();
-    py.clear();
-    pz.clear();
-    E.clear();
-    mass.clear();
-    vx.clear();
-    vy.clear();
-    vz.clear();
-    init.clear();
-    scat.clear();
-    virt.clear();
-    lund.clear();
-  }
-
-  void add_all() {
-    for (int idx = 0; idx < glb_pyJets.n; ++idx) {
-      add(idx);
-    }
-  }
-
-  void add(int32_t idx) {
-    n_part += 1;
-    int ks = glb_pyJets.K[0][idx];
-    int kf = glb_pyJets.K[1][idx];
-    index.push_back(idx + 1);
-    status.push_back(ks);
-    type.push_back(kf);
-    parent.push_back(glb_pyJets.K[2][idx]);
-    daughter.push_back(glb_pyJets.K[3][idx]);
-    charge.push_back(PYCHGE(kf)/3);
-    TLorentzVector ptmp{glb_pyJets.P[0][idx], glb_pyJets.P[1][idx],
-                        glb_pyJets.P[2][idx], glb_pyJets.P[3][idx]};
-    mom.push_back(ptmp.Vect().Mag());
-    px.push_back(ptmp.X());
-    py.push_back(ptmp.Y());
-    pz.push_back(ptmp.Z());
-    E.push_back(ptmp.E());
-    mass.push_back({glb_pyJets.P[4][idx]});
-    // mm -> cm
-    vx.push_back(glb_pyJets.V[0][idx]/10);
-    vy.push_back(glb_pyJets.V[1][idx]/10);
-    vz.push_back(glb_pyJets.V[2][idx]/10);
-    if (glb_event.scatIdx == idx) {
-      init.push_back(false);
-      scat.push_back(true);
-      virt.push_back(false);
-      lund.push_back(false);
-    } else if (ks == 21) {
-      init.push_back(true);
-      scat.push_back(false);
-      virt.push_back(false);
-      lund.push_back(false);
-    } else if (ks >= 1 && ks <= 10 && ks != 2 && ks != 3) {
-      init.push_back(false);
-      scat.push_back(false);
-      virt.push_back(false);
-      lund.push_back(true);
-    } else {
-      init.push_back(false);
-      scat.push_back(false);
-      virt.push_back(true);
-      lund.push_back(false);
-    }
-  }
-};
-part_buffer glb_parts;
-
-void conf_generator(const configuration& conf) {
-  LOG_INFO("conf", "Configuring generator");
-
-  // User control switch MSEL (manual section 9.2, also 8.3.1)
-  // MSEL=2 additionally enables the elastic and diffactive components of the
-  // VMD and GVMD parts
-  glb_gen.PhotoType = 2;
-
-  // Let LUND decay some of the semi-stable particles
-  glb_gen.DecayPiK = 0;
-  glb_gen.DecayLamKs = 1;
-
-  // RADGEN options, currently disabled
-  glb_gen.radgenGenLUT = 0;
-  glb_gen.radgenUseLUT = 0;
-  sprintf(glb_gen.GenRad, "NO\n");
-
-  // activate Fermi momentum? (experimental)
-  glb_gen.enableFM = 1;
-  glb_gen.FMNBins = 1000; // Number of bins for FM discretization
-  glb_gen.FMCutOff = 1.5; // Max FM (in GeV)
-
-  // Parameterizations for F2Pythia and R
-  // (I believe these are only relevant when running with RADGEN)
-  sprintf(glb_gen.FStruct, "F2PY");
-  sprintf(glb_gen.R, "1990");
-
-  // PYTHIA model to run with:
-  //  * RAD for realistic DIS
-  //  * GAM for real photon beam
-  sprintf(glb_gen.PyModel, "RAD");
-}
-void conf_mc(const configuration& conf) {
-  LOG_INFO("conf", "Configuring MC");
-
-  glb_mc.RunNo = conf.get<int>("run");
-  glb_mc.NEvents = conf.get<int>("events");
-  glb_mc.TarA = conf.get<int>("target/A");
-  glb_mc.TarA = conf.get<int>("target/Z");
-  glb_mc.TargetNucleon[0] = conf.get<char>("target/nucleon");
-  glb_mc.EneBeam = conf.get<double>("beam/energy");
-  glb_mc.Q2Min = conf.get<double>("Q2/min");
-  glb_mc.Q2Max = conf.get<double>("Q2/max");
-  glb_mc.xMin = conf.get<double>("x/min");
-  glb_mc.xMax = conf.get<double>("x/max");
-  glb_mc.yMin = conf.get<double>("y/min");
-  glb_mc.yMax = conf.get<double>("y/max");
-
-  LOG_INFO("conf", "Will generate " + std::to_string(glb_mc.NEvents) +
-                       " events for run " + std::to_string(glb_mc.RunNo));
-
-  // hardcoded electron beam for now
-  sprintf(glb_mc.BeamParType, "ELEC");
-  // hardcoded Beam and target Polarization
-  // needed because RADGEN can also run with polarized generators
-  glb_mc.PBValue = 0;
-  glb_mc.PTValue = 0;
-  sprintf(glb_mc.PBeam, "L\n");
-  sprintf(glb_mc.PTarget, "L\n");
-
-  // reconsider fermi motion when running for A=1 options
-  if (glb_mc.TarA == 1) {
-    glb_gen.enableFM = 0;
-  }
-}
-// init the RNG
-void init_random() {
-  LOG_INFO("init",
-           "initializing the random generator with seed (run number): " +
-               std::to_string(glb_mc.RunNo));
-  int dummy1{0};
-  int dummy2{0};
-  const char* command{" "};
-  RNDMQ(dummy1, dummy2, glb_mc.RunNo, command);
-}
-// init the event structure
-void init_event() {
-  LOG_INFO("init", "Initializing the event common blocks");
-  // beam particle
-  if (!strncmp(glb_mc.BeamParType, "POSI", 4)) {
-    glb_event.beamLType = -11;
-  } else if (!strncmp(glb_mc.BeamParType, "ELEC", 4)) {
-    glb_event.beamLType = 11;
-  } else if (!strncmp(glb_mc.BeamParType, "GAM", 3)) {
-    glb_event.beamLType = 22;
-  } else {
-    throw exception("Unrecognized beam particle. This should NEVER HAPPEN!");
-  }
-
-  // Init FM engine
-  if (glb_gen.enableFM) {
-    INIT_FM(glb_mc.TarZ, glb_mc.TarA, glb_gen.FMCutOff, glb_gen.FMNBins);
-  }
+int run_pythia(const configuration& conf, const std::string& output) {
 
-  glb_event.evNum = 0;
-  glb_event.nEvGen = 0;
-  glb_event.nEvGenPrev = 0;
-  glb_event.nErr = 0;
+  // get the run info
+  const size_t run = conf.get<size_t>("run");
+  const size_t events = conf.get<size_t>("events");
 
-  glb_event.weight = 1.f;
-  glb_event.Q2Gen = 0.f;
-  glb_event.nuGen = 0.f;
-  glb_event.xGen = 0.f;
-  glb_event.yGen = 0.f;
-  glb_event.W2Gen = 0.f;
-  glb_event.thetaGen = 0.f;
-  glb_event.phiGen = 0.f;
-  glb_event.EGen = 0.f;
-  glb_event.pGen = 0.f;
-  glb_event.pxGen = 0.f;
-  glb_event.pyGen = 0.f;
-  glb_event.pzGen = 0.f;
-  glb_event.vxGen = 0.f;
-  glb_event.vyGen = 0.f;
-  glb_event.vzGen = 0.f;
-}
-void init_radgen() {
-  // TODO
-}
-void init_gmc() {
-  LOG_INFO("init", "Initializing GMC/PYTHIA");
-  int ok;
-  GMC_INIT(ok);
-  if (!ok) {
-    throw exception("Failed to initialized GMC/PYTHIA");
+  // setup output files: 
+  //  * root file 
+  std::shared_ptr<TFile> ofile =
+      std::make_shared<TFile>((output + ".root").c_str(), "recreate");
+  ep_event_out ev_writer{ofile, "ep_event"};
+  //  * lund file (optional)
+  std::ofstream olund;
+  if (conf.get<bool>("write_lund")) {
+    olund.open(output + ".lund");
   }
-}
-void conf_output(const configuration& conf, const std::string& output) {
-  LOG_INFO("conf", "Configuring output");
-  if (conf.get<bool>("root")) {
-    glb_out.root =
-        new TFile{(output + ".root").c_str(), "recreate"};
-  }
-  if (conf.get<bool>("text")) {
-    glb_out.text = std::make_unique<std::ofstream>((output + ".lund").c_str());
-  }
-  tassert(glb_out.root || glb_out.text, "Need at least one output format");
-}
 
-void clear_event_kinematics() {
-  glb_event.scatIdx = -1;
-  glb_event.Q2Gen = 0.f;
-  glb_event.nuGen = 0.f;
-  glb_event.xGen = 0.f;
-  glb_event.yGen = 0.f;
-  glb_event.W2Gen = 0.f;
-  glb_event.pxGen = 0.f;
-  glb_event.pyGen = 0.f;
-  glb_event.pzGen = 0.f;
-  glb_event.vxGen = 0.f;
-  glb_event.vyGen = 0.f;
-  glb_event.vzGen = 0.f;
-}
-
-void events() {
-  TTree* tree;
-  if (glb_out.root) {
-    glb_out.root->cd();
-    tree = new TTree("pythia6m", "pythia6m");
-    tree->Branch("event", &glb_event.evNum);
-    tree->Branch("evgen", &glb_event.nEvGen);
-    tree->Branch("xsec", &glb_pyPars.pari[0]);
-    tree->Branch("process", &glb_event.process);
-    tree->Branch("Q2", &glb_event.Q2Gen);
-    tree->Branch("nu", &glb_event.nuGen);
-    tree->Branch("x", &glb_event.xGen);
-    tree->Branch("y", &glb_event.yGen);
-    tree->Branch("W2", &glb_event.W2Gen);
-    tree->Branch("n_part", &glb_parts.n_part);
-    tree->Branch("index", &glb_parts.index);
-    tree->Branch("status", &glb_parts.type);
-    tree->Branch("type", &glb_parts.type);
-    tree->Branch("parent", &glb_parts.parent);
-    tree->Branch("daughter", &glb_parts.daughter);
-    tree->Branch("charge", &glb_parts.charge);
-    tree->Branch("mom", &glb_parts.mom);
-    tree->Branch("px", &glb_parts.px);
-    tree->Branch("py", &glb_parts.py);
-    tree->Branch("pz", &glb_parts.pz);
-    tree->Branch("E", &glb_parts.E);
-    tree->Branch("mass", &glb_parts.mass);
-    tree->Branch("vx", &glb_parts.vx);
-    tree->Branch("vy", &glb_parts.vy);
-    tree->Branch("vz", &glb_parts.vz);
-    tree->Branch("init", &glb_parts.init);
-    tree->Branch("scat", &glb_parts.scat);
-    tree->Branch("virt", &glb_parts.virt);
-    tree->Branch("lund", &glb_parts.lund);
-  }
+  // setup pythia6m
+  pythia6m_generator pygen{conf, "generator"};
 
-  // event loop
-  for (int iev{0}; iev < glb_mc.NEvents; ++iev) {
-    // clear the event kinematics
-    clear_event_kinematics();
-    glb_parts.clear();
+  // generate our events
+  progress_meter progress{events};
+  for (int iev = 0; iev < events; ++iev) {
+    progress.update();
 
     // generate one event
-    int nEvGen{0};
-    GMC_EVENT(nEvGen);
-
-    // update counters
-    glb_event.nEvGen += nEvGen;
-
-    //
-    // insert good event check here if needed
-    // 
-
-    // we now have a good event
-    glb_event.nEvGenPrev = glb_event.nEvGen;
-    glb_event.evNum = iev;
-
-    // talk to the human
-    if (!((iev + 1) % 5000)) {
-      LOG_INFO(__func__, "Generated Events: " +
-                             std::to_string(glb_event.nEvGen) +
-                             ", Cross Section: " +
-                             std::to_string(glb_pyPars.pari[0]) +
-                             ", Good Events: " +
-                             std::to_string(iev + 1));
-    }
-
-    // get all the particles
-    glb_parts.add_all();
-
-    // write output
-    if (glb_out.root && tree) {
-      tree->Fill();
+    const auto ev = pygen.generate<ep_event>();
+
+    // add additonal event selection here if needed
+    
+    // push event to the output writer
+    ev_writer.push(ev);
+
+    // also write lund file if needed
+    if (olund) {
+      ev.print(olund, [](const particle& p) {
+        return true; /*add selection criteria for output lines here */
+      });
     }
-    if (glb_out.text) {
-      char buf[1024];
-      int nparts = 0;
-      for (int ii = 0; ii < glb_parts.n_part; ++ii) {
-        if (glb_parts.lund[ii] || glb_parts.scat[ii]) {
-          nparts += 1;
-        }
-      }
-      sprintf(buf, "%4i %4i %4i %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e\n",
-              nparts, glb_mc.TarA, glb_mc.TarZ, 0., 0., glb_event.xGen,
-              glb_event.yGen, glb_event.W2Gen, glb_event.Q2Gen,
-              glb_event.nuGen);
-      *glb_out.text << buf;
-      for (int ii = 0; ii < glb_parts.n_part; ++ii) {
-        if (glb_parts.lund[ii] || glb_parts.scat[ii]) {
-          sprintf(buf,
-                  "%4i %4i %4i %4i %4i %4i %12.6e %12.6e %12.6e %12.6e %12.6e "
-                  "%12.6e %12.6e %12.6e\n",
-                  glb_parts.index[ii], glb_parts.charge[ii], 1, glb_parts.type[ii], glb_parts.parent[ii],
-                  glb_parts.daughter[ii], glb_parts.px[ii], glb_parts.py[ii],
-                  glb_parts.pz[ii], glb_parts.E[ii], glb_parts.mass[ii],
-                  glb_parts.vx[ii], glb_parts.vy[ii], glb_parts.vz[ii]);
-          *glb_out.text << buf;
-        }
-      }
-    }
-  }
-  if (glb_out.root) {
-    tree->Write();
-    // close the TFile
-    glb_out.root->Close();
   }
-  // print the run stats
-  GMC_RUNSTAT();
-
-  LOG_INFO("main",
-           "TOTAL Generated events: " + std::to_string(glb_event.nEvGen));
-  LOG_INFO("main", "TOTAL PYTHIA Cross Section [ubarn]: " +
-                       std::to_string(glb_pyPars.pari[0]));
-}
-
-int run_pythia(const configuration& conf, const std::string& output) {
-  LOG_INFO("main", "PYTHIA6M");
-  // configuration
-  conf_generator(conf);
-  conf_mc(conf);
-  conf_output(conf, output);
-  // initialization
-  init_random();
-  init_event();
-  init_gmc();
-
-  events();
 
+  // that's all!
   return 0;
 }
 
diff --git a/pythia6m/util/CMakeLists.txt b/pythia6m/core/CMakeLists.txt
similarity index 95%
rename from pythia6m/util/CMakeLists.txt
rename to pythia6m/core/CMakeLists.txt
index ef00ec205f339602d1c015c650d4dcb5d03a2b45..68977db5bf809b31ea8424306e584ef916dd4e75 100644
--- a/pythia6m/util/CMakeLists.txt
+++ b/pythia6m/core/CMakeLists.txt
@@ -1,7 +1,7 @@
 ################################################################################
 ## CMAKE Settings
 ################################################################################
-set (LIBRARY "pythia6m_util")
+set (LIBRARY "pythia6m_core")
 set (TARGETS ${TARGETS} ${LIBRARY} PARENT_SCOPE)
 
 ################################################################################
@@ -44,4 +44,4 @@ install(TARGETS ${LIBRARY}
   RUNTIME DESTINATION "${INSTALL_BIN_DIR}" COMPONENT bin
   LIBRARY DESTINATION "${INSTALL_LIB_DIR}" COMPONENT lib
   ARCHIVE DESTINATION "${INSTALL_LIB_DIR}" COMPONENT lib
-  PUBLIC_HEADER DESTINATION "${INSTALL_INCLUDE_DIR}/util" COMPONENT dev)
+  PUBLIC_HEADER DESTINATION "${INSTALL_INCLUDE_DIR}/core" COMPONENT dev)
diff --git a/pythia6m/util/assert.hh b/pythia6m/core/assert.hh
similarity index 88%
rename from pythia6m/util/assert.hh
rename to pythia6m/core/assert.hh
index 7b6022d55184b4856c8c449fcfb8a63d7abb6787..ac800c953cc6854a04c6aacc3fa8a601a61d77ae 100644
--- a/pythia6m/util/assert.hh
+++ b/pythia6m/core/assert.hh
@@ -1,8 +1,8 @@
-#ifndef PYTHIA6M_UTIL_ASSERT_LOADED
-#define PYTHIA6M_UTIL_ASSERT_LOADED
+#ifndef PYTHIA6M_CORE_ASSERT_LOADED
+#define PYTHIA6M_CORE_ASSERT_LOADED
 
-#include <pythia6m/util/logger.hh>
-#include <pythia6m/util/exception.hh>
+#include <pythia6m/core/logger.hh>
+#include <pythia6m/core/exception.hh>
 #include <string>
 
 // =============================================================================
diff --git a/pythia6m/util/configuration.cc b/pythia6m/core/configuration.cc
similarity index 99%
rename from pythia6m/util/configuration.cc
rename to pythia6m/core/configuration.cc
index 191ea2ee18a3a6657efea9786fe424afbc8fb29b..48cb1d16a2fc7808caefee54d902d1c258778b79 100644
--- a/pythia6m/util/configuration.cc
+++ b/pythia6m/core/configuration.cc
@@ -1,6 +1,6 @@
 #include "configuration.hh"
 
-#include <pythia6m/util/logger.hh>
+#include <pythia6m/core/logger.hh>
 
 using boost::property_tree::ptree_bad_path;
 using boost::property_tree::ptree_error;
diff --git a/pythia6m/util/configuration.hh b/pythia6m/core/configuration.hh
similarity index 98%
rename from pythia6m/util/configuration.hh
rename to pythia6m/core/configuration.hh
index c519851c5e4bb04437d09a71defa7282cb49817e..b2e9b3d6784fcda61171c4f68c67adada237eaa8 100644
--- a/pythia6m/util/configuration.hh
+++ b/pythia6m/core/configuration.hh
@@ -1,5 +1,5 @@
-#ifndef PYTHIA6M_UTIL_CONFIGURATION_LOADED
-#define PYTHIA6M_UTIL_CONFIGURATION_LOADED
+#ifndef PYTHIA6M_CORE_CONFIGURATION_LOADED
+#define PYTHIA6M_CORE_CONFIGURATION_LOADED
 
 #include <boost/lexical_cast.hpp>
 #include <boost/optional.hpp>
@@ -7,9 +7,9 @@
 #include <boost/property_tree/json_parser.hpp>
 #include <boost/property_tree/ptree.hpp>
 #include <map>
-#include <pythia6m/util/exception.hh>
-#include <pythia6m/util/interval.hh>
-#include <pythia6m/util/stringify.hh>
+#include <pythia6m/core/exception.hh>
+#include <pythia6m/core/interval.hh>
+#include <pythia6m/core/stringify.hh>
 #include <string>
 
 namespace pythia6m {
diff --git a/pythia6m/util/exception.hh b/pythia6m/core/exception.hh
similarity index 85%
rename from pythia6m/util/exception.hh
rename to pythia6m/core/exception.hh
index 04baf292248d697e1e142d697c57654646a1947d..8455573f71c9a560d5031ac8f8700a48d6022f8e 100644
--- a/pythia6m/util/exception.hh
+++ b/pythia6m/core/exception.hh
@@ -1,5 +1,5 @@
-#ifndef PYTHIA6M_UTIL_EXCEPTION_LOADED
-#define PYTHIA6M_UTIL_EXCEPTION_LOADED
+#ifndef PYTHIA6M_CORE_EXCEPTION_LOADED
+#define PYTHIA6M_CORE_EXCEPTION_LOADED
 
 #include <exception>
 #include <string>
diff --git a/pythia6m/util/framework.cc b/pythia6m/core/framework.cc
similarity index 95%
rename from pythia6m/util/framework.cc
rename to pythia6m/core/framework.cc
index 057b73c00094438c518d8eab5b1a27f6d66b1933..b829346b5ddd6d5905eabf364e1394eac457e022 100644
--- a/pythia6m/util/framework.cc
+++ b/pythia6m/core/framework.cc
@@ -3,10 +3,10 @@
 #include <exception>
 #include <cstdlib>
 
-#include <pythia6m/util/exception.hh>
-#include <pythia6m/util/configuration.hh>
-#include <pythia6m/util/logger.hh>
-#include <pythia6m/util/stringify.hh>
+#include <pythia6m/core/exception.hh>
+#include <pythia6m/core/configuration.hh>
+#include <pythia6m/core/logger.hh>
+#include <pythia6m/core/stringify.hh>
 
 #include <TSystem.h>
 
@@ -44,33 +44,8 @@ framework::framework(int argc, char* argv[],
   LOG_INFO("pythia6m", "Starting pythia6m framework");
   LOG_INFO("pythia6m", "Configuration file: " + args_["conf"].as<std::string>());
 
-  // get our run info and number of events to be generated
-  const int run = get_option<int>("run");
-  const int events = get_option<int>("events");
-
   // output file name
-
-  // path
-  output_ = args_["out"].as<std::string>();
-  if (output_.back() != '/') {
-    output_ += '/';
-  }
-
-  // what MC program are we running?
-  output_ += conf_.get<std::string>("type");
-
-  // add optional tag
-  auto tag = conf_.get_optional<std::string>("tag");
-  if (tag && tag->size()) {
-    output_ += "." + *tag;
-  }
-  // add the run info and number of generated events
-  char info[1024];
-  sprintf(info, ".run%05i-%i", run, events);
-  output_ += info;
-
-  // communicate file name to user
-  LOG_INFO("pythia6m", "Output files will be written to: " + output_ + ".*");
+  set_output();
 
   // write the a copy of the configuration file to the output
   ptree settings;
@@ -182,6 +157,34 @@ ptree framework::get_settings() const {
   }
   return settings;
 }
+// =============================================================================
+// Implementation: framework::set_output
+// =============================================================================
+void framework::set_output() {
+  // get our run info and number of events to be generated
+  const int run = get_option<int>("run");
+  const int events = get_option<int>("events");
+
+  output_ = args_["out"].as<std::string>();
+  if (output_.back() != '/') {
+    output_ += '/';
+  }
+
+  // what MC program are we running?
+  output_ += conf_.get<std::string>("type");
+
+  // add optional tag
+  auto tag = conf_.get_optional<std::string>("tag");
+  if (tag && tag->size()) {
+    output_ += "." + *tag;
+  }
+  // add the run info and number of generated events
+  char info[1024];
+  sprintf(info, ".run%05i-%i", run, events);
+  output_ += info;
+  LOG_INFO("pythia6m", "Output files will be written to: " + output_ + ".*");
+}
+
 // =============================================================================
 // Implementation: framework::root_suppress_signals
 // Suppress the ROOT signal handlers, as they can cause undefined behavior and
diff --git a/pythia6m/util/framework.hh b/pythia6m/core/framework.hh
similarity index 92%
rename from pythia6m/util/framework.hh
rename to pythia6m/core/framework.hh
index a832b23a7b91aca7904927f71fe0f47a41a40037..fa3e48d8c9313fc1d3eb5afe944420c3c5891409 100644
--- a/pythia6m/util/framework.hh
+++ b/pythia6m/core/framework.hh
@@ -5,7 +5,7 @@
 #include <vector>
 #include <string>
 
-#include <pythia6m/util/configuration.hh>
+#include <pythia6m/core/configuration.hh>
 
 #include <boost/program_options.hpp>
 
@@ -31,8 +31,8 @@ class framework_parse_error;
 #define MAKE_PYTHIA6M_FRAMEWORK(function)                                      \
   int main(int argc, char* argv[]) {                                           \
     try {                                                                      \
-      pythia6m::framework pythia6m{argc, argv, (function)};                    \
-      pythia6m.run();                                                          \
+      ::pythia6m::framework pyfw{argc, argv, (function)};                      \
+      pyfw.run();                                                              \
       return 0;                                                                \
     } catch (...) {                                                            \
       return 1;                                                                \
@@ -53,9 +53,9 @@ public:
 private:
   boost::program_options::variables_map parse_arguments(int argc,
                                                         char* argv[]) const;
-                                                          
-  // get our settings ptree used to initialize the configuration
-  ptree get_settings() const;
+
+  ptree get_settings() const; // get settings ptree to initialize the conf
+  void set_output();          // set the output file name
 
   // suppress the ROOT signal handler
   int root_suppress_signals() const;
@@ -69,7 +69,6 @@ private:
   configuration conf_;
   std::string output_;
   pythia6m_function_type pythia6m_function_;
-
 };
 } // ns pythia6m
 
diff --git a/pythia6m/util/interval.hh b/pythia6m/core/interval.hh
similarity index 100%
rename from pythia6m/util/interval.hh
rename to pythia6m/core/interval.hh
diff --git a/pythia6m/util/logger.cc b/pythia6m/core/logger.cc
similarity index 100%
rename from pythia6m/util/logger.cc
rename to pythia6m/core/logger.cc
diff --git a/pythia6m/util/logger.hh b/pythia6m/core/logger.hh
similarity index 97%
rename from pythia6m/util/logger.hh
rename to pythia6m/core/logger.hh
index a05d4887a3c13e30bc61524ce20bccd30969b7f6..83a7d92c9eebfcb292a600143ad9acd509045e0a 100644
--- a/pythia6m/util/logger.hh
+++ b/pythia6m/core/logger.hh
@@ -1,5 +1,5 @@
-#ifndef PYTHIA6M_UTIL_LOGGER_LOADED
-#define PYTHIA6M_UTIL_LOGGER_LOADED
+#ifndef PYTHIA6M_CORE_LOGGER_LOADED
+#define PYTHIA6M_CORE_LOGGER_LOADED
 
 #include <ctime>
 #include <ostream>
@@ -85,7 +85,7 @@ void log(const std::string& mtitle, const std::string& mtext,
 
 // =============================================================================
 // log_handler class designed for global usage,
-// threading seutil
+// threading secure
 // =============================================================================
 namespace pythia6m {
 class log_handler {
@@ -127,7 +127,7 @@ private:
 
 // =============================================================================
 // global logger function
-// threading seutil
+// threading secure
 // =============================================================================
 namespace pythia6m {
 template <log_level level>
diff --git a/pythia6m/core/particle.hh b/pythia6m/core/particle.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b537cb305a1cdf60ebc865295c328f2eefb34565
--- /dev/null
+++ b/pythia6m/core/particle.hh
@@ -0,0 +1,151 @@
+#ifndef PYTHIA6M_CORE_EVENT_LOADED
+#define PYTHIA6M_CORE_EVENT_LOADED
+
+#include <TLorentzVector.h>
+#include <cmath>
+#include <functional>
+#include <ostream>
+#include <pythia6m/pythia/pythia6.h>
+#include <vector>
+#include <cstring>
+
+// pythia 6 particle data, created from the pythia6 common blocks
+namespace pythia6m {
+
+// single particle info
+struct particle {
+  int32_t index;
+  int32_t status;
+  int32_t type;
+  int32_t parent;
+  int32_t daughter1; // this is the first daughter, or zero
+  int32_t daughter2; // this is the last daughter, or zero
+  int32_t charge;
+  TLorentzVector vertex;
+  TLorentzVector mom;
+  double mass;
+  double z;
+  double Phperp;
+  double t;
+
+  particle() {}
+  particle(const int32_t idx)
+      : index{idx}
+      , status{glb_pyJets.K[0][idx]}
+      , type{glb_pyJets.K[1][idx]}
+      , parent{glb_pyJets.K[2][idx]}
+      , daughter1{glb_pyJets.K[3][idx]}
+      , daughter1{glb_pyJets.K[4][idx]}
+      , charge{PYCHGE(type) / 3}
+      , vertex{glb_pyJets.V[0][idx] / 10, glb_pyJets.V[1][idx] / 10,
+               glb_pyJets.V[2][idx] / 10, glb_pyJets.V[3][idx] / 10}
+      , mom{glb_pyJets.P[0][idx], glb_pyJets.P[1][idx], glb_pyJets.P[2][idx],
+            glb_pyJets.P[3][idx]}
+      , mass{glb_pyJets.P[4][idx]} {}
+
+  // initial state particle?
+  constexpr bool is_init() const { return status == 21; }
+  // lund (final state) particle?
+  constexpr bool is_lund() const {
+    return status >= 1 && status <= 10 && status != 2 && status != 3;
+  }
+  // virtual particle?
+  constexpr bool is_virt() const { return !is_init() && !is_lund(); }
+
+  // print a LUND particle record for this particle
+  void print(std::ostream& os) {
+    char buf[1024];
+    snprintf(
+        buf, 1024, "%4i %4i %4i %4i %4i %4i %12.6e %12.6e %12.6e %12.6e %12.6e "
+                   "%12.6e %12.6e %12.6e\n",
+        index, charge, 1 /* active */, type, parent, daughter1, mom.Px(),
+        mom.Py(), mom.Pz(), mom.E(), mass, vertex.X(), vertex.Y(), vertex.Z());
+    os << buf;
+  }
+};
+
+// main event class, mainly a wrapper around a vector<particle>
+struct event {
+  int32_t evnum;
+  int32_t evgen;
+  double xsec;
+  int32_t process;
+  std::vector<particle> particles;
+
+  event()
+      : evnum{glb_event.evNum}
+      , evgen{glb_event.nEvGen}
+      , xsec{glb_pyPars.pari[0]}
+      , process{glb_event.process} {
+    for (int idx = 0; idx < glb_pyJets.n; ++idx) {
+      particles.push_back(idx);
+    }
+  }
+
+  constexpr const particle& beam() const {
+    return particles[0];
+  }
+  constexpr const particle& target() const {
+    return particles[1];
+  }
+
+  // print a LUND event record for all final state particles
+  void print(std::ostream& os(),
+             std::function<bool(const particle&)> selector) const {
+    // check what particles we want to write
+    std::vector<int> selected; // indices of the selected particles
+    for (const auto& part : particles) {
+      if (part.is_lund() && selector(part)) {
+        selected.push_back(part.idx);
+      }
+    }
+    // output buffer
+    char buf[1024];
+    // write the header (only number of events needed)
+    snprintf(buf, 1024,
+             "%4i %4i %4i %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e\n",
+             selected.size(), 0, 0, 0, 0, 0, 0, 0, 0, 0);
+    os << buf;
+    // write the selected particle lines
+    for (int i : selected) {
+      particles[i].print(os);
+    }
+  }
+};
+
+// e+p leptoproduction 
+struct ep_event : event {
+  TLorentzVector q;
+  double s;
+  double Q2;
+  double nu;
+  double x;
+  double y;
+  double W2;
+
+  ep_event() : event(),
+  q {beam.mom-particles[glb_event.scatIdx].mom}
+  
+  
+  {}
+
+  constexpr const particle& scat() const {
+    return particles[glb_event.scatIdx];
+  }
+
+};
+struct ep_event {
+  TLorentzVector beam;
+  TLorentzVector target;
+  TLorentzVector scat;
+  TLorentzVector q;
+
+// gamma+p photoproduction
+struct photo_event
+
+
+
+};
+
+
+#endif
diff --git a/pythia6m/util/progress_meter.hh b/pythia6m/core/progress_meter.hh
similarity index 93%
rename from pythia6m/util/progress_meter.hh
rename to pythia6m/core/progress_meter.hh
index f8b9e8917377f7e668ffbfc1605226f0464af282..8f5822ec61ec98d75d81a189e5b5f2ce793b6458 100644
--- a/pythia6m/util/progress_meter.hh
+++ b/pythia6m/core/progress_meter.hh
@@ -1,5 +1,5 @@
-#ifndef PYTHIA6M_UTIL_PROGRESS_METER_LOADED
-#define PYTHIA6M_UTIL_PROGRESS_METER_LOADED
+#ifndef PYTHIA6M_CORE_PROGRESS_METER_LOADED
+#define PYTHIA6M_CORE_PROGRESS_METER_LOADED
 
 #include <cmath>
 #include <cstdint>
diff --git a/pythia6m/util/stringify.hh b/pythia6m/core/stringify.hh
similarity index 95%
rename from pythia6m/util/stringify.hh
rename to pythia6m/core/stringify.hh
index f2def1aab7ab6b3b927abe1514794f7fc27a5fe4..36c9c04b0592759723c4bb55b2a5ab8e4eceb01d 100644
--- a/pythia6m/util/stringify.hh
+++ b/pythia6m/core/stringify.hh
@@ -1,5 +1,5 @@
-#ifndef PYTHIA6M_UTIL_STRINGIFY_LOADED
-#define PYTHIA6M_UTIL_STRINGIFY_LOADED
+#ifndef PYTHIA6M_CORE_STRINGIFY_LOADED
+#define PYTHIA6M_CORE_STRINGIFY_LOADED
 
 #include <string>
 #include <utility>
@@ -7,7 +7,7 @@
 
 #include <boost/lexical_cast.hpp>
 
-#include <pythia6m/util/type_traits.hh>
+#include <pythia6m/core/type_traits.hh>
 
 namespace pythia6m {
 
diff --git a/pythia6m/util/type_traits.hh b/pythia6m/core/type_traits.hh
similarity index 97%
rename from pythia6m/util/type_traits.hh
rename to pythia6m/core/type_traits.hh
index 71bac487a556947e7417669d2d0ba541e777fbbc..c7bd54a2a3b1375be5a9e251f6029164ee0921f5 100644
--- a/pythia6m/util/type_traits.hh
+++ b/pythia6m/core/type_traits.hh
@@ -1,5 +1,5 @@
-#ifndef PYTHIA6M_UTIL_TYPE_TRAITS_LOADED
-#define PYTHIA6M_UTIL_TYPE_TRAITS_LOADED
+#ifndef PYTHIA6M_CORE_TYPE_TRAITS_LOADED
+#define PYTHIA6M_CORE_TYPE_TRAITS_LOADED
 
 #include <type_traits>
 
diff --git a/pythia6m/interface/CMakeLists.txt b/pythia6m/interface/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..43d892f7c6361b51d466a09006ca089eff89e173
--- /dev/null
+++ b/pythia6m/interface/CMakeLists.txt
@@ -0,0 +1,47 @@
+################################################################################
+## CMAKE Settings
+################################################################################
+set (LIBRARY "pythia6m_interface")
+set (TARGETS ${TARGETS} ${LIBRARY} PARENT_SCOPE)
+
+################################################################################
+## Sources and install headers
+################################################################################
+file (GLOB SOURCES "*.cc")
+file (GLOB HEADERS "*.hh")
+
+################################################################################
+## Include directories
+################################################################################
+include_directories("${PROJECT_SOURCE_DIR}")
+include_directories("${CMAKE_CURRENT_SOURCE_DIR}")
+
+################################################################################
+## External Libraries
+################################################################################
+## ROOT
+find_package(ROOT REQUIRED)
+include_directories(${ROOT_INCLUDE_DIR})
+## boost
+find_package(Boost COMPONENTS program_options filesystem system REQUIRED)
+include_directories(${Boost_INCLUDE_DIRS})
+
+################################################################################
+## Compile and Link
+################################################################################
+add_library(${LIBRARY} ${SOURCES})
+#target_link_libraries(${LIBRARY} ${ROOT_LIBRARIES} ${Boost_LIBRARIES})
+set_target_properties(${LIBRARY} PROPERTIES 
+  VERSION ${PYTHIA6M_VERSION} 
+  SOVERSION ${PYTHIA6M_SOVERSION}
+  PUBLIC_HEADER "${HEADERS}")
+
+################################################################################
+## Export and Install
+################################################################################
+install(TARGETS ${LIBRARY}
+  EXPORT ${PROJECT_NAME}-targets
+  RUNTIME DESTINATION "${INSTALL_BIN_DIR}" COMPONENT bin
+  LIBRARY DESTINATION "${INSTALL_LIB_DIR}" COMPONENT lib
+  ARCHIVE DESTINATION "${INSTALL_LIB_DIR}" COMPONENT lib
+  PUBLIC_HEADER DESTINATION "${INSTALL_INCLUDE_DIR}/interface" COMPONENT dev)
diff --git a/pythia6m/interface/event.cc b/pythia6m/interface/event.cc
new file mode 100644
index 0000000000000000000000000000000000000000..217ce7689864033d9bf78b7595b29d9f2f91d642
--- /dev/null
+++ b/pythia6m/interface/event.cc
@@ -0,0 +1,125 @@
+#include "event.hh"
+#include <cstring>
+#include <iostream>
+#include <pythia6m/gmc/gen_set.h>
+#include <pythia6m/gmc/gmc.h>
+#include <pythia6m/gmc/mc_set.h>
+#include <pythia6m/gmc/mcevnt.h>
+#include <pythia6m/pythia6/pythia6.h>
+
+namespace pythia6m {
+
+particle::particle(const int32_t idx)
+    : index{idx} //
+    , status{glb_pyJets.K[0][idx]}
+    , type{glb_pyJets.K[1][idx]}
+    , parent{glb_pyJets.K[2][idx]}
+    , daughter1{glb_pyJets.K[3][idx]}
+    , daughter2{glb_pyJets.K[4][idx]}
+    , charge{PYCHGE(type) / 3}
+    , vertex{glb_pyJets.V[0][idx] / 10, glb_pyJets.V[1][idx] / 10,
+             glb_pyJets.V[2][idx] / 10, glb_pyJets.V[3][idx] / 10}
+    , mom{glb_pyJets.P[0][idx], glb_pyJets.P[1][idx], glb_pyJets.P[2][idx],
+          glb_pyJets.P[3][idx]}
+    , mass{glb_pyJets.P[4][idx]} {}
+void particle::print(std::ostream& os) const {
+  // write a LUND particle line. index is stored as the Fortran index
+  // (consistent with parent/daughter indices)
+  char buf[1024];
+  snprintf(
+      buf, 1024, "%4i %4i %4i %4i %4i %4i %12.6e %12.6e %12.6e %12.6e %12.6e "
+                 "%12.6e %12.6e %12.6e\n",
+      index + 1, charge, is_lund() /*only LUND particles labeled as active*/,
+      type, parent, daughter1, mom.Px(), mom.Py(), mom.Pz(), mom.E(), mass,
+      vertex.X(), vertex.Y(), vertex.Z());
+  os << buf;
+}
+
+event::event()
+    : evnum{glb_event.evNum}
+    , evgen{glb_event.nEvGen}
+    , xsec{glb_pyPars.pari[0]}
+    , process{glb_event.process}
+    , s{(particles[0].mom + particles[1].mom).M2()} {
+  for (int idx = 0; idx < glb_pyJets.n; ++idx) {
+    particles.push_back(idx);
+  }
+}
+event::event(event&& rhs)
+    : evnum{rhs.evnum}
+    , evgen{rhs.evgen}
+    , xsec{rhs.xsec}
+    , process{rhs.process}
+    , s{rhs.s}
+    , particles{std::move(rhs.particles)} {}
+
+event& event::operator=(event&& rhs) {
+  evnum = rhs.evnum;
+  evgen = rhs.evgen;
+  xsec = rhs.xsec;
+  process = rhs.process;
+  s = rhs.s;
+  particles = std::move(rhs.particles);
+  return *this;
+}
+void event::print(std::ostream& os,
+                  std::function<bool(const particle&)> selector) const {
+  // check what particles we want to write
+  std::vector<int> selected; // indices of the selected particles
+  for (const auto& part : particles) {
+    if (part.is_lund() && selector(part)) {
+      selected.push_back(part.index);
+    }
+  }
+  // output buffer
+  char buf[1024];
+  // write the header (only number of events needed)
+  snprintf(buf, 1024,
+           "%4zu %4i %4i %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e\n",
+           selected.size(), 0, 0, 0., 0., 0., 0., 0., 0., 0.);
+  os << buf;
+  // write the selected particle lines
+  for (int i : selected) {
+    particles[i].print(os);
+  }
+}
+
+ep_event::ep_event()
+    : event()
+    , q{beam().mom - particles[glb_event.scatIdx].mom}
+    , Q2{glb_event.Q2Gen}
+    , nu{glb_event.nuGen}
+    , x{glb_event.xGen}
+    , y{glb_event.yGen}
+    , W2{glb_event.W2Gen}
+    , scat_index{glb_event.scatIdx} {}
+ep_event::ep_event(ep_event&& rhs)
+    : event{std::move(rhs)}
+    , q{rhs.q}
+    , Q2{rhs.Q2}
+    , nu{rhs.nu}
+    , x{rhs.x}
+    , y{rhs.y}
+    , W2{rhs.W2}
+    , scat_index{glb_event.scatIdx} {}
+ep_event& ep_event::operator=(ep_event&& rhs) {
+  q = rhs.q;
+  Q2 = rhs.Q2;
+  nu = rhs.nu;
+  x = rhs.x;
+  y = rhs.y;
+  W2 = rhs.W2;
+  scat_index = rhs.scat_index;
+  static_cast<event&>(*this) = std::move(rhs);
+  return *this;
+}
+
+photo_event::photo_event() : event() {}
+
+photo_event::photo_event(photo_event&& rhs) : event(std::move(rhs)) {}
+photo_event& photo_event::operator=(photo_event&& rhs) {
+  static_cast<event&>(*this) = std::move(rhs);
+  return *this;
+}
+
+} // ns pythia6m
diff --git a/pythia6m/interface/event.hh b/pythia6m/interface/event.hh
new file mode 100644
index 0000000000000000000000000000000000000000..57784f6eeb8e5e6433f6177b8848d52b76e7f706
--- /dev/null
+++ b/pythia6m/interface/event.hh
@@ -0,0 +1,111 @@
+#ifndef PYTHIA6M_CORE_EVENT_LOADED
+#define PYTHIA6M_CORE_EVENT_LOADED
+
+#include <Math/Vector4D.h>
+#include <cmath>
+#include <functional>
+#include <ostream>
+#include <vector>
+
+// pythia 6 particle data, created from the pythia6 common blocks
+namespace pythia6m {
+
+// single particle info
+struct particle {
+  int32_t index; // Fortran index (starts at 1)
+  int32_t status;
+  int32_t type;
+  int32_t parent;
+  int32_t daughter1; // this is the first daughter, or zero
+  int32_t daughter2; // this is the last daughter, or zero
+  int32_t charge;
+  ROOT::Math::XYZTVector vertex;
+  ROOT::Math::XYZTVector mom;
+  double mass;
+
+  particle() = default;
+  
+  // create a particle from using the global pythia particle index
+  particle(const int32_t idx);
+
+  // initial state particle?
+  constexpr bool is_init() const { return status == 21; }
+  // lund (final state) particle?
+  constexpr bool is_lund() const {
+    return status >= 1 && status <= 10 && status != 2 && status != 3;
+  }
+  // virtual particle?
+  constexpr bool is_virt() const { return !is_init() && !is_lund(); }
+
+  // print a LUND particle record for this particle
+  void print(std::ostream& os) const;
+};
+
+// main event class, mainly a wrapper around a vector<particle>
+struct event {
+  int32_t evnum;
+  int32_t evgen;
+  double xsec;
+  int32_t process;
+  double s;
+  std::vector<particle> particles;
+
+  // create event using the pythia global event info
+  event();
+
+  // copy constructors
+  event(const event&) = default;
+  event(event&& rhs);
+
+  // assignment
+  event& operator=(const event&) = default;
+  event& operator=(event&& rhs);
+
+  constexpr const particle& beam() const {
+    return particles[0];
+  }
+  constexpr const particle& target() const {
+    return particles[1];
+  }
+
+  // print a LUND event record for all final state particles
+  void print(std::ostream& os,
+             std::function<bool(const particle&)> selector) const;
+};
+
+// e+p leptoproduction 
+struct ep_event : event {
+  ROOT::Math::XYZTVector q;
+  double Q2;
+  double nu;
+  double x;
+  double y;
+  double W2;
+  int32_t scat_index;
+
+  ep_event();
+  ep_event(const ep_event&) = default;
+  ep_event(ep_event&& rhs);
+  ep_event& operator=(const ep_event&) = default;
+  ep_event& operator=(ep_event&& rhs);
+
+  const particle& scat() const { return particles[scat_index]; }
+};
+// gamma+p photoproduction
+struct photo_event : event {
+
+  photo_event();
+  photo_event(const photo_event&) = default;
+  photo_event(photo_event&& rhs);
+  photo_event& operator=(const photo_event&) = default;
+  photo_event& operator=(photo_event&& rhs);
+
+  double t(int32_t index) const {
+    return (beam().mom - particles[index].mom).M2();
+  }
+};
+
+
+} // ns pythia6m
+
+#endif
diff --git a/pythia6m/interface/event_out.hh b/pythia6m/interface/event_out.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9c6b956b83d7fa251a12f95c018cbe94ec94e468
--- /dev/null
+++ b/pythia6m/interface/event_out.hh
@@ -0,0 +1,170 @@
+#include <Math/Vector4D.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <memory>
+#include <pythia6m/core/assert.hh>
+#include <pythia6m/core/logger.hh>
+#include <pythia6m/interface/event.hh>
+#include <string>
+#include <vector>
+
+namespace pythia6m {
+
+// store event to ROOT file
+class event_out {
+public:
+  event_out(std::shared_ptr<TFile> f, const std::string& name) : file_{f} {
+    LOG_INFO("event_out", "Initializing ROOT output stream");
+    tassert(file_, "invalid file pointer");
+    file_->cd();
+    tree_ = new TTree(name.c_str(), name.c_str());
+    tassert(tree_, "Failed to inialize tree " + name);
+    create_branches();
+  }
+  ~event_out() { tree_->AutoSave(); }
+
+  void push(const event& e) {
+    clear();
+    evnum = e.evnum;
+    evgen = e.evgen;
+    xsec = e.xsec;
+    process = e.process;
+    s = e.s;
+    for (const auto& part : e.particles) {
+      add(part);
+    }
+    tree_->Fill();
+  }
+
+protected:
+  void clear() {
+    n_part = 0;
+    index.clear();
+    status.clear();
+    type.clear();
+    parent.clear();
+    daughter1.clear();
+    daughter2.clear();
+    charge.clear();
+    vertex.clear();
+    mom.clear();
+    mass.clear();
+    init.clear();
+    lund.clear();
+  }
+  void add(const particle& part) {
+    n_part += 1;
+    index.push_back(part.index + 1);
+    status.push_back(part.status);
+    type.push_back(part.type);
+    parent.push_back(part.parent);
+    daughter1.push_back(part.daughter1);
+    daughter2.push_back(part.daughter2);
+    charge.push_back(part.charge);
+    vertex.push_back(part.vertex);
+    mom.push_back(part.mom);
+    mass.push_back(part.mass);
+    init.push_back(part.is_init());
+    lund.push_back(part.is_lund());
+  }
+
+  std::shared_ptr<TFile> file_;
+  TTree* tree_;
+
+private:
+  void create_branches() {
+    tree_->Branch("event", &evnum);
+    tree_->Branch("evgen", &evgen);
+    tree_->Branch("xsec", &xsec);
+    tree_->Branch("process", &process);
+    tree_->Branch("s", &s);
+    tree_->Branch("n_part", &n_part);
+    tree_->Branch("index", &index);
+    tree_->Branch("status", &status);
+    tree_->Branch("type", &type);
+    tree_->Branch("parent", &daughter1);
+    tree_->Branch("daughter1", &daughter1);
+    tree_->Branch("daughter2", &daughter2);
+    tree_->Branch("charge", &charge);
+    tree_->Branch("vertex", &vertex);
+    tree_->Branch("mom", &mom);
+    tree_->Branch("mass", &mass);
+    tree_->Branch("init", &init);
+    tree_->Branch("lund", &lund);
+  }
+
+  // event data
+  int32_t evnum;
+  int32_t evgen;
+  double xsec;
+  int32_t process;
+  double s;
+
+  // particle data
+  int32_t n_part;
+  std::vector<int32_t> index;
+  std::vector<int32_t> status;
+  std::vector<int32_t> type;
+  std::vector<int32_t> parent;
+  std::vector<int32_t> daughter1; // this is the first daughter, or zero
+  std::vector<int32_t> daughter2; // this is the last daughter, or zero
+  std::vector<int32_t> charge;
+  std::vector<ROOT::Math::XYZTVector> vertex;
+  std::vector<ROOT::Math::XYZTVector> mom;
+  std::vector<double> mass;
+  std::vector<bool> init;
+  std::vector<bool> lund;
+};
+
+// derived class for ep-event
+class ep_event_out : public event_out {
+public:
+  ep_event_out(std::shared_ptr<TFile> f, const std::string& name)
+      : event_out(std::move(f), name) {
+    create_branches();
+  }
+
+  void push(const ep_event& e) {
+    Q2 = e.Q2;
+    nu = e.nu;
+    x = e.x;
+    y = e.y;
+    W2 = e.W2;
+    scat_index = e.scat_index;
+    event_out::push(e);
+  }
+
+private:
+  void create_branches() {
+    tree_->Branch("Q2", &Q2);
+    tree_->Branch("nu", &nu);
+    tree_->Branch("x", &x);
+    tree_->Branch("y", &y);
+    tree_->Branch("W2", &W2);
+    tree_->Branch("scat_index", &scat_index);
+  }
+
+  double Q2;
+  double nu;
+  double x;
+  double y;
+  double W2;
+  int32_t scat_index;
+};
+
+class photo_event_out : public event_out {
+public:
+  photo_event_out(std::shared_ptr<TFile> f, const std::string& name)
+      : event_out(std::move(f), name) {
+    create_branches();
+  }
+
+  void push(const photo_event& e) { event_out::push(e); }
+
+private:
+  void create_branches() {
+    ; // nothing here
+  }
+};
+
+} // ns pythia6m
diff --git a/pythia6m/interface/pythia6m_generator.cc b/pythia6m/interface/pythia6m_generator.cc
new file mode 100644
index 0000000000000000000000000000000000000000..e022c28b36464cf1b838661096812ca8b3f6d45e
--- /dev/null
+++ b/pythia6m/interface/pythia6m_generator.cc
@@ -0,0 +1,181 @@
+#include "pythia6m_generator.hh"
+
+
+#include <cstring>
+#include <pythia6m/core/exception.hh>
+#include <pythia6m/core/logger.hh>
+#include <pythia6m/fmotion/fmotion.h>
+#include <pythia6m/gmc/gen_set.h>
+#include <pythia6m/gmc/gmc.h>
+#include <pythia6m/gmc/mc_set.h>
+#include <pythia6m/gmc/mcevnt.h>
+#include <pythia6m/gmc/random.h>
+#include <pythia6m/interface/event.hh>
+#include <pythia6m/pythia6/pythia6.h>
+
+// globals from GMC
+GeneratorSettings glb_gen;
+MCSettings glb_mc;
+EventData glb_event;
+
+namespace pythia6m {
+pythia6m_generator::pythia6m_generator(const configuration& conf, const string_path& path)
+    : conf_(conf, path) {
+  LOG_INFO("PYTHIA6M", "Initializing the generator");
+  // configuration
+  conf_generator();
+  conf_gmc();
+  // run and events (stored in the top level path of the configuration)
+  glb_mc.RunNo = conf.get<int>("run");
+  glb_mc.NEvents = conf.get<int>("events");
+  LOG_INFO("conf", "Will generate " + std::to_string(glb_mc.NEvents) +
+                       " events for run " + std::to_string(glb_mc.RunNo));
+
+  // initialization
+  init_random();
+  init_event();
+  init_gmc();
+}
+
+pythia6m_generator::~pythia6m_generator() { 
+  GMC_RUNSTAT();
+  LOG_INFO("PYTHIA6M",
+           "TOTAL Generated events: " + std::to_string(glb_event.nEvGen));
+  LOG_INFO("PYTHIA6M", "TOTAL PYTHIA Cross Section [ubarn]: " +
+                           std::to_string(glb_pyPars.pari[0]));
+}
+
+// =============================================================================
+// Pythia6m private: configure generator
+// =============================================================================
+void pythia6m_generator::conf_generator() const {
+  LOG_INFO("conf", "Configuring generator");
+
+  // User control switch MSEL (manual section 9.2, also 8.3.1)
+  // MSEL=2 additionally enables the elastic and diffactive components of the
+  // VMD and GVMD parts
+  glb_gen.PhotoType = 2;
+
+  // Let LUND decay some of the semi-stable particles
+  glb_gen.DecayPiK = 0;
+  glb_gen.DecayLamKs = 1;
+
+  // RADGEN options, currently disabled
+  glb_gen.radgenGenLUT = 0;
+  glb_gen.radgenUseLUT = 0;
+  sprintf(glb_gen.GenRad, "NO\n");
+
+  // activate Fermi momentum? (experimental)
+  glb_gen.enableFM = 1;
+  glb_gen.FMNBins = 1000; // Number of bins for FM discretization
+  glb_gen.FMCutOff = 1.5; // Max FM (in GeV)
+
+  // Parameterizations for F2Pythia and R
+  // (I believe these are only relevant when running with RADGEN)
+  sprintf(glb_gen.FStruct, "F2PY");
+  sprintf(glb_gen.R, "1990");
+
+  // PYTHIA model to run with:
+  //  * RAD for realistic DIS
+  //  * GAM for real photon beam
+  sprintf(glb_gen.PyModel, "RAD");
+}
+// =============================================================================
+// Pythia6m private: configure GMC
+// =============================================================================
+void pythia6m_generator::conf_gmc() const {
+  LOG_INFO("conf", "Configuring MC");
+
+  glb_mc.TarA = conf_.get<int>("target/A");
+  glb_mc.TarA = conf_.get<int>("target/Z");
+  glb_mc.TargetNucleon[0] = conf_.get<char>("target/nucleon");
+  glb_mc.EneBeam = conf_.get<double>("beam/energy");
+  glb_mc.Q2Min = conf_.get<double>("Q2/min");
+  glb_mc.Q2Max = conf_.get<double>("Q2/max");
+  glb_mc.xMin = conf_.get<double>("x/min");
+  glb_mc.xMax = conf_.get<double>("x/max");
+  glb_mc.yMin = conf_.get<double>("y/min");
+  glb_mc.yMax = conf_.get<double>("y/max");
+
+  // hardcoded electron beam for now
+  sprintf(glb_mc.BeamParType, "ELEC");
+  // hardcoded Beam and target Polarization
+  // needed because RADGEN can also run with polarized generators
+  glb_mc.PBValue = 0;
+  glb_mc.PTValue = 0;
+  sprintf(glb_mc.PBeam, "L\n");
+  sprintf(glb_mc.PTarget, "L\n");
+
+  // reconsider fermi motion when running for A=1 options
+  if (glb_mc.TarA == 1) {
+    glb_gen.enableFM = 0;
+  }
+}
+// =============================================================================
+// Pythia6m private: init event data
+// =============================================================================
+void pythia6m_generator::init_event() const {
+  LOG_INFO("init", "Initializing the event common blocks");
+  // beam particle
+  if (!strncmp(glb_mc.BeamParType, "POSI", 4)) {
+    glb_event.beamLType = -11;
+  } else if (!strncmp(glb_mc.BeamParType, "ELEC", 4)) {
+    glb_event.beamLType = 11;
+  } else if (!strncmp(glb_mc.BeamParType, "GAM", 3)) {
+    glb_event.beamLType = 22;
+  } else {
+    throw exception("Unrecognized beam particle. This should NEVER HAPPEN!");
+  }
+
+  // Init FM engine
+  if (glb_gen.enableFM) {
+    INIT_FM(glb_mc.TarZ, glb_mc.TarA, glb_gen.FMCutOff, glb_gen.FMNBins);
+  }
+
+  glb_event.evNum = 0;
+  glb_event.nEvGen = 0;
+  glb_event.nEvGenPrev = 0;
+  glb_event.nErr = 0;
+
+  glb_event.weight = 1.f;
+  glb_event.Q2Gen = 0.f;
+  glb_event.nuGen = 0.f;
+  glb_event.xGen = 0.f;
+  glb_event.yGen = 0.f;
+  glb_event.W2Gen = 0.f;
+  glb_event.thetaGen = 0.f;
+  glb_event.phiGen = 0.f;
+  glb_event.EGen = 0.f;
+  glb_event.pGen = 0.f;
+  glb_event.pxGen = 0.f;
+  glb_event.pyGen = 0.f;
+  glb_event.pzGen = 0.f;
+  glb_event.vxGen = 0.f;
+  glb_event.vyGen = 0.f;
+  glb_event.vzGen = 0.f;
+}
+// =============================================================================
+// Pythia6m private: init GMC
+// =============================================================================
+void pythia6m_generator::init_gmc() const {
+  LOG_INFO("init", "Initializing GMC/PYTHIA");
+  int ok;
+  GMC_INIT(ok);
+  if (!ok) {
+    throw exception("Failed to initialized GMC/PYTHIA");
+  }
+}
+// =============================================================================
+// Pythia6m private: init RNG
+// =============================================================================
+void pythia6m_generator::init_random() const {
+  LOG_INFO("init",
+           "initializing the random generator with seed (run number): " +
+               std::to_string(glb_mc.RunNo));
+  int dummy1{0};
+  int dummy2{0};
+  const char* command{" "};
+  RNDMQ(dummy1, dummy2, glb_mc.RunNo, command);
+}
+
+} // ns pythia6m
diff --git a/pythia6m/interface/pythia6m_generator.hh b/pythia6m/interface/pythia6m_generator.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0a9a67b68689893843578382afb4b144452f25ff
--- /dev/null
+++ b/pythia6m/interface/pythia6m_generator.hh
@@ -0,0 +1,81 @@
+#ifndef PYTHIA6M_CORE_PYTHIA6M_GENERATOR_LOADED
+#define PYTHIA6M_CORE_PYTHIA6M_GENERATOR_LOADED
+
+#include <cstring>
+#include <pythia6m/core/configuration.hh>
+#include <pythia6m/fmotion/fmotion.h>
+#include <pythia6m/gmc/gen_set.h>
+#include <pythia6m/gmc/gmc.h>
+#include <pythia6m/gmc/mc_set.h>
+#include <pythia6m/gmc/mcevnt.h>
+#include <pythia6m/gmc/random.h>
+#include <pythia6m/interface/event.hh>
+#include <pythia6m/pythia6/pythia6.h>
+
+namespace pythia6m {
+
+// =============================================================================
+// pythia6m_generator MC master class
+// =============================================================================
+class pythia6m_generator {
+public:
+  pythia6m_generator(const configuration& conf, const string_path& path);
+  ~pythia6m_generator();
+
+  template <class Event> Event generate() const;
+  configuration& conf() { return conf_; }
+  const configuration& conf() const { return conf_; }
+
+private:
+  void conf_generator() const;
+  void conf_gmc() const;
+  void init_event() const;
+  void init_gmc() const;
+  void init_random() const;
+  // void init_radgen(); // TODO
+
+  void clear_event_kinematics() const;
+
+  configuration conf_;
+};
+
+// =============================================================================
+// implementation/pythia6m_generator: generate a new good event
+// =============================================================================
+template <class Event> Event pythia6m_generator::generate() const {
+  // start out clean
+  clear_event_kinematics();
+  int nEvGen{0};
+
+  // generate a good event
+  GMC_EVENT(nEvGen);
+
+  // update counts
+  glb_event.nEvGen += nEvGen;
+  glb_event.evNum += 1;
+
+  // get all event info (created from the globals, event type specified by
+  // user)
+  Event ev;
+  return ev;
+}
+// =============================================================================
+// implementation/pythia6m_generator private: clear_event_kinematics
+// =============================================================================
+inline void pythia6m_generator::clear_event_kinematics() const {
+  glb_event.scatIdx = -1;
+  glb_event.Q2Gen = 0.f;
+  glb_event.nuGen = 0.f;
+  glb_event.xGen = 0.f;
+  glb_event.yGen = 0.f;
+  glb_event.W2Gen = 0.f;
+  glb_event.pxGen = 0.f;
+  glb_event.pyGen = 0.f;
+  glb_event.pzGen = 0.f;
+  glb_event.vxGen = 0.f;
+  glb_event.vyGen = 0.f;
+  glb_event.vzGen = 0.f;
+}
+
+} // ns pythia6m
+#endif