initial commit

29 files changed, 3757 insertions, 0 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000..46aca8b
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,65 @@
+#----------------------------------------------------------------------------
+# Setup the project
+cmake_minimum_required(VERSION 2.6 FATAL_ERROR)
+project(G4BeamTest)
+
+#----------------------------------------------------------------------------
+# Find Geant4 package, activating all available UI and Vis drivers by default
+# You can set WITH_GEANT4_UIVIS to OFF via the command line or ccmake/cmake-gui
+# to build a batch mode only executable
+#
+option(WITH_GEANT4_UIVIS "Build example with Geant4 UI and Vis drivers" ON)
+if(WITH_GEANT4_UIVIS)
+  find_package(Geant4 REQUIRED ui_all vis_all)
+else()
+  find_package(Geant4 REQUIRED)
+endif()
+
+#----------------------------------------------------------------------------
+# Setup Geant4 include directories and compile definitions
+#
+include(${Geant4_USE_FILE})
+
+#----------------------------------------------------------------------------
+# Locate sources and headers for this project
+#
+include_directories(${PROJECT_SOURCE_DIR}/include 
+                    ${Geant4_INCLUDE_DIR})
+file(GLOB sources ${PROJECT_SOURCE_DIR}/src/*.cxx)
+file(GLOB headers ${PROJECT_SOURCE_DIR}/include/*.h)
+
+#----------------------------------------------------------------------------
+# Add the executable, and link it to the Geant4 libraries
+#
+add_executable(G4BeamTest G4BeamTest.cc ${sources} ${headers})
+target_link_libraries(G4BeamTest ${Geant4_LIBRARIES} )
+
+#----------------------------------------------------------------------------
+# Copy all scripts to the build directory, i.e. the directory in which we
+# build G4BeamTest. This is so that we can run the executable directly because it
+# relies on these scripts being in the current working directory.
+#
+set(G4BeamTest_SCRIPTS
+    G4BeamTest.out
+    G4BeamTest.in
+    optPhoton.mac
+    gui.mac
+    icons.mac
+    run.png
+    vis.mac
+  )
+
+foreach(_script ${G4BeamTest_SCRIPTS})
+  configure_file(
+    ${PROJECT_SOURCE_DIR}/${_script}
+    ${PROJECT_BINARY_DIR}/${_script}
+    COPYONLY
+    )
+endforeach()
+
+#----------------------------------------------------------------------------
+# Install the executable to 'bin' directory under CMAKE_INSTALL_PREFIX
+#
+install(TARGETS G4BeamTest DESTINATION bin)
+
+
diff --git a/G4BeamTest.cxx b/G4BeamTest.cxx
new file mode 100644
index 0000000..afe05dd
--- /dev/null
+++ b/G4BeamTest.cxx
@@ -0,0 +1,141 @@
+#ifdef G4MULTITHREADED
+#include "G4MTRunManager.hh"
+#else
+#include "G4RunManager.hh"
+#endif
+
+#include "G4UImanager.hh"
+
+#include "BtPhysicsList.hh"
+#include "BtDetectorConstruction.hh"
+#include "BtActionInitialization.hh"
+
+#ifdef G4VIS_USE
+#include "G4VisExecutive.hh"
+#endif
+
+#ifdef G4UI_USE
+#include "G4UIExecutive.hh"
+#endif
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+namespace {
+  void PrintUsage() {
+    G4cerr << " Usage: " << G4endl;
+    G4cerr << " OpNovice [-m macro ] [-u UIsession] [-t nThreads] [-r seed] "
+           << G4endl;
+    G4cerr << "   note: -t option is available only for multi-threaded mode."
+           << G4endl;
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+int main(int argc,char** argv)
+{
+  // Evaluate arguments
+  //
+  if ( argc > 9 ) {
+    PrintUsage();
+    return 1;
+  }
+
+  G4String macro;
+  G4String session;
+#ifdef G4MULTITHREADED
+  G4int nThreads = 0;
+#endif
+
+  G4long myseed = 345354;
+  for ( G4int i=1; i<argc; i=i+2 ) {
+     if      ( G4String(argv[i]) == "-m" ) macro   = argv[i+1];
+     else if ( G4String(argv[i]) == "-u" ) session = argv[i+1];
+     else if ( G4String(argv[i]) == "-r" ) myseed  = atoi(argv[i+1]);
+#ifdef G4MULTITHREADED
+     else if ( G4String(argv[i]) == "-t" ) {
+                    nThreads = G4UIcommand::ConvertToInt(argv[i+1]);
+    }
+#endif
+    else {
+      PrintUsage();
+      return 1;
+    }
+  }
+
+  // Choose the Random engine
+  //
+  G4Random::setTheEngine(new CLHEP::RanecuEngine);
+
+  // Construct the default run manager
+  //
+#ifdef G4MULTITHREADED
+  G4MTRunManager * runManager = new G4MTRunManager;
+  if ( nThreads > 0 ) runManager->SetNumberOfThreads(nThreads);
+#else
+  G4RunManager * runManager = new G4RunManager;
+#endif
+
+  // Seed the random number generator manually
+  G4Random::setTheSeed(myseed);
+
+  // Set mandatory initialization classes
+  //
+  // Detector construction
+  runManager-> SetUserInitialization(new btDetectorConstruction());
+  // Physics list
+  runManager-> SetUserInitialization(new btPhysicsList());
+  // User action initialization
+  runManager->SetUserInitialization(new btActionInitialization());
+
+  // Initialize G4 kernel
+  //
+  runManager->Initialize();
+
+#ifdef G4VIS_USE
+  // Initialize visualization
+  //
+  G4VisManager* visManager = new G4VisExecutive;
+  // G4VisExecutive can take a verbosity argument - see /vis/verbose guidance.
+  // G4VisManager* visManager = new G4VisExecutive("Quiet");
+  visManager->Initialize();
+#endif
+
+  // Get the pointer to the User Interface manager
+  //
+  G4UImanager* UImanager = G4UImanager::GetUIpointer(); 
+   
+  if ( macro.size() ) {
+     // Batch mode
+     G4String command = "/control/execute ";
+     UImanager->ApplyCommand(command+macro);
+  }
+  else // Define UI session for interactive mode
+  {
+#ifdef G4UI_USE
+     G4UIExecutive * ui = new G4UIExecutive(argc,argv,session);
+#ifdef G4VIS_USE
+     UImanager->ApplyCommand("/control/execute vis.mac");
+#else
+     UImanager->ApplyCommand("/control/execute OpNovice.in");
+#endif
+     if (ui->IsGUI())
+        UImanager->ApplyCommand("/control/execute gui.mac");
+     ui->SessionStart();
+     delete ui;
+#endif
+  }
+
+  // Job termination
+  // Free the store: user actions, physics_list and detector_description are
+  //                 owned and deleted by the run manager, so they should not
+  //                 be deleted in the main() program !
+
+#ifdef G4VIS_USE
+  delete visManager;
+#endif
+  delete runManager;
+
+  return 0;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/I3G4TankResponse.cxx b/I3G4TankResponse.cxx
new file mode 100644
index 0000000..cc34d1e
--- /dev/null
+++ b/I3G4TankResponse.cxx
@@ -0,0 +1,374 @@
+/**
+ * Copyright (C) 2009
+ * The IceCube collaboration
+ * ID: $Id: I3G4TankResponse.cxx 158929 2017-10-20 00:58:24Z cweaver $
+ *
+ * @file I3G4TankResponse.cxx
+ * @version $Revision: 158929 $
+ * @date $Date: 2017-10-20 01:58:24 +0100 (Fri, 20 Oct 2017) $
+ * @author Tilo Waldenmaier <tilo.waldenmaier@desy.de>, Thomas Melzig <thomas.melzig@desy.de>
+ *
+ * $LastChangedBy: cweaver $
+ */
+
+
+#include <topsimulator/interface/I3IceTopResponseFactory.h>
+#include <topsimulator/GeoFunctions.h>
+#include <g4-tankresponse/g4classes/G4IceTopTank.h>
+#include <g4-tankresponse/g4classes/G4Interface.h>
+#include <g4-tankresponse/I3G4TankResponse.h>
+//#include <pmt-simulator/I3PMTConstants.h>
+#include <phys-services/I3GSLRandomService.h>
+#include <dataclasses/physics/I3Particle.h>
+#include <dataclasses/I3Direction.h>
+#include <simclasses/I3MCPE.h>
+
+#include <cmath>
+
+const double I3G4TankResponse::PHOTONS_PER_VEM = 32880.;
+const double I3G4TankResponse::WEIGHTED_PHOTONS_PER_VEM = 1435.;
+const double I3G4TankResponse::MEV_PER_VEM = 186.6 * I3Units::MeV;
+const double I3G4TankResponse::VEM_THRESHOLD = 2000.0;
+
+I3G4TankResponse::I3G4TankResponse(I3Configuration& config, const I3Context& context, const TankKey& tankKey):
+  I3TankResponse(config, context, tankKey),  tauZirco_(7.94 * I3Units::m / I3Constants::c),
+  tauTyvek_(12.6 * I3Units::m / I3Constants::c),
+  decayTime_(0), safetyMargin_(0.3*I3Units::m), chargeScale_(1),
+  cherenkovMethod_(true), cherenkovWeighting_(false),
+  particlesBeforeThreshold_(0), particlesAfterThreshold_(0),
+  vemCounter_(0), visMacroFile_(""), pePerVEM_(), randomServiceName_(""),
+  g4Interface_(NULL), g4Tank_(NULL)
+{
+  AddParameter("TimeConstantZirco",
+               "Time constant for tanks with zirconium coating.",
+               tauZirco_);
+  AddParameter("TimeConstantTyvek",
+               "Time constant for tanks with a tykec liner.",
+               tauTyvek_);
+  AddParameter("VisMacro", "Geant4 visualization macro file",
+               visMacroFile_);
+  AddParameter("RandomServiceName", "Name of random service",
+               randomServiceName_);
+  AddParameter("SafetyMargin",
+               "Check if a particle hits the tank within safety margin.",
+               safetyMargin_);
+  AddParameter("ChargeScale",
+               "Scale all charges by this factor",
+               chargeScale_);
+  AddParameter("CherenkovMethod",
+               "Use the number of created cherenkov photons instead of energy deposition",
+               cherenkovMethod_);
+  AddParameter("CherenkovWeighting",
+               "Use position dependent weighting of the number of created cherenkov photons",
+               cherenkovWeighting_);
+}
+
+
+I3G4TankResponse::~I3G4TankResponse()
+{
+  if (G4Interface::GetInstance()) {
+    delete g4Interface_;
+  }
+}
+
+
+void I3G4TankResponse::Configure()
+{
+  log_info("Configuring Tank Response:");
+  
+  GetParameter("TimeConstantZirco", tauZirco_);
+  log_info(" + Time Constant (Zirconium): %.1f ns", tauZirco_ / I3Units::ns);
+  
+  GetParameter("TimeConstantTyvek", tauTyvek_);
+  log_info(" + Time Constant (Tyvek)    : %.1f ns", tauTyvek_ / I3Units::ns);
+  
+  GetParameter("VisMacro", visMacroFile_);
+#ifdef G4VIS_USE
+  log_info(" + Visualization macro      : %s", (visMacroFile_.empty() ? "DISABLED" : visMacroFile_.c_str()));
+#else
+  if (!visMacroFile_.empty()) {
+    log_fatal("VisMacro non-empty, but no Geant4 visualization in this build");
+  }
+#endif
+
+  GetParameter("SafetyMargin", safetyMargin_);
+  log_info(" + Safety margin            : %.1f m", safetyMargin_ / I3Units::m);
+  
+  GetParameter("ChargeScale", chargeScale_);
+  log_info(" + Charge scale             : %.2f", chargeScale_);
+  
+  GetParameter("CherenkovMethod", cherenkovMethod_);
+  log_info(" + Cherenkov method         : %s", cherenkovMethod_ ? "PHOTON COUNTING" : "ENERGY DEPOSIT");
+  
+  GetParameter("CherenkovWeighting", cherenkovWeighting_);
+  log_info(" + Cherenkov weighting      : %s", cherenkovWeighting_ ? "ENABLED" : "DISABLED");
+  
+  log_info(" + Saturation threshold     : %.1f VEM", VEM_THRESHOLD);
+  
+  GetParameter("RandomServiceName", randomServiceName_);
+  
+  // Look for a random generator service
+  if (randomServiceName_.empty())
+  {
+    randomService_ = I3RandomServicePtr(new I3GSLRandomService(0));
+    log_info("+ Random service           : I3GSLRandomService  (default)");
+  }
+  else
+  {
+    randomService_ = GetContext().Get<I3RandomServicePtr>(randomServiceName_);
+    log_info("+ Random service           : %s  (EXTERNAL)",  randomServiceName_.c_str());
+  }
+  if(!randomService_) log_fatal("Missing random service!");
+}
+
+
+double I3G4TankResponse::GetX() const
+{
+  return (g4Tank_ ? g4Tank_->GetX_I3() : NAN);
+}
+
+
+double I3G4TankResponse::GetY() const
+{
+  return (g4Tank_ ? g4Tank_->GetY_I3() : NAN);
+}
+
+
+double I3G4TankResponse::GetZ() const
+{
+  return (g4Tank_ ? g4Tank_->GetZ_I3() : NAN);
+}
+
+
+double I3G4TankResponse::GetTankRadius() const
+{
+  return (g4Tank_ ? g4Tank_->GetTankRadius_I3() : NAN);
+}
+
+
+double I3G4TankResponse::GetTankHeight() const
+{
+  return (g4Tank_ ? g4Tank_->GetTankHeight_I3() : NAN);
+}
+
+
+double I3G4TankResponse::GetSnowHeight() const
+{
+  return (g4Tank_ ? g4Tank_->GetSnowHeight_I3() : NAN);
+}
+
+
+void I3G4TankResponse::Initialize(const I3Geometry& geometry,
+                                  const I3Calibration& calib,
+                                  const I3DetectorStatus& status)
+{
+  // Get the stations geometry
+  const I3StationGeoMap& stationMap = geometry.stationgeo;
+  
+  // Get the dom status
+  const std::map<OMKey, I3DOMStatus>& domStatusMap = status.domStatus;
+  
+  // Get the VEM calibration
+  const std::map<OMKey, I3VEMCalibration>& vemCalMap = calib.vemCal;
+  
+  I3StationGeoMap::const_iterator station_iter = stationMap.find(tankKey_.string);
+  if (station_iter==stationMap.end()) {
+    log_fatal("Station %d doesn't exist in geometry!", tankKey_.string);
+    return;
+  }
+  
+  unsigned int tankID = tankKey_.tank == TankKey::TankA ? 0 : 1;
+  const I3TankGeo& tankGeo = station_iter->second.at(tankID);
+  
+  // Set decay time for specific tank coating
+  switch (tankGeo.tanktype) 
+  {
+    case I3TankGeo::Tyvek_Lined:
+      decayTime_ = tauTyvek_;
+      break;
+    case I3TankGeo::Zirconium_Lined:
+      decayTime_ = tauZirco_;
+      break;
+    default:
+      log_fatal_stream("Unknown type of tank " << tankKey_ << "! Is your GCD-file up-to-date?");
+      return;
+  }
+  
+  // Loop over all DOMs in the tank
+  pePerVEM_.clear();
+  I3Vector<OMKey>::const_iterator dom_iter;
+  for (dom_iter = tankGeo.omKeyList_.begin();
+       dom_iter != tankGeo.omKeyList_.end(); dom_iter++)
+  {
+    // Get the DOMStatus of this DOM to see if it is in the configuration
+    std::map<OMKey, I3DOMStatus>::const_iterator status_iter = domStatusMap.find(*dom_iter);
+    if(status_iter==domStatusMap.end()) continue;
+    
+    // Check if PMT is powered up
+    if(!(status_iter->second.pmtHV>0)) {
+      log_warn_stream("HV of " << *dom_iter << " is down!");
+    }
+    
+    // Get the VEM calibration
+    std::map<OMKey, I3VEMCalibration>::const_iterator vem_iter = vemCalMap.find(*dom_iter);
+    if (vem_iter==vemCalMap.end()) {
+      log_warn_stream("Missing VEM calibration of module " << *dom_iter << ". Skipping it!");
+      continue;
+    }
+    // Fill the map of pePerVEM values for each DOM
+    // The actual pePerVEM value is scaled down by the mean PE charge (0.85) to account for the
+    // asymmetric single photoelectron distribution --> This only works together with the
+    // pmt-simulator option "TreatIceTopAsInIce"=True!!!!
+    //pePerVEM_[*dom_iter] = (vem_iter->second.pePerVEM / vem_iter->second.corrFactor) /
+    //                        I3PMTConstants::MEAN_NORMALIZED_PE;
+    pePerVEM_[*dom_iter] = (vem_iter->second.pePerVEM / vem_iter->second.corrFactor) / MEAN_NORMALIZED_PE;
+  }
+  
+  g4Interface_ = G4Interface::GetInstance();
+  if (!g4Interface_) g4Interface_ = new G4Interface(visMacroFile_);
+  
+  g4Tank_ = new G4IceTopTank(tankKey_, geometry);
+  
+  g4Interface_->InstallTank(g4Tank_);
+}
+
+
+void I3G4TankResponse::BeginEvent(const I3Particle& primary)
+{
+  g4Interface_->InitializeEvent();
+  // Get the primary direction
+  I3Direction dir = primary.GetDir();
+  double x_dir = dir.GetX();
+  double y_dir = dir.GetY();
+  double z_dir = dir.GetZ();
+  // Get the vector pointing from tank to particle
+  double x_diff = primary.GetX() - g4Tank_->GetX_I3();
+  double y_diff = primary.GetY() - g4Tank_->GetY_I3();
+  double z_diff = primary.GetZ() - g4Tank_->GetZ_I3();
+  // Get the closest distance of this tank to primary track
+  double d = sqrt(pow(-x_dir * y_diff + x_diff * y_dir, 2) +
+                  pow( x_dir * z_diff - x_diff * z_dir, 2) +
+                  pow(-y_dir * z_diff + y_diff * z_dir, 2));
+  
+  double distToCore = 20.0 * I3Units::m;
+  if (d < distToCore) {
+    log_debug_stream("+ Tank " << tankKey_ << " lies closer than "
+                     << distToCore << " m to the primary: " << d << " m");
+  }
+  
+  // Set the VEM and particle counter to zero
+  vemCounter_ = 0.0;
+  particlesBeforeThreshold_ = 0;
+  particlesAfterThreshold_ = 0;
+}
+
+
+bool I3G4TankResponse::TrackParticle(const ExtendedI3Particle& particle, HitHistoCollection& hitHC, HitHistoCollection& cherHitCollection)
+{
+  if (GeoFunctions::IntersectCylinder(GetX(), GetY(), GetZ(), GetVirtualTankHeight(),
+                                      GetVirtualTankRadius(), particle))
+  {
+    if (vemCounter_ < VEM_THRESHOLD) {
+      particlesBeforeThreshold_++;
+    } else {
+      particlesAfterThreshold_++;
+      return true;
+    }
+    
+    g4Interface_->InjectParticle(particle);
+    
+    // Loop over all active DOMs in tank
+    double sum_vem = 0.;
+    int num_vem = 0;
+    std::map<OMKey, double>::const_iterator pePerVEM_iter;
+    for (pePerVEM_iter=pePerVEM_.begin(); pePerVEM_iter!=pePerVEM_.end(); ++pePerVEM_iter)
+    {
+      const OMKey& omKey = pePerVEM_iter->first;
+      ExtendedI3Particle iceTrack(particle); // copy constructor preserves ParticleID
+      iceTrack.SetTime(particle.GetTime() + g4Tank_->GetTime_I3(omKey));
+      iceTrack.SetShape(I3Particle::MCTrack);
+
+      double vem_mean = 0.;
+      double cher_mean = 0.;
+      if (cherenkovMethod_) {
+        if (cherenkovWeighting_) {
+          cher_mean = g4Tank_->GetNumCherenkovWeight(omKey);
+          vem_mean = cher_mean / WEIGHTED_PHOTONS_PER_VEM;
+        } else {
+          cher_mean = g4Tank_->GetNumCherenkov(omKey);
+          vem_mean = cher_mean / PHOTONS_PER_VEM;
+        }
+      } else {
+        cher_mean = g4Tank_->GetEDep_I3(omKey);
+        vem_mean = cher_mean / MEV_PER_VEM;
+      }
+      sum_vem += vem_mean;
+      if (vem_mean > 0) num_vem++;
+      
+      // Calculate mean number of photoelectrons according to actual VEM calibration
+      double npe_mean = vem_mean * pePerVEM_iter->second * chargeScale_;
+      
+      // Dial npe according to Poisson distribution
+      int num_hits = randomService_->Poisson(npe_mean);
+
+      cherHitCollection.GetHitHisto(omKey).Fill(iceTrack.GetStartTime(), int(cher_mean), iceTrack);
+      GenerateHits(num_hits, iceTrack, hitHC.GetHitHisto(omKey));
+      //log_warn("DOM: %s, VEM: %.2f, ",omKey.str().c_str(), vem_mean);
+    }
+    if (num_vem > 0) vemCounter_ += sum_vem / num_vem;
+    else return false;
+
+    return true;
+  }
+  log_trace("No cylinder intersection. Tank position (%f, %f, %f), Particle position (%f, %f, %f) and direction (%f, %f, %f). Height: %f. Radius: %f.",
+	    GetX(), GetY(), GetZ(),
+	    particle.GetPos().GetX(), particle.GetPos().GetY(), particle.GetPos().GetZ(),
+	    particle.GetDir().GetX(), particle.GetDir().GetY(), particle.GetDir().GetZ(),
+	    GetVirtualTankHeight(), GetVirtualTankRadius());
+  return false;
+}
+
+
+void I3G4TankResponse::GenerateHits(int npe, const ExtendedI3Particle& p, HitHisto& hitHisto)
+{
+  if (npe == 0)
+    return;
+
+  // draw individual npe with time delay from exponential distribution
+  const double t0 = p.GetStartTime();
+
+  if (npe == 1) {
+    // no hokus pokus for single pe
+    hitHisto.Fill(t0 + randomService_->Exp(decayTime_), 1, p);
+  } else {
+    // fill times into a vector and pass that to HitHisto
+    std::vector<double> times(npe);
+    for (int j = 0; j < npe; ++j)
+      times[j] = t0 + randomService_->Exp(decayTime_);
+    hitHisto.Fill(times, p);
+  }
+}
+
+
+void I3G4TankResponse::EndEvent(HitHistoCollection &hitHC, HitHistoCollection& cherHitCollection)
+{
+  g4Interface_->TerminateEvent();
+  if (particlesAfterThreshold_ > 0) {
+    double scalingFactor = static_cast<double>(particlesAfterThreshold_) /
+                           static_cast<double>(particlesBeforeThreshold_) + 1.0;
+
+    std::map<OMKey, double>::const_iterator pePerVEM_iter;
+    for (pePerVEM_iter = pePerVEM_.begin(); pePerVEM_iter != pePerVEM_.end(); ++pePerVEM_iter)
+    {
+      hitHC.GetHitHisto(pePerVEM_iter->first).Scale(scalingFactor);
+      cherHitCollection.GetHitHisto(pePerVEM_iter->first).Scale(scalingFactor);
+    }
+    log_trace_stream(tankKey_ << " reached VEM-threshold, particles before: "
+                              << particlesBeforeThreshold_ << ", after: "
+                              << particlesAfterThreshold_ << ", scaling-factor: "
+                              << scalingFactor);
+  }
+}
+
+
+I3_SERVICE_FACTORY(I3IceTopResponseFactory<I3G4TankResponse>);
diff --git a/include/G4BeamTestDetectorConstruction.h b/include/G4BeamTestDetectorConstruction.h
new file mode 100644
index 0000000..ce53256
--- /dev/null
+++ b/include/G4BeamTestDetectorConstruction.h
@@ -0,0 +1,45 @@
+#ifndef _TOPSIMULATOR_G4BEAMTESTDETECTORCONSTRUCTION_H_
+#define _TOPSIMULATOR_G4BEAMTESTDETECTORCONSTRUCTION_H_
+
+#include <G4VUserDetectorConstruction.hh>
+#include <G4ThreeVector.hh>
+
+// class G4BeamTestTank;
+
+class G4BeamTestDetectorConstruction: public G4VUserDetectorConstruction
+{
+ public:
+  G4BeamTestDetectorConstruction();
+
+  ~G4BeamTestDetectorConstruction();
+
+  G4VPhysicalVolume* Construct();
+
+  void SetVerboseLevel(G4int level) {verboseLevel_=level;}
+
+  /* void InstallTank(G4BeamTestTank* tank) {tankList_.push_back(tank);} */
+
+  const G4ThreeVector& GetWorldOrigin() const {return origin_;}
+
+ private:
+
+  void CreateMaterials();
+
+  void CreateAir();
+  /* void CreateIce(); */
+  /* void CreateSnow(); */
+  void CreateWater();
+  void CreatePlastic();
+  /* void CreateTyvek(); */
+  /* void CreatePerlite(); */
+  void CreateGlassSphere();
+  void CreateEffectiveDOMMaterial();
+
+  G4ThreeVector origin_;
+
+  G4int verboseLevel_;
+
+  // std::vector<G4BeamTestTank*> tankList_;
+};
+
+#endif
diff --git a/include/G4BeamTestEMPhysics.h b/include/G4BeamTestEMPhysics.h
new file mode 100644
index 0000000..797719e
--- /dev/null
+++ b/include/G4BeamTestEMPhysics.h
@@ -0,0 +1,64 @@
+/*
+ * copyright  (C) 2010
+ * The Icecube Collaboration
+ *
+ * $Id: G4BeamTestEMPhysics.h 154687 2017-04-05 15:46:57Z jgonzalez $
+ *
+ * @version $Revision: 154687 $
+ * @date $LastChangedDate: 2017-04-05 16:46:57 +0100 (Wed, 05 Apr 2017) $
+ * @author Fabian Kislat <fabian.kislat@desy.de>, Peter Nießen (tanktop)
+ * Last changed by: $LastChangedBy: jgonzalez $
+ */
+
+#ifndef G4TANKRESPONSE_G4BEAMTESTEMPHYSICS_H_INCLUDED
+#define G4TANKRESPONSE_G4BEAMTESTEMPHYSICS_H_INCLUDED
+
+#include <globals.hh>
+#include <G4VPhysicsConstructor.hh>
+#include <G4PhotoElectricEffect.hh>
+#include <G4ComptonScattering.hh>
+#include <G4GammaConversion.hh>
+#include <G4eMultipleScattering.hh>
+#include <G4eIonisation.hh>
+#include <G4eBremsstrahlung.hh>
+#include <G4eplusAnnihilation.hh>
+
+/**
+   @class G4BeamTestEMPhysics
+   @brief Electromagnetic physics. Used only if Geant4 version is earlier than 4.10.
+
+   This class implements the electromagnetic interactions
+   - Photoelectric effect
+   - Compton scattering
+   - Gamma conversion
+   - Multiple scattering
+   - Ionisation/Bremsstrahlung for electrons
+   - Positron annihilation
+*/
+class G4BeamTestEMPhysics : public G4VPhysicsConstructor {
+public:
+  G4BeamTestEMPhysics();
+  ~G4BeamTestEMPhysics();
+
+  void ConstructParticle();
+  void ConstructProcess();
+
+private:
+  // Gamma physics
+  G4PhotoElectricEffect photoEffect;
+  G4ComptonScattering comptonEffect;
+  G4GammaConversion pairProduction;
+
+  // Electron physics
+  G4eMultipleScattering electronMultipleScattering;
+  G4eIonisation electronIonisation;
+  G4eBremsstrahlung electronBremsStrahlung;
+  
+  //Positron physics
+  G4eMultipleScattering positronMultipleScattering;
+  G4eIonisation positronIonisation; 
+  G4eBremsstrahlung positronBremsStrahlung;  
+  G4eplusAnnihilation annihilation;
+};
+
+#endif // G4TANKRESPONSE_G4BEAMTESTEMPHYSICS_H_INCLUDED
diff --git a/include/G4BeamTestGeneralPhysics.h b/include/G4BeamTestGeneralPhysics.h
new file mode 100644
index 0000000..012bf40
--- /dev/null
+++ b/include/G4BeamTestGeneralPhysics.h
@@ -0,0 +1,35 @@
+/*
+ * copyright  (C) 2010
+ * The Icecube Collaboration
+ *
+ * $Id: G4BeamTestGeneralPhysics.h 149388 2016-08-18 21:50:04Z jgonzalez $
+ *
+ * @version $Revision: 149388 $
+ * @date $LastChangedDate: 2016-08-18 22:50:04 +0100 (Thu, 18 Aug 2016) $
+ * @author Fabian Kislat <fabian.kislat@desy.de> Last changed by: $LastChangedBy: jgonzalez $
+ */
+
+#ifndef G4TANKRESPONSE_G4BEAMTESTGENERALPHYSICS_H_INCLUDED
+#define G4TANKRESPONSE_G4BEAMTESTGENERALPHYSICS_H_INCLUDED
+
+#include <globals.hh>
+#include <G4VPhysicsConstructor.hh>
+#include <G4Decay.hh>
+
+/**
+ * Implementation of G4VPhysicsConstructor for decay processes and geantino.
+ */
+class G4BeamTestGeneralPhysics : public G4VPhysicsConstructor
+{
+public:
+  G4BeamTestGeneralPhysics();
+  virtual ~G4BeamTestGeneralPhysics();
+
+  void ConstructParticle();
+  void ConstructProcess();
+
+private:
+  G4Decay decay;
+};
+
+#endif // G4TANKRESPONSE_G4BEAMTESTGENERALPHYSICS_H_INCLUDED
diff --git a/include/G4BeamTestHadronPhysics.h b/include/G4BeamTestHadronPhysics.h
new file mode 100644
index 0000000..0fd9868
--- /dev/null
+++ b/include/G4BeamTestHadronPhysics.h
@@ -0,0 +1,338 @@
+/*
+ * copyright  (C) 2010
+ * The Icecube Collaboration
+ *
+ * $Id: G4BeamTestHadronPhysics.h 154687 2017-04-05 15:46:57Z jgonzalez $
+ *
+ * @version $Revision: 154687 $
+ * @date $LastChangedDate: 2017-04-05 16:46:57 +0100 (Wed, 05 Apr 2017) $
+ * @author Fabian Kislat <fabian.kislat@desy.de> Last changed by: $LastChangedBy: jgonzalez $
+ */
+
+#ifndef G4TANKRESPONSE_G4BEAMTESTHADRONPHYSICS_H_INCLUDED
+#define G4TANKRESPONSE_G4BEAMTESTHADRONPHYSICS_H_INCLUDED
+
+#include <vector>
+
+#include <globals.hh>
+
+#include <G4VPhysicsConstructor.hh>
+
+#include <G4hMultipleScattering.hh>
+#include <G4hIonisation.hh>
+
+#include <G4HadronElasticProcess.hh>
+#include <G4HadronFissionProcess.hh>
+#include <G4HadronCaptureProcess.hh>
+
+#include <G4PionPlusInelasticProcess.hh>
+#include <G4PionMinusInelasticProcess.hh>
+#include <G4KaonPlusInelasticProcess.hh>
+#include <G4KaonZeroSInelasticProcess.hh>
+#include <G4KaonZeroLInelasticProcess.hh>
+#include <G4KaonMinusInelasticProcess.hh>
+#include <G4ProtonInelasticProcess.hh>
+#include <G4AntiProtonInelasticProcess.hh>
+#include <G4NeutronInelasticProcess.hh>
+#include <G4AntiNeutronInelasticProcess.hh>
+#include <G4LambdaInelasticProcess.hh>
+#include <G4AntiLambdaInelasticProcess.hh>
+#include <G4SigmaPlusInelasticProcess.hh>
+#include <G4SigmaMinusInelasticProcess.hh>
+#include <G4AntiSigmaPlusInelasticProcess.hh>
+#include <G4AntiSigmaMinusInelasticProcess.hh>
+#include <G4XiZeroInelasticProcess.hh>
+#include <G4XiMinusInelasticProcess.hh>
+#include <G4AntiXiZeroInelasticProcess.hh>
+#include <G4AntiXiMinusInelasticProcess.hh>
+#include <G4DeuteronInelasticProcess.hh>
+#include <G4TritonInelasticProcess.hh>
+#include <G4AlphaInelasticProcess.hh>
+#include <G4OmegaMinusInelasticProcess.hh>
+#include <G4AntiOmegaMinusInelasticProcess.hh>
+
+// Low-energy Models
+#include <G4LElastic.hh>   
+#include <G4LFission.hh>
+#include <G4LCapture.hh>
+
+#include <G4LEPionPlusInelastic.hh>
+#include <G4LEPionMinusInelastic.hh>
+#include <G4LEKaonPlusInelastic.hh>
+#include <G4LEKaonZeroSInelastic.hh>
+#include <G4LEKaonZeroLInelastic.hh>
+#include <G4LEKaonMinusInelastic.hh>
+#include <G4LEProtonInelastic.hh>
+#include <G4LEAntiProtonInelastic.hh>
+#include <G4LENeutronInelastic.hh>
+#include <G4LEAntiNeutronInelastic.hh>
+#include <G4LELambdaInelastic.hh>
+#include <G4LEAntiLambdaInelastic.hh>
+#include <G4LESigmaPlusInelastic.hh>
+#include <G4LESigmaMinusInelastic.hh>
+#include <G4LEAntiSigmaPlusInelastic.hh>
+#include <G4LEAntiSigmaMinusInelastic.hh>
+#include <G4LEXiZeroInelastic.hh>
+#include <G4LEXiMinusInelastic.hh>
+#include <G4LEAntiXiZeroInelastic.hh>
+#include <G4LEAntiXiMinusInelastic.hh>
+#include <G4LEDeuteronInelastic.hh>
+#include <G4LETritonInelastic.hh>
+#include <G4LEAlphaInelastic.hh>
+#include <G4LEOmegaMinusInelastic.hh>
+#include <G4LEAntiOmegaMinusInelastic.hh>
+
+// High-energy Models
+
+#include <G4HEPionPlusInelastic.hh>
+#include <G4HEPionMinusInelastic.hh>
+#include <G4HEKaonPlusInelastic.hh>
+#include <G4HEKaonZeroInelastic.hh>
+#include <G4HEKaonZeroInelastic.hh>
+#include <G4HEKaonMinusInelastic.hh>
+#include <G4HEProtonInelastic.hh>
+#include <G4HEAntiProtonInelastic.hh>
+#include <G4HENeutronInelastic.hh>
+#include <G4HEAntiNeutronInelastic.hh>
+#include <G4HELambdaInelastic.hh>
+#include <G4HEAntiLambdaInelastic.hh>
+#include <G4HESigmaPlusInelastic.hh>
+#include <G4HESigmaMinusInelastic.hh>
+#include <G4HEAntiSigmaPlusInelastic.hh>
+#include <G4HEAntiSigmaMinusInelastic.hh>
+#include <G4HEXiZeroInelastic.hh>
+#include <G4HEXiMinusInelastic.hh>
+#include <G4HEAntiXiZeroInelastic.hh>
+#include <G4HEAntiXiMinusInelastic.hh>
+#include <G4HEOmegaMinusInelastic.hh>
+#include <G4HEAntiOmegaMinusInelastic.hh>
+
+// Stopping processes
+#include <G4AntiProtonAnnihilationAtRest.hh>
+#include <G4AntiNeutronAnnihilationAtRest.hh>
+
+#ifdef TRIUMF_STOP_PIMINUS
+#include <G4PionMinusAbsorptionAtRest.hh>
+#else
+#include <G4PiMinusAbsorptionAtRest.hh>
+#endif
+#ifdef TRIUMF_STOP_KMINUS
+#include <G4KaonMinusAbsorption.hh>
+#else
+#include <G4KaonMinusAbsorptionAtRest.hh>
+#endif
+
+// quark gluon string model with chips afterburner.
+#include <G4TheoFSGenerator.hh>
+#include <G4ExcitationHandler.hh>
+#include <G4PreCompoundModel.hh>
+#include <G4GeneratorPrecompoundInterface.hh>
+#include <G4QGSModel.hh>
+#include <G4QGSParticipants.hh>
+#include <G4QGSMFragmentation.hh>
+#include <G4ExcitedStringDecay.hh>
+
+/**
+   @brief Hadron physics. Used only if Geant4 version is earlier than 4.10.
+   @author GEANT4/Peter Niessen
+   @date Sat Jul 24 23:53:47 EDT 2004
+
+   Many processes. You're encouraged to check the source.
+*/
+class G4BeamTestHadronPhysics : public G4VPhysicsConstructor {
+public: 
+
+  /**
+   * The constructor
+   */
+  G4BeamTestHadronPhysics();
+
+  /**
+   * The virtual destructor
+   */
+  virtual ~G4BeamTestHadronPhysics();
+
+  /**
+   * This method will be invoked in the Construct() method. 
+   * each particle type will be instantiated
+   */
+  virtual void ConstructParticle();
+ 
+  /**
+   * This method will be invoked in the Construct() method.
+   * each physics process will be instantiated and
+   * registered to the process manager of each particle type 
+   */
+  virtual void ConstructProcess();
+
+protected:
+
+  // Elastic Process
+  G4HadronElasticProcess elasticProcess_;
+  G4LElastic *elasticModel_;
+  
+  // Pi + 
+  G4PionPlusInelasticProcess pionPlusInelastic_;
+  G4LEPionPlusInelastic *lePionPlusModel_;
+  G4hMultipleScattering pionPlusMult_;
+  G4hIonisation pionPlusIonisation_;
+
+  // Pi -
+  G4PionMinusInelasticProcess pionMinusInelastic_;
+  G4LEPionMinusInelastic *lePionMinusModel_;
+  G4hMultipleScattering pionMinusMult_;
+  G4hIonisation pionMinusIonisation_;
+#ifdef TRIUMF_STOP_PIMINUS
+  G4PionMinusAbsorptionAtRest pionMinusAbsorption_;
+#else
+  G4PiMinusAbsorptionAtRest pionMinusAbsorption_;
+#endif
+
+  // pi+ and pi-
+  G4TheoFSGenerator theoModel_;
+  G4ExcitationHandler excitationHandler_;
+  G4PreCompoundModel *preEquilib_;
+  G4GeneratorPrecompoundInterface cascade_;
+  G4QGSModel< G4QGSParticipants > stringModel_;
+  G4QGSMFragmentation fragmentation_;
+  G4ExcitedStringDecay *stringDecay_;
+
+  // K + 
+  G4KaonPlusInelasticProcess kaonPlusInelastic_;
+  G4LEKaonPlusInelastic *leKaonPlusModel_;
+  G4HEKaonPlusInelastic *heKaonPlusModel_;
+  G4hMultipleScattering kaonPlusMult_;
+  G4hIonisation kaonPlusIonisation_;
+	
+  // K -
+  G4KaonMinusInelasticProcess kaonMinusInelastic_;
+  G4LEKaonMinusInelastic *leKaonMinusModel_;
+  G4HEKaonMinusInelastic *heKaonMinusModel_;
+  G4hMultipleScattering kaonMinusMult_;
+  G4hIonisation kaonMinusIonisation_;
+#ifdef TRIUMF_STOP_KMINUS
+  G4KaonMinusAbsorption kaonMinusAbsorption_;
+#else
+  G4PiMinusAbsorptionAtRest kaonMinusAbsorption_;
+#endif
+
+  // K0L
+  G4KaonZeroLInelasticProcess kaonZeroLInelastic_;
+  G4LEKaonZeroLInelastic *leKaonZeroLModel_;
+  G4HEKaonZeroInelastic *heKaonZeroLModel_;
+
+  // K0S
+  G4KaonZeroSInelasticProcess kaonZeroSInelastic_;
+  G4LEKaonZeroSInelastic *leKaonZeroSModel_;
+  G4HEKaonZeroInelastic *heKaonZeroSModel_;
+
+  // Proton
+  G4ProtonInelasticProcess protonInelastic_;
+  G4LEProtonInelastic *leProtonModel_;
+  G4HEProtonInelastic *heProtonModel_;
+  G4hMultipleScattering protonMult_;
+  G4hIonisation protonIonisation_;
+ 
+  // anti-proton
+  G4AntiProtonInelasticProcess antiProtonInelastic_;
+  G4LEAntiProtonInelastic *leAntiProtonModel_;
+  G4HEAntiProtonInelastic *heAntiProtonModel_;
+  G4hMultipleScattering antiProtonMult_;
+  G4hIonisation antiProtonIonisation_;
+  G4AntiProtonAnnihilationAtRest antiProtonAnnihilation_;
+    
+  // neutron
+  G4NeutronInelasticProcess neutronInelastic_;
+  G4LENeutronInelastic *leNeutronModel_;
+  G4HENeutronInelastic *heNeutronModel_;
+  G4HadronFissionProcess neutronFission_;
+  G4LFission *neutronFissionModel_;
+  G4HadronCaptureProcess neutronCapture_;
+  G4LCapture *neutronCaptureModel_;
+
+
+  // anti-neutron
+  G4AntiNeutronInelasticProcess antiNeutronInelastic_;
+  G4LEAntiNeutronInelastic *leAntiNeutronModel_;
+  G4HEAntiNeutronInelastic *heAntiNeutronModel_;
+  G4AntiNeutronAnnihilationAtRest antiNeutronAnnihilation_;
+   
+  // Lambda
+  G4LambdaInelasticProcess lambdaInelastic_;
+  G4LELambdaInelastic *leLambdaModel_;
+  G4HELambdaInelastic *heLambdaModel_;
+  
+  // AntiLambda
+  G4AntiLambdaInelasticProcess antiLambdaInelastic_;
+  G4LEAntiLambdaInelastic *leAntiLambdaModel_;
+  G4HEAntiLambdaInelastic *heAntiLambdaModel_;
+  
+  // SigmaMinus
+  G4SigmaMinusInelasticProcess sigmaMinusInelastic_;
+  G4LESigmaMinusInelastic *leSigmaMinusModel_;
+  G4HESigmaMinusInelastic *heSigmaMinusModel_;
+  G4hMultipleScattering sigmaMinusMult_;
+  G4hIonisation sigmaMinusIonisation_;
+  
+  // AntiSigmaMinus
+  G4AntiSigmaMinusInelasticProcess antiSigmaMinusInelastic_;
+  G4LEAntiSigmaMinusInelastic *leAntiSigmaMinusModel_;
+  G4HEAntiSigmaMinusInelastic *heAntiSigmaMinusModel_;
+  G4hMultipleScattering antiSigmaMinusMult_;
+  G4hIonisation antiSigmaMinusIonisation_;
+   
+  // SigmaPlus
+  G4SigmaPlusInelasticProcess sigmaPlusInelastic_;
+  G4LESigmaPlusInelastic *leSigmaPlusModel_;
+  G4HESigmaPlusInelastic *heSigmaPlusModel_;
+  G4hMultipleScattering sigmaPlusMult_;
+  G4hIonisation sigmaPlusIonisation_;
+  
+  // AntiSigmaPlus
+  G4AntiSigmaPlusInelasticProcess antiSigmaPlusInelastic_;
+  G4LEAntiSigmaPlusInelastic *leAntiSigmaPlusModel_;
+  G4HEAntiSigmaPlusInelastic *heAntiSigmaPlusModel_;
+  G4hMultipleScattering antiSigmaPlusMult_;
+  G4hIonisation antiSigmaPlusIonisation_;
+  
+  // XiZero
+  G4XiZeroInelasticProcess xiZeroInelastic_;
+  G4LEXiZeroInelastic *leXiZeroModel_;
+  G4HEXiZeroInelastic *heXiZeroModel_;
+  
+  // AntiXiZero
+  G4AntiXiZeroInelasticProcess  antiXiZeroInelastic_;
+  G4LEAntiXiZeroInelastic*  leAntiXiZeroModel_;
+  G4HEAntiXiZeroInelastic*  heAntiXiZeroModel_;
+  
+  // XiMinus
+  G4XiMinusInelasticProcess xiMinusInelastic_;
+  G4LEXiMinusInelastic *leXiMinusModel_;
+  G4HEXiMinusInelastic *heXiMinusModel_;
+  G4hMultipleScattering xiMinusMult_;
+  G4hIonisation xiMinusIonisation_;
+
+  // AntiXiMinus
+  G4AntiXiMinusInelasticProcess antiXiMinusInelastic_;
+  G4LEAntiXiMinusInelastic *leAntiXiMinusModel_;
+  G4HEAntiXiMinusInelastic *heAntiXiMinusModel_;
+  G4hMultipleScattering antiXiMinusMult_;
+  G4hIonisation antiXiMinusIonisation_;
+  
+  // OmegaMinus
+  G4OmegaMinusInelasticProcess omegaMinusInelastic_;
+  G4LEOmegaMinusInelastic *leOmegaMinusModel_;
+  G4HEOmegaMinusInelastic *heOmegaMinusModel_;
+  G4hMultipleScattering omegaMinusMult_;
+  G4hIonisation omegaMinusIonisation_;
+   
+  // AntiOmegaMinus
+  G4AntiOmegaMinusInelasticProcess antiOmegaMinusInelastic_;
+  G4LEAntiOmegaMinusInelastic *leAntiOmegaMinusModel_;
+  G4HEAntiOmegaMinusInelastic *heAntiOmegaMinusModel_;
+  G4hMultipleScattering antiOmegaMinusMult_;
+  G4hIonisation antiOmegaMinusIonisation_;
+};
+
+
+#endif // G4TANKRESPONSE_G4BEAMTESTHADRONPHYSICS_H_INCLUDED
diff --git a/include/G4BeamTestIonPhysics.h b/include/G4BeamTestIonPhysics.h
new file mode 100644
index 0000000..276b0d5
--- /dev/null
+++ b/include/G4BeamTestIonPhysics.h
@@ -0,0 +1,100 @@
+/*
+ * copyright  (C) 2010
+ * The Icecube Collaboration
+ *
+ * $Id: G4BeamTestIonPhysics.h 154687 2017-04-05 15:46:57Z jgonzalez $
+ *
+ * @version $Revision: 154687 $
+ * @date $LastChangedDate: 2017-04-05 16:46:57 +0100 (Wed, 05 Apr 2017) $
+ * @author Fabian Kislat <fabian.kislat@desy.de>, Peter Nießen (tanktop)
+ * Last changed by: $LastChangedBy: jgonzalez $
+ */
+
+#ifndef G4TANKRESPONSE_G4BEAMTESTIONPHYSICS_H_INCLUDED
+#define G4TANKRESPONSE_G4BEAMTESTIONPHYSICS_H_INCLUDED
+
+#include <globals.hh>
+#include <G4ios.hh>
+#include <G4VPhysicsConstructor.hh>
+#include <G4HadronElasticProcess.hh>
+#include <G4LElastic.hh>
+#include <G4DeuteronInelasticProcess.hh>
+#include <G4LEDeuteronInelastic.hh>
+#include <G4TritonInelasticProcess.hh>
+#include <G4LETritonInelastic.hh>
+#include <G4AlphaInelasticProcess.hh>
+#include <G4LEAlphaInelastic.hh>
+#include <G4hIonisation.hh>
+#include <G4ionIonisation.hh>
+#include <G4hMultipleScattering.hh>
+
+/**
+   @class G4BeamTestIonPhysics
+   @brief Ion physics. Used only if Geant4 version is earlier than 4.10.
+   @author GEANT4/Peter Niessen
+   @date Sun Jul 25 00:24:42 EDT 2004
+
+   The ion physics. Check the source for details.
+*/
+class G4BeamTestIonPhysics : public G4VPhysicsConstructor {
+
+public:
+
+  /**
+   * The constructor
+   */
+  G4BeamTestIonPhysics();
+
+  /**
+   * The virtual destructor
+   */
+  virtual ~G4BeamTestIonPhysics();
+
+  /**
+   * This method will be invoked in the Construct() method. 
+   * each particle type will be instantiated
+   */
+  virtual void ConstructParticle();
+ 
+  /**
+   * This method will be invoked in the Construct() method.
+   * each physics process will be instantiated and
+   * registered to the process manager of each particle type 
+   */
+  virtual void ConstructProcess();
+
+protected:
+  // Elastic Process
+  G4HadronElasticProcess      elasticProcess_;
+  G4LElastic*                 elasticModel_;
+
+  // Generic Ion physics
+  G4hMultipleScattering        ionMultipleScattering_;
+  G4ionIonisation             ionIonisation_;
+
+  // Deuteron physics
+  G4hMultipleScattering        deuteronMultipleScattering_;
+  G4hIonisation               deuteronIonisation_;
+  G4DeuteronInelasticProcess  deuteronProcess_;
+  G4LEDeuteronInelastic*      deuteronModel_;
+
+  // Triton physics
+  G4hMultipleScattering        tritonMultipleScattering_;
+  G4hIonisation               tritonIonisation_;
+  G4TritonInelasticProcess    tritonProcess_;
+  G4LETritonInelastic*        tritonModel_;
+  
+  // Alpha physics
+  G4hMultipleScattering        alphaMultipleScattering_;
+  G4hIonisation               alphaIonisation_;
+  G4AlphaInelasticProcess     alphaProcess_;
+  G4LEAlphaInelastic*         alphaModel_;
+
+  // He3 physics
+  G4hMultipleScattering        he3MultipleScattering_;
+  G4hIonisation               he3Ionisation_;
+
+};
+
+
+#endif // G4TANKRESPONSE_G4BEAMTESTIONPHYSICS_H_INCLUDED
diff --git a/include/G4BeamTestMuonPhysics.h b/include/G4BeamTestMuonPhysics.h
new file mode 100644
index 0000000..de09d73
--- /dev/null
+++ b/include/G4BeamTestMuonPhysics.h
@@ -0,0 +1,71 @@
+/*
+ * copyright  (C) 2010
+ * The Icecube Collaboration
+ *
+ * $Id: G4BeamTestMuonPhysics.h 154687 2017-04-05 15:46:57Z jgonzalez $
+ *
+ * @version $Revision: 154687 $
+ * @date $LastChangedDate: 2017-04-05 16:46:57 +0100 (Wed, 05 Apr 2017) $
+ * @author Fabian Kislat <fabian.kislat@desy.de>, Peter Nießen (tanktop)
+ * Last changed by: $LastChangedBy: jgonzalez $
+ */
+
+#ifndef G4TANKRESPONSE_G4BEAMTESTMUONPHYSICS_H_INCLUDED
+#define G4TANKRESPONSE_G4BEAMTESTMUONPHYSICS_H_INCLUDED
+
+#include <globals.hh>
+#include <G4VPhysicsConstructor.hh>
+#include <G4MuMultipleScattering.hh>
+#include <G4MuBremsstrahlung.hh>
+#include <G4MuPairProduction.hh>
+#include <G4MuIonisation.hh>
+#include <G4hIonisation.hh>
+#include <G4MuonMinusCaptureAtRest.hh>
+
+/**
+   @class G4BeamTestMuonPhysics
+   @brief Muon/tau Physics. Used only if Geant4 version is earlier than 4.10.
+
+   This class implements the physics processes for the muons. For
+   muons, it contains
+   - Ionisation
+   - Multiple scattering
+   - Bremsstrahlung
+   - Pair production
+   - Capture at rest (mu-)
+   For taus, it does
+   - Multiple scattering
+   - Ionisation
+*/
+class G4BeamTestMuonPhysics : public G4VPhysicsConstructor
+{
+public: 
+  G4BeamTestMuonPhysics();
+  ~G4BeamTestMuonPhysics();
+  
+  void ConstructParticle();
+  void ConstructProcess();
+
+private:
+  // Muon physics
+  G4MuIonisation         muPlusIonisation_;
+  G4MuMultipleScattering   muPlusMultipleScattering_;
+  G4MuBremsstrahlung     muPlusBremsstrahlung_;
+  G4MuPairProduction     muPlusPairProduction_;
+
+  G4MuIonisation         muMinusIonisation_;
+  G4MuMultipleScattering   muMinusMultipleScattering_;
+  G4MuBremsstrahlung     muMinusBremsstrahlung_;
+  G4MuPairProduction     muMinusPairProduction_;
+
+  G4MuonMinusCaptureAtRest muMinusCaptureAtRest_;
+
+  // Tau physics
+  G4MuMultipleScattering   tauPlusMultipleScattering_;
+  G4hIonisation          tauPlusIonisation_;
+
+  G4MuMultipleScattering   tauMinusMultipleScattering_;
+  G4hIonisation          tauMinusIonisation_;
+};
+
+#endif // G4TANKRESPONSE_G4BEAMTESTMUONPHYSICS_H_INCLUDED
diff --git a/include/G4BeamTestPhysicsList.h b/include/G4BeamTestPhysicsList.h
new file mode 100644
index 0000000..e253e33
--- /dev/null
+++ b/include/G4BeamTestPhysicsList.h
@@ -0,0 +1,35 @@
+/**
+ * Copyright (C) 2009
+ * The IceCube collaboration
+ * ID: $Id: G4BeamTestPhysicsList.h 149388 2016-08-18 21:50:04Z jgonzalez $
+ *
+ * @file G4BeamTestPhysicsList.h
+ * @version $Rev: 149388 $
+ * @date $Date: 2016-08-18 22:50:04 +0100 (Thu, 18 Aug 2016) $
+ * @author Tilo Waldenmaier, Thomas Melzig, Fabian Kislat
+ */
+
+#ifndef _TOPSIMULATOR_G4BEAMTESTPHYSICSLIST_H
+#define _TOPSIMULATOR_G4BEAMTESTPHYSICSLIST_H
+
+#include <G4VModularPhysicsList.hh>
+
+/**
+ * Implementation of G4VModularPhysicsList. The top-level physics list. It combines all the other physics lists:
+ *
+ *  \li G4BeamTestGeneralPhysics
+ *  \li G4BeamTestHadronPhysics
+ *  \li G4BeamTestIonPhysics
+ *  \li G4BeamTestMuonPhysics
+ *  \li G4BeamTestEMPhysics
+ */
+class G4BeamTestPhysicsList: public G4VModularPhysicsList
+{
+public:
+  G4BeamTestPhysicsList();
+  ~G4BeamTestPhysicsList();
+private:
+  void SetCuts();
+};
+
+#endif
diff --git a/include/G4BeamTestRunManager.h b/include/G4BeamTestRunManager.h
new file mode 100644
index 0000000..c142b38
--- /dev/null
+++ b/include/G4BeamTestRunManager.h
@@ -0,0 +1,46 @@
+/**
+ * Copyright (C) 2009
+ * The IceCube collaboration
+ * ID: $Id: G4BeamTestRunManager.h 149391 2016-08-18 21:58:04Z jgonzalez $
+ *
+ * @file G4BeamTestRunManager.h
+ * @version $Rev: 149391 $
+ * @date $Date: 2016-08-18 22:58:04 +0100 (Thu, 18 Aug 2016) $
+ * @author Tilo Waldenmaier
+ */
+
+
+#ifndef TOPSIMULATOR_G4BEAMTESTRUNMANAGER_H
+#define TOPSIMULATOR_G4BEAMTESTRUNMANAGER_H
+
+#include <G4RunManager.hh>
+
+class G4ParticleGun;
+
+/**
+ * Implementation of G4RunManager
+ */
+class G4BeamTestRunManager: public G4RunManager
+{
+ public:
+  G4BeamTestRunManager();
+
+  static G4BeamTestRunManager* GetIceTopRunManager() {return (G4IceTopRunManager*)GetRunManager();}
+
+  // Disable BeamOn
+  void BeamOn(G4int n_event,const char* macroFile=0,G4int n_select=-1);
+
+  void InitializeRun();
+  void InjectParticle(G4ParticleGun* particleGun);
+  void TerminateRun();
+
+ protected:
+  G4Event* GenerateEvent(G4int i_event);
+
+ private:
+  // This method is an exact copy of UpdateScoring which is private in the G4RunManager
+  void Update_Scoring();
+
+};
+
+#endif
diff --git a/include/G4BeamTestTank.h b/include/G4BeamTestTank.h
new file mode 100644
index 0000000..b08e223
--- /dev/null
+++ b/include/G4BeamTestTank.h
@@ -0,0 +1,120 @@
+/**
+ * Copyright (C) 2009
+ * The IceCube collaboration
+ * ID: $Id: G4BeamTestTank.h 149431 2016-08-19 17:38:06Z jgonzalez $
+ *
+ * @file G4BeamTestTank.h
+ * @version $Rev: 149431 $
+ * @date $Date: 2016-08-19 18:38:06 +0100 (Fri, 19 Aug 2016) $
+ * @author Tilo Waldenmaier, Thomas Melzig
+ */
+
+
+#ifndef _TOPSIMULATOR_G4BEAMTESTTANK_H_
+#define _TOPSIMULATOR_G4BEAMTESTTANK_H_
+
+class OMKey;
+struct TankKey;
+class I3Geometry;
+class I3Position;
+
+class G4LogicalVolume;
+class G4VPhysicalVolume;
+class G4TankIceSD;
+
+#include <G4Types.hh>
+#include <G4ThreeVector.hh>
+#include <map>
+
+#include <icetray/I3Logging.h>
+
+/**
+ * This class constructs the physical volume of single tanks, serves as bridge to the sensitive detector in the tank, and computes the location of the point that define the snow surface. It also stores the basic tank dimensions (height, radius, snow).
+ *
+ * The sensitive detectors are not PMTs. Instead, the ice itself is the sensitive detector through G4TankIceSD. In this way, time is saved by not tracking photons.
+ *
+ * The points that are used to interpolate the snow surface are chosen so the snow surface is essentially flat on top of the tanks. This is done by determining the three vertices of a 10-meter equilateral triangle centered at the tank. The orientation of the triangle depends on whether the tank is part of a station with one or two tanks (if there are more tanks in a station, a fatal error occurs). In the single-tank case, the first vertex lies along the line pointing toward the first DOM. In the two-tank case, it lies along the line pointing toward it's neighboring tank.
+ *
+ */
+class G4BeamTestTank
+{
+ public:
+  G4BeamTestTank(const TankKey& tankKey, const I3Geometry& geometry);
+
+  ~G4BeamTestTank();
+
+  const G4ThreeVector& GetPos() const {return position_;}
+
+  /// Position of center of the tank
+  I3Position GetPos_I3();
+  double GetX_I3();
+  double GetY_I3();
+  double GetZ_I3();
+
+  /// Delaunay points around tank
+  const G4ThreeVector& GetDelaunayPoint1() const {return delaunayPoint1_;}
+  const G4ThreeVector& GetDelaunayPoint2() const {return delaunayPoint2_;}
+  const G4ThreeVector& GetDelaunayPoint3() const {return delaunayPoint3_;}
+
+  /// Tank dimensions in Geant4 units
+  G4double GetTankHeight_G4() {return tankHeight_;}
+  G4double GetTankRadius_G4() {return outerRadius_;}
+  G4double GetSnowHeight_G4() {return snowHeight_;}
+
+  /// Tank dimensions in IceCube units and reference system
+  double   GetTankHeight_I3();
+  double   GetTankRadius_I3();
+  double   GetSnowHeight_I3();
+
+  /// Energy deposit for a given OM, in Geant4 units (from G4TankIceSD)
+  G4double GetEDep_G4(const OMKey& omKey);
+  /// Emission time for a given OM, in Geant4 units (from G4TankIceSD)
+  G4double GetTime_G4(const OMKey& omKey);
+
+  /// Energy deposit for a given OM, in Geant4 units (from G4TankIceSD)
+  double   GetEDep_I3(const OMKey& omKey);
+  /// Emission time for a given OM, in Geant4 units (from G4TankIceSD)
+  double   GetTime_I3(const OMKey& omKey);
+
+  /// Number of Cherenkovs for a given OM (from G4TankIceSD)
+  double GetNumCherenkov(const OMKey& omKey);
+  /// Number of Cherenkovs weighted by emission point (from G4TankIceSD)
+  double GetNumCherenkovWeight(const OMKey& omKey);
+
+  const TankKey& GetTankKey() const {return tankKey_;}
+
+  /// Method to call when constructing detector
+  G4VPhysicalVolume* InstallTank(G4VPhysicalVolume* mother, const G4ThreeVector& origin);
+
+ private:
+
+  G4double tankThickness_;
+  G4double tankHeight_;
+  G4double fillHeight_;
+  G4double innerRadius_;
+  G4double outerRadius_;
+
+  G4double snowHeight_;
+  G4double perliteHeight_;
+
+  G4double glassOuterRadius_;
+  G4double glassThickness_;
+
+  G4ThreeVector position_;
+
+  G4ThreeVector delaunayPoint1_;
+  G4ThreeVector delaunayPoint2_;
+  G4ThreeVector delaunayPoint3_;
+
+  std::map<OMKey, G4ThreeVector> relDomPositions_;
+
+  G4LogicalVolume* tankLog_;
+  G4TankIceSD* iceSD_;
+
+  const TankKey& tankKey_;
+  const I3Geometry& geometry_;
+
+  SET_LOGGER("G4BeamTestTank");
+};
+
+#endif
diff --git a/include/G4BeamTestUserSteppingAction.h b/include/G4BeamTestUserSteppingAction.h
new file mode 100644
index 0000000..9807bbf
--- /dev/null
+++ b/include/G4BeamTestUserSteppingAction.h
@@ -0,0 +1,31 @@
+#ifndef G4BEAMTESTUSERSTEPPINGACTION_H_INCLUDED
+#define G4BEAMTESTUSERSTEPPINGACTION_H_INCLUDED
+/**
+ * Copyright (C) 2011
+ * The IceCube collaboration
+ * ID: $Id$
+ *
+ * @file G4BeamTestUserSteppingAction.h
+ * @version $Revision$
+ * @date $Date$
+ * @author Thomas Melzig
+ *
+ * $LastChangedBy$
+ */
+
+
+#include "G4UserSteppingAction.hh"
+
+/**
+ * Implementation of G4UserSteppingAction. This class kills gammas below threshold (set by G4BeamTestTank).
+ */
+class G4BeamTestUserSteppingAction : public G4UserSteppingAction {
+
+  public:
+    G4BeamTestUserSteppingAction();
+   ~G4BeamTestUserSteppingAction() {}
+
+    void UserSteppingAction(const G4Step*);
+};
+
+#endif  // G4BEAMTESTUSERSTEPPINGACTION_H_INCLUDED
diff --git a/include/G4BeamTestUserTrackingAction.h b/include/G4BeamTestUserTrackingAction.h
new file mode 100644
index 0000000..bb6d155
--- /dev/null
+++ b/include/G4BeamTestUserTrackingAction.h
@@ -0,0 +1,32 @@
+#ifndef G4BEAMTESTUSERTRACKINGACTION_H_INCLUDED
+#define G4BEAMTESTUSERTRACKINGACTION_H_INCLUDED
+/**
+ * Copyright (C) 2011
+ * The IceCube collaboration
+ * ID: $Id$
+ *
+ * @file G4BeamTestUserTrackingAction.h
+ * @version $Revision$
+ * @date $Date$
+ * @author Thomas Melzig
+ *
+ * $LastChangedBy$
+ */
+
+
+#include "G4UserTrackingAction.hh"
+
+/**
+ * Implementation of G4UserTrackingAction. This class kills gammas below threshold (set by G4BeamTestTank).
+ */
+class G4BeamTestUserTrackingAction : public G4UserTrackingAction {
+
+  public:
+    G4BeamTestUserTrackingAction();
+   ~G4BeamTestUserTrackingAction() {}
+
+    void PreUserTrackingAction(const G4Track*);
+    void PostUserTrackingAction(const G4Track*);
+};
+
+#endif  // G4BEAMTESTUSERTRACKINGACTION_H_INCLUDED
diff --git a/include/G4Interface.h b/include/G4Interface.h
new file mode 100644
index 0000000..c9bec07
--- /dev/null
+++ b/include/G4Interface.h
@@ -0,0 +1,68 @@
+/**
+ * Copyright (C) 2009
+ * The IceCube collaboration
+ * ID: $Id: G4Interface.h 149388 2016-08-18 21:50:04Z jgonzalez $
+ *
+ * @file G4Interface.h
+ * @version $Rev: 149388 $
+ * @date $Date: 2016-08-18 22:50:04 +0100 (Thu, 18 Aug 2016) $
+ * @author Tilo Waldenmaier
+ */
+
+#ifndef _TOPSIM_G4INTERFACE_H_
+#define _TOPSIM_G4INTERFACE_H_
+
+#include <g4-tankresponse/g4classes/G4BeamTestRunManager.h>
+#include <icetray/I3Logging.h>
+
+#ifdef G4VIS_USE
+class G4VisManager;
+#endif
+
+class I3Particle;
+class G4BeamTestTank;
+class G4BeamTestDetectorConstruction;
+
+/**
+ * Top-level class to handle Geant4. All global things are initialized here (run manager, visualization manager, detector construction, physics list and user actions).
+ */
+
+class G4Interface
+{
+ public:
+  G4Interface(const std::string& visMacro="");
+  ~G4Interface();
+
+  //  Static method which returns the singleton pointer to this class
+  static G4Interface* GetInstance() {return g4Interface_;}
+
+  /// Add a tank to the geometry. Should not be called after initialized.
+  void InstallTank(G4BeamTestTank* tank);
+
+  /// Initialize event. Most Geant4 global things are initialized the first time this is called.
+  void InitializeEvent();
+  /// To be called after simulating each IceTray event.
+  void TerminateEvent();
+  /// Simulate a single particle (InitializeEvent must be called first)
+  void InjectParticle(const I3Particle& particle);
+
+ private:
+  void Initialize();
+
+  static G4Interface* g4Interface_;
+
+  G4BeamTestRunManager runManager_;
+
+#ifdef G4VIS_USE
+  G4VisManager* visManager_;
+#endif
+
+  G4BeamTestDetectorConstruction* detector_;
+  bool initialized_;
+  bool eventInitialized_;
+  std::string visMacro_;
+
+  SET_LOGGER("G4Interface");
+};
+
+#endif
diff --git a/include/G4TankIceSD.h b/include/G4TankIceSD.h
new file mode 100644
index 0000000..75375c1
--- /dev/null
+++ b/include/G4TankIceSD.h
@@ -0,0 +1,69 @@
+/**
+ * Copyright (C) 2009
+ * The IceCube collaboration
+ * ID: $Id: G4TankIceSD.h 149388 2016-08-18 21:50:04Z jgonzalez $
+ *
+ * @file G4TankIceSD.h
+ * @version $Rev: 149388 $
+ * @date $Date: 2016-08-18 22:50:04 +0100 (Thu, 18 Aug 2016) $
+ * @author Tilo Waldenmaier, Thomas Melzig
+ */
+
+
+#ifndef _G4TANKRESPONSE_G4TankIceSD_H
+#define _G4TANKRESPONSE_G4TankIceSD_H
+
+#include <icetray/OMKey.h>
+
+#include <G4VSensitiveDetector.hh>
+#include <G4ThreeVector.hh>
+
+#include <map>
+
+class G4Step;
+class G4HCofThisEvent;
+class G4TouchableHistory;
+
+
+/**
+ * An "ice sensitive detector". This sensitive detector is meant to be associated with the ice logical volume in a tank.
+ *
+ * This class keeps track of the energy losses and number Cherenkov photons produced in the ice of each tank.
+ * The Cherenkov photons are counted in two ways. One is a simple count and the other is a weighted count
+ * where the weight depends on the distance from the photon emission point an the OM.
+ */
+class G4TankIceSD : public G4VSensitiveDetector
+{
+ public:
+   G4TankIceSD(G4String name, const std::map<OMKey, G4ThreeVector>& domPositions);
+  ~G4TankIceSD();
+
+  /// Methods called by Geant4 framework
+  void Initialize(G4HCofThisEvent *HCE);
+  G4bool ProcessHits(G4Step *aStep, G4TouchableHistory *ROhist);
+  void EndOfEvent(G4HCofThisEvent *HCE);
+
+  /// Get total energy deposit for a given OM (same for both OMs in a tank)
+  G4double GetEDep(const OMKey& omKey) {return sumEdep_[omKey];}
+  /// Get average emission time weighted by deposited energy (same for both OMs in a tank)
+  G4double GetTime(const OMKey& omKey) {return cogTime_[omKey];}
+  /// Get number of Cherenkov photons for a given OM (same for both OMs in a tank)
+  G4double GetNumCherenkov(const OMKey& omKey) {return cherenkovCounter_[omKey];}
+  /// Get number of Cherenkov photons for a given OM weighted relative to emission point
+  G4double GetNumCherenkovWeight(const OMKey& omKey) {return cherenkovCounterWeight_[omKey];}
+
+ private:
+  //ExN04TrackerHitsCollection *trackerCollection;
+  const std::map<OMKey, G4ThreeVector> domPositions_;
+
+  /// Cherenkov production. See technical note
+  G4double GetCerenkovNumber(G4Step* aStep);
+  G4double GetProbability(G4double radius, G4double height);
+
+  std::map<OMKey, G4double> sumEdep_;
+  std::map<OMKey, G4double> cogTime_;
+  std::map<OMKey, G4double> cherenkovCounter_;
+  std::map<OMKey, G4double> cherenkovCounterWeight_;
+};
+
+#endif
diff --git a/src/G4BeamTestDetectorConstruction.cxx b/src/G4BeamTestDetectorConstruction.cxx
new file mode 100644
index 0000000..6eeaa54
--- /dev/null
+++ b/src/G4BeamTestDetectorConstruction.cxx
@@ -0,0 +1,193 @@
+#include <G4NistManager.hh>
+#include <G4Material.hh>
+
+#include <G4LogicalVolume.hh>
+#include <G4PVPlacement.hh>
+
+#include <G4Box.hh>
+
+#include <G4VisAttributes.hh>
+#include <G4UserLimits.hh>
+
+#include "G4BeamTestDetectorConstruction.h"
+#include "G4BeamTestTank.h"
+
+
+G4BeamTestDetectorConstruction::G4BeamTestDetectorConstruction():
+  origin_(NAN, NAN, NAN), verboseLevel_(0)/* , tankList_(0) */
+{
+}
+
+
+G4BeamTestDetectorConstruction::~G4BeamTestDetectorConstruction()
+{
+}
+
+
+G4VPhysicalVolume* G4BeamTestDetectorConstruction::Construct()
+{
+  /* if(tankList_.empty()) return NULL; */
+
+  CreateMaterials();
+
+  // Determine bottom z position of snow volume
+  G4double zSnowBottom(tankList_.at(0)->GetPos().z());
+  BOOST_FOREACH(G4BeamTestTank* tank, tankList_)
+  {
+    // z position of bottom of tank
+    G4double z = tank->GetPos().z() - 0.5*tank->GetTankHeight_G4();
+    zSnowBottom = std::min(z, zSnowBottom);
+  }
+
+  // Subtract safety margin
+  zSnowBottom -= 1.0*CLHEP::m;
+
+  // Triangulate snow surface
+  G4Delaunay delaunay;
+  BOOST_FOREACH(G4BeamTestTank* tank, tankList_)
+  {
+    // z position of snow surface
+    G4double z = tank->GetPos().z() + 0.5 * tank->GetTankHeight_G4() + tank->GetSnowHeight_G4();
+
+    delaunay.AddPoint(tank->GetDelaunayPoint1().x(),
+                      tank->GetDelaunayPoint1().y(),
+                      z - zSnowBottom);
+    delaunay.AddPoint(tank->GetDelaunayPoint2().x(),
+                      tank->GetDelaunayPoint2().y(),
+                      z - zSnowBottom);
+    delaunay.AddPoint(tank->GetDelaunayPoint3().x(),
+                      tank->GetDelaunayPoint3().y(),
+                      z - zSnowBottom);
+  }
+
+  // Create tesselated snow volume
+  G4TessellatedSolid* solidSnow = new G4TessellatedSolid("solid_snow");
+  delaunay.BuildSolid(solidSnow, 50.0*CLHEP::m, 100.0*CLHEP::m);
+
+  // Determine World dimensions
+  G4double xHalfLength = 0.5 * (delaunay.GetXmax() - delaunay.GetXmin());
+  G4double yHalfLength = 0.5 * (delaunay.GetYmax() - delaunay.GetYmin());
+  G4double zHalfLength = 0.5 * (delaunay.GetZmax() + 20.0 * CLHEP::m); // 20 m of atmosphere
+
+  // World origin in IceCube coordinates
+  origin_.set(delaunay.GetOrigin().x(), delaunay.GetOrigin().y(), zSnowBottom + zHalfLength);
+
+  // Create world volume
+  G4Box* world_box = new G4Box("solid_world", xHalfLength, yHalfLength, zHalfLength);
+  G4LogicalVolume* worldLog =
+      new G4LogicalVolume(world_box, G4Material::GetMaterial("Air"), "log_world", 0, 0, 0);
+  G4VPhysicalVolume* worldPhys =
+      new G4PVPlacement(0, G4ThreeVector(), worldLog, "world", 0, false, 0);
+
+  // Snow layer
+  G4LogicalVolume* snowLog =
+      new G4LogicalVolume(solidSnow, G4Material::GetMaterial("Snow"), "log_snow", 0, 0, 0);
+  G4VPhysicalVolume* snowPhys =
+      new G4PVPlacement(0, G4ThreeVector(0, 0, -zHalfLength), snowLog, "snow", worldLog, false, 0);
+
+  // Instantiation of a set of visualization attributes with cyan colour
+  G4VisAttributes * snowVisAtt = new G4VisAttributes(G4Colour(0., 1., 1.));
+  // Assignment of the visualization attributes to the logical volume
+  snowLog->SetVisAttributes(snowVisAtt);
+
+  // Install tanks
+  BOOST_FOREACH(G4BeamTestTank* tank, tankList_)
+  {
+    tank->InstallTank(snowPhys, origin_);
+  }
+
+  // User limits (energy cutoffs)
+  // Do not create photons or electrons below cherenkov threshold
+  // See also corresponding UserSpecialCuts in Physicslist !!!!
+  G4UserLimits* energyLimit = new G4UserLimits();
+  energyLimit->SetUserMinEkine(280.0 * CLHEP::keV);  // Cherenkov threshold of electrons in ice TODO(shivesh)
+  worldLog->SetUserLimits(energyLimit);
+  snowLog->SetUserLimits(energyLimit);
+  
+  return worldPhys;
+}
+
+/*****************************************************************/
+
+void G4BeamTestDetectorConstruction::CreateMaterials()
+{
+  CreateAir();
+  /* CreateIce(); */
+  /* CreateSnow(); */
+  CreateWater();
+  CreatePlastic();
+  /* CreatePerlite(); */
+  CreateGlassSphere();
+  CreateEffectiveDOMMaterial();
+
+  //if(verboseLevel_>0) G4cout << *G4Material::GetMaterialTable() << G4endl;
+}
+
+/*****************************************************************/
+
+void G4BeamTestDetectorConstruction::CreateAir()
+{
+  G4NistManager* nistManager = G4NistManager::Instance();
+  nistManager->ConstructNewGasMaterial("Air","G4_AIR"/* , (273.15 - 40.0)*CLHEP::kelvin, 670.0E-3*CLHEP::bar */);
+}
+
+/*****************************************************************/
+
+void G4BeamTestDetectorConstruction::CreateWater()
+{
+  G4NistManager* nistManager = G4NistManager::Instance();
+  // G4Material* ice = new G4Material("Water", 1.0 * CLHEP::g / CLHEP::cm3, 2, kStateLiquid);
+  // ice->AddElement(nistManager->FindOrBuildElement("H"), 2);
+  // ice->AddElement(nistManager->FindOrBuildElement("O"), 1);
+  nistManager->ConstructNewGasMaterial("Water","G4_WATER");
+}
+
+/*****************************************************************/
+
+void G4BeamTestDetectorConstruction::CreatePlastic()
+{
+  G4NistManager* nistManager = G4NistManager::Instance();
+  /*
+  G4Material* plastic = new G4Material("Plastic", 1.5 * CLHEP::g / CLHEP::cm3, 3, kStateSolid);
+  plastic->AddElement(nistManager->FindOrBuildElement("H"), 20);
+  plastic->AddElement(nistManager->FindOrBuildElement("C"), 10);
+  plastic->AddElement(nistManager->FindOrBuildElement("O"), 5);
+  */
+  
+  // POM
+  G4Material* plastic = new G4Material("Plastic", 1.425 * CLHEP::g / CLHEP::cm3, 3, kStateSolid);
+  plastic->AddElement(nistManager->FindOrBuildElement("H"), 2);
+  plastic->AddElement(nistManager->FindOrBuildElement("C"), 1);
+  plastic->AddElement(nistManager->FindOrBuildElement("O"), 1);
+  
+  //nistManager->FindOrBuildMaterial("G4_POLYOXYMETHYLENE");
+  
+}
+
+/*****************************************************************/
+
+void G4BeamTestDetectorConstruction::CreateGlassSphere()
+{
+  // Elemental composition not exact for detailed line-up look at
+  // docushare collection 690 -> XRF and ICP Analysis of Pressure Sphere Glass
+
+  // 20 lbs. weight = 9072 g
+  // 6.5" outer radius & 0.5" thickness = 4024 cm3
+  G4NistManager* nistManager = G4NistManager::Instance();
+  G4Material* glass = new G4Material("Glass", 2.254 * CLHEP::g / CLHEP::cm3, 2, kStateSolid);
+  glass->AddElement(nistManager->FindOrBuildElement("Si"), 1);
+  glass->AddElement(nistManager->FindOrBuildElement("O"),  2);
+}
+
+/*****************************************************************/
+
+void G4BeamTestDetectorConstruction::CreateEffectiveDOMMaterial()
+{
+  // Mass of the a complete DOM: 12050 g
+  // Mass without glass sphere: 2978 g
+  // 6" inner glass radius: volume = 14830 cm3
+  G4NistManager* nistManager = G4NistManager::Instance();
+  G4Material* glass = new G4Material("effectiveDOM", 0.2 * CLHEP::g / CLHEP::cm3, 2, kStateSolid);
+  glass->AddElement(nistManager->FindOrBuildElement("Si"), 1);
+  glass->AddElement(nistManager->FindOrBuildElement("O"),  2);
+}
diff --git a/src/G4BeamTestEMPhysics.cxx b/src/G4BeamTestEMPhysics.cxx
new file mode 100644
index 0000000..d36207a
--- /dev/null
+++ b/src/G4BeamTestEMPhysics.cxx
@@ -0,0 +1,71 @@
+/*
+ * copyright  (C) 2010
+ * The Icecube Collaboration
+ *
+ * $Id: G4BeamTestEMPhysics.cxx 86420 2012-03-20 16:00:37Z kislat $
+ *
+ * @version $Revision: 86420 $
+ * @date $LastChangedDate: 2012-03-20 16:00:37 +0000 (Tue, 20 Mar 2012) $
+ * @author Fabian Kislat <fabian.kislat@desy.de>, Peter Nießen (tanktop)
+ * Last changed by: $LastChangedBy: kislat $
+ */
+
+
+#include <globals.hh>
+#include <G4ParticleDefinition.hh>
+#include <G4ParticleTable.hh>
+#include <G4Gamma.hh>
+#include <G4Electron.hh>
+#include <G4Positron.hh>
+#include <G4NeutrinoE.hh>
+#include <G4AntiNeutrinoE.hh>
+#include <G4ProcessManager.hh>
+
+#include "G4BeamTestEMPhysics.h"
+
+
+G4BeamTestEMPhysics::G4BeamTestEMPhysics()
+  : G4VPhysicsConstructor("standard")
+{}
+
+
+G4BeamTestEMPhysics::~G4BeamTestEMPhysics()
+{}
+
+
+void G4BeamTestEMPhysics::ConstructParticle()
+{
+  // gamma
+  G4Gamma::GammaDefinition();
+ 
+  // electron
+  G4Electron::ElectronDefinition();
+  G4Positron::PositronDefinition();
+  G4NeutrinoE::NeutrinoEDefinition();
+  G4AntiNeutrinoE::AntiNeutrinoEDefinition();
+}
+
+
+void G4BeamTestEMPhysics::ConstructProcess()
+{
+  G4ProcessManager *pManager = 0;
+  
+  // Gamma Physics
+  pManager = G4Gamma::Gamma()->GetProcessManager();
+  pManager->AddDiscreteProcess(&photoEffect);
+  pManager->AddDiscreteProcess(&comptonEffect);
+  pManager->AddDiscreteProcess(&pairProduction);
+
+  // Electron Physics
+  pManager = G4Electron::Electron()->GetProcessManager();
+  pManager->AddProcess(&electronMultipleScattering, -1, 1, 1);
+  pManager->AddProcess(&electronIonisation,         -1, 2, 2);
+  pManager->AddProcess(&electronBremsStrahlung,     -1, 3, 3);  
+
+  //Positron Physics
+  pManager = G4Positron::Positron()->GetProcessManager();
+  pManager->AddProcess(&positronMultipleScattering, -1, 1, 1);
+  pManager->AddProcess(&positronIonisation,         -1, 2, 2);
+  pManager->AddProcess(&positronBremsStrahlung,     -1, 3, 3);  
+  pManager->AddProcess(&annihilation,                0,-1, 4);  
+}
diff --git a/src/G4BeamTestGeneralPhysics.cxx b/src/G4BeamTestGeneralPhysics.cxx
new file mode 100644
index 0000000..2def9b4
--- /dev/null
+++ b/src/G4BeamTestGeneralPhysics.cxx
@@ -0,0 +1,62 @@
+/*
+ * copyright  (C) 2010
+ * The Icecube Collaboration
+ *
+ * $Id: G4BeamTestGeneralPhysics.cxx 152849 2017-01-20 21:44:25Z jgonzalez $
+ *
+ * @version $Revision: 152849 $
+ * @date $LastChangedDate: 2017-01-20 21:44:25 +0000 (Fri, 20 Jan 2017) $
+ * @author Fabian Kislat <fabian.kislat@desy.de> Last changed by: $LastChangedBy: jgonzalez $
+ */
+
+#include "G4BeamTestGeneralPhysics.h"
+
+#include <iomanip>   
+
+#include <globals.hh>
+#include <G4ios.hh>
+#include <G4ParticleDefinition.hh>
+#include <G4ProcessManager.hh>
+// Bosons
+#include <G4ChargedGeantino.hh>
+#include <G4Geantino.hh>
+#include <G4Version.hh>
+
+G4BeamTestGeneralPhysics::G4BeamTestGeneralPhysics()
+  :  G4VPhysicsConstructor("general")
+{}
+
+
+G4BeamTestGeneralPhysics::~G4BeamTestGeneralPhysics()
+{}
+
+
+void G4BeamTestGeneralPhysics::ConstructParticle()
+{
+  // pseudo-particles
+  G4Geantino::GeantinoDefinition();
+  G4ChargedGeantino::ChargedGeantinoDefinition();  
+}
+
+
+void G4BeamTestGeneralPhysics::ConstructProcess()
+{
+  //AddTransportation();
+
+// Decay processes are set somewhere else now
+#if G4VERSION_NUMBER < 1000
+  // Add Decay Process
+  theParticleIterator->reset();
+  while ((*theParticleIterator)()) {
+    G4ParticleDefinition *particle = theParticleIterator->value();
+    G4ProcessManager *pmanager = particle->GetProcessManager();
+    if (decay.IsApplicable(*particle)) { 
+      pmanager->AddProcess(&decay);
+      // set ordering for PostStepDoIt and AtRestDoIt
+      pmanager->SetProcessOrdering(&decay, idxPostStep);
+      pmanager->SetProcessOrdering(&decay, idxAtRest);
+    }
+  }
+#endif
+}
+
diff --git a/src/G4BeamTestHadronPhysics.cxx b/src/G4BeamTestHadronPhysics.cxx
new file mode 100644
index 0000000..5150b67
--- /dev/null
+++ b/src/G4BeamTestHadronPhysics.cxx
@@ -0,0 +1,425 @@
+/*
+ * copyright  (C) 2010
+ * The Icecube Collaboration
+ *
+ * $Id: G4BeamTestHadronPhysics.cxx 152814 2017-01-19 21:34:52Z jgonzalez $
+ *
+ * @version $Revision: 152814 $
+ * @date $LastChangedDate: 2017-01-19 21:34:52 +0000 (Thu, 19 Jan 2017) $
+ * @author Fabian Kislat <fabian.kislat@desy.de> Last changed by: $LastChangedBy: jgonzalez $
+ *
+ * Copied from tanktop
+ */
+
+#include <iomanip>   
+
+#include <globals.hh>
+#include <G4ios.hh>
+#include <G4ParticleDefinition.hh>
+#include <G4ParticleTable.hh>
+#include <G4MesonConstructor.hh>
+#include <G4BaryonConstructor.hh>
+#include <G4ShortLivedConstructor.hh>
+#include <G4ProcessManager.hh>
+
+#include "G4BeamTestHadronPhysics.h"
+
+G4BeamTestHadronPhysics::G4BeamTestHadronPhysics()
+  :  G4VPhysicsConstructor("hadron") {
+}
+
+/********************************************************************/
+
+G4BeamTestHadronPhysics::~G4BeamTestHadronPhysics()
+{
+  delete stringDecay_;
+}
+
+/********************************************************************/
+
+void G4BeamTestHadronPhysics::ConstructParticle()
+{
+
+  //  Construct all mesons
+  G4MesonConstructor mesonConstructor;
+  mesonConstructor.ConstructParticle();
+
+  //  Construct all barions
+  G4BaryonConstructor baryonConstructor;
+  baryonConstructor.ConstructParticle();
+
+  //  Construct resonaces and quarks
+  G4ShortLivedConstructor shortLivedConstructor;
+  shortLivedConstructor.ConstructParticle();  
+
+}
+
+/********************************************************************/
+
+void G4BeamTestHadronPhysics::ConstructProcess()
+{
+  G4ProcessManager *pManager = 0;
+
+  // G4cout << "" << G4endl;
+  // G4cout << "You are using the G4BeamTestHadronPhysics" << G4endl;
+  // G4cout << " - Note that this hadronic physics list is not optimized for any particular usage" << G4endl;
+  // G4cout << " - If you wish to have a starting point tailored for a particular area of work," << G4endl;
+  // G4cout << "   please use one of the available physics lists by use-case." << G4endl;
+  // G4cout << "   More information can also be found from the Geant4 HyperNews." << G4endl;
+  // G4cout << "" << G4endl;
+
+  // Elastic Process
+  elasticModel_ = new G4LElastic();
+  elasticProcess_.RegisterMe(elasticModel_);
+
+  // pi+ and pi-
+  preEquilib_ = new G4PreCompoundModel(&excitationHandler_);
+  cascade_.SetDeExcitation(preEquilib_);  
+  theoModel_.SetTransport(&cascade_);
+  theoModel_.SetHighEnergyGenerator(&stringModel_);
+  stringDecay_ = new G4ExcitedStringDecay(&fragmentation_);
+  stringModel_.SetFragmentationModel(stringDecay_);
+  theoModel_.SetMinEnergy(15*CLHEP::GeV);
+  theoModel_.SetMaxEnergy(100*CLHEP::TeV);
+
+  // PionPlus
+  pManager = G4PionPlus::PionPlus()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  lePionPlusModel_ = new G4LEPionPlusInelastic();
+  pionPlusInelastic_.RegisterMe(lePionPlusModel_);
+  pionPlusInelastic_.RegisterMe(&theoModel_);
+  pManager->AddDiscreteProcess(&pionPlusInelastic_);
+
+  pManager->AddProcess(&pionPlusIonisation_, ordInActive,2, 2);
+
+  pManager->AddProcess(&pionPlusMult_);
+  pManager->SetProcessOrdering(&pionPlusMult_, idxAlongStep, 1);
+  pManager->SetProcessOrdering(&pionPlusMult_, idxPostStep, 1);
+
+  // PionMinus
+  pManager = G4PionMinus::PionMinus()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  lePionMinusModel_ = new G4LEPionMinusInelastic();
+  pionMinusInelastic_.RegisterMe(lePionMinusModel_);
+  pionMinusInelastic_.RegisterMe(&theoModel_);
+  pManager->AddDiscreteProcess(&pionMinusInelastic_);
+
+  pManager->AddProcess(&pionMinusIonisation_, ordInActive,2, 2);
+
+  pManager->AddProcess(&pionMinusMult_);
+  pManager->SetProcessOrdering(&pionMinusMult_, idxAlongStep, 1);
+  pManager->SetProcessOrdering(&pionMinusMult_, idxPostStep, 1);
+
+  pManager->AddRestProcess(&pionMinusAbsorption_, ordDefault);
+
+  // KaonPlus
+  pManager = G4KaonPlus::KaonPlus()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  leKaonPlusModel_ = new G4LEKaonPlusInelastic();
+  heKaonPlusModel_ = new G4HEKaonPlusInelastic();
+  kaonPlusInelastic_.RegisterMe(leKaonPlusModel_);
+  kaonPlusInelastic_.RegisterMe(&theoModel_);
+  pManager->AddDiscreteProcess(&kaonPlusInelastic_);
+
+  pManager->AddProcess(&kaonPlusIonisation_, ordInActive,2, 2);
+
+  pManager->AddProcess(&kaonPlusMult_);
+  pManager->SetProcessOrdering(&kaonPlusMult_, idxAlongStep, 1);
+  pManager->SetProcessOrdering(&kaonPlusMult_, idxPostStep, 1);
+
+  // KaonMinus
+  pManager = G4KaonMinus::KaonMinus()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  leKaonMinusModel_ = new G4LEKaonMinusInelastic();
+  heKaonMinusModel_ = new G4HEKaonMinusInelastic();
+  kaonMinusInelastic_.RegisterMe(leKaonMinusModel_);
+  kaonMinusInelastic_.RegisterMe(heKaonMinusModel_);
+  pManager->AddDiscreteProcess(&kaonMinusInelastic_);
+
+  pManager->AddProcess(&kaonMinusIonisation_, ordInActive,2, 2);
+
+  pManager->AddProcess(&kaonMinusMult_);
+  pManager->SetProcessOrdering(&kaonMinusMult_, idxAlongStep, 1);
+  pManager->SetProcessOrdering(&kaonMinusMult_, idxPostStep, 1);
+
+  pManager->AddRestProcess(&kaonMinusAbsorption_, ordDefault);
+
+  // KaonZeroL
+  pManager = G4KaonZeroLong::KaonZeroLong()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  leKaonZeroLModel_ = new G4LEKaonZeroLInelastic();
+  heKaonZeroLModel_ = new G4HEKaonZeroInelastic();
+  kaonZeroLInelastic_.RegisterMe(leKaonZeroLModel_);
+  kaonZeroLInelastic_.RegisterMe(heKaonZeroLModel_);
+  pManager->AddDiscreteProcess(&kaonZeroLInelastic_);
+ 
+  // KaonZeroS
+  pManager = G4KaonZeroShort::KaonZeroShort()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  leKaonZeroSModel_ = new G4LEKaonZeroSInelastic();
+  heKaonZeroSModel_ = new G4HEKaonZeroInelastic();
+  kaonZeroSInelastic_.RegisterMe(leKaonZeroSModel_);
+  kaonZeroSInelastic_.RegisterMe(heKaonZeroSModel_);
+  pManager->AddDiscreteProcess(&kaonZeroSInelastic_);
+
+  // Proton
+  pManager = G4Proton::Proton()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  leProtonModel_ = new G4LEProtonInelastic();
+  heProtonModel_ = new G4HEProtonInelastic();
+  protonInelastic_.RegisterMe(leProtonModel_);
+  protonInelastic_.RegisterMe(&theoModel_);
+  pManager->AddDiscreteProcess(&protonInelastic_);
+
+  pManager->AddProcess(&protonIonisation_, ordInActive,2, 2);
+
+  pManager->AddProcess(&protonMult_);
+  pManager->SetProcessOrdering(&protonMult_, idxAlongStep, 1);
+  pManager->SetProcessOrdering(&protonMult_, idxPostStep, 1);
+
+  // anti-Proton
+  pManager = G4AntiProton::AntiProton()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  leAntiProtonModel_ = new G4LEAntiProtonInelastic();
+  heAntiProtonModel_ = new G4HEAntiProtonInelastic();
+  antiProtonInelastic_.RegisterMe(leAntiProtonModel_);
+  antiProtonInelastic_.RegisterMe(heAntiProtonModel_);
+  pManager->AddDiscreteProcess(&antiProtonInelastic_);
+
+  pManager->AddProcess(&antiProtonIonisation_, ordInActive,2, 2);
+
+  pManager->AddProcess(&antiProtonMult_);
+  pManager->SetProcessOrdering(&antiProtonMult_, idxAlongStep, 1);
+  pManager->SetProcessOrdering(&antiProtonMult_, idxPostStep, 1);
+
+  pManager->AddRestProcess(&antiProtonAnnihilation_);
+
+  // Neutron
+  pManager = G4Neutron::Neutron()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  leNeutronModel_ = new G4LENeutronInelastic();
+  heNeutronModel_ = new G4HENeutronInelastic();
+  neutronInelastic_.RegisterMe(leNeutronModel_);
+  neutronInelastic_.RegisterMe(&theoModel_);
+  pManager->AddDiscreteProcess(&neutronInelastic_);
+  
+  neutronFissionModel_ = new G4LFission();
+  neutronFission_.RegisterMe(neutronFissionModel_);
+  pManager->AddDiscreteProcess(&neutronFission_);
+
+  neutronCaptureModel_ = new G4LCapture();
+  neutronCapture_.RegisterMe(neutronCaptureModel_);
+  pManager->AddDiscreteProcess(&neutronCapture_);
+
+  // AntiNeutron
+  pManager = G4AntiNeutron::AntiNeutron()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  leAntiNeutronModel_ = new G4LEAntiNeutronInelastic();
+  heAntiNeutronModel_ = new G4HEAntiNeutronInelastic();
+  antiNeutronInelastic_.RegisterMe(leAntiNeutronModel_);
+  antiNeutronInelastic_.RegisterMe(heAntiNeutronModel_);
+  pManager->AddDiscreteProcess(&antiNeutronInelastic_);
+    
+  pManager->AddRestProcess(&antiNeutronAnnihilation_);
+
+  // Lambda
+  pManager = G4Lambda::Lambda()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  leLambdaModel_ = new G4LELambdaInelastic();
+  heLambdaModel_ = new G4HELambdaInelastic();
+  lambdaInelastic_.RegisterMe(leLambdaModel_);
+  lambdaInelastic_.RegisterMe(heLambdaModel_);
+  pManager->AddDiscreteProcess(&lambdaInelastic_);
+  
+  // AntiLambda
+  pManager = G4AntiLambda::AntiLambda()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  leAntiLambdaModel_ = new G4LEAntiLambdaInelastic();
+  heAntiLambdaModel_ = new G4HEAntiLambdaInelastic();
+  antiLambdaInelastic_.RegisterMe(leAntiLambdaModel_);
+  antiLambdaInelastic_.RegisterMe(heAntiLambdaModel_);
+  pManager->AddDiscreteProcess(&antiLambdaInelastic_);
+    
+  // SigmaMinus
+  pManager = G4SigmaMinus::SigmaMinus()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  leSigmaMinusModel_ = new G4LESigmaMinusInelastic();
+  heSigmaMinusModel_ = new G4HESigmaMinusInelastic();
+  sigmaMinusInelastic_.RegisterMe(leSigmaMinusModel_);
+  sigmaMinusInelastic_.RegisterMe(heSigmaMinusModel_);
+  pManager->AddDiscreteProcess(&sigmaMinusInelastic_);
+
+  pManager->AddProcess(&sigmaMinusIonisation_, ordInActive,2, 2);
+
+  pManager->AddProcess(&sigmaMinusMult_);
+  pManager->SetProcessOrdering(&sigmaMinusMult_, idxAlongStep, 1);
+  pManager->SetProcessOrdering(&sigmaMinusMult_, idxPostStep, 1);
+
+  // anti-SigmaMinus
+  pManager = G4AntiSigmaMinus::AntiSigmaMinus()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  leAntiSigmaMinusModel_ = new G4LEAntiSigmaMinusInelastic();
+  heAntiSigmaMinusModel_ = new G4HEAntiSigmaMinusInelastic();
+  antiSigmaMinusInelastic_.RegisterMe(leAntiSigmaMinusModel_);
+  antiSigmaMinusInelastic_.RegisterMe(heAntiSigmaMinusModel_);
+  pManager->AddDiscreteProcess(&antiSigmaMinusInelastic_);
+
+  pManager->AddProcess(&antiSigmaMinusIonisation_, ordInActive,2, 2);
+
+  pManager->AddProcess(&antiSigmaMinusMult_);
+  pManager->SetProcessOrdering(&antiSigmaMinusMult_, idxAlongStep, 1);
+  pManager->SetProcessOrdering(&antiSigmaMinusMult_, idxPostStep, 1);
+
+  // SigmaPlus
+  pManager = G4SigmaPlus::SigmaPlus()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  leSigmaPlusModel_ = new G4LESigmaPlusInelastic();
+  heSigmaPlusModel_ = new G4HESigmaPlusInelastic();
+  sigmaPlusInelastic_.RegisterMe(leSigmaPlusModel_);
+  sigmaPlusInelastic_.RegisterMe(heSigmaPlusModel_);
+  pManager->AddDiscreteProcess(&sigmaPlusInelastic_);
+
+  pManager->AddProcess(&sigmaPlusIonisation_, ordInActive,2, 2);
+
+  pManager->AddProcess(&sigmaPlusMult_);
+  pManager->SetProcessOrdering(&sigmaPlusMult_, idxAlongStep, 1);
+  pManager->SetProcessOrdering(&sigmaPlusMult_, idxPostStep, 1);
+
+  // anti-SigmaPlus
+  pManager = G4AntiSigmaPlus::AntiSigmaPlus()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  leAntiSigmaPlusModel_ = new G4LEAntiSigmaPlusInelastic();
+  heAntiSigmaPlusModel_ = new G4HEAntiSigmaPlusInelastic();
+  antiSigmaPlusInelastic_.RegisterMe(leAntiSigmaPlusModel_);
+  antiSigmaPlusInelastic_.RegisterMe(heAntiSigmaPlusModel_);
+  pManager->AddDiscreteProcess(&antiSigmaPlusInelastic_);
+
+  pManager->AddProcess(&antiSigmaPlusIonisation_, ordInActive,2, 2);
+
+  pManager->AddProcess(&antiSigmaPlusMult_);
+  pManager->SetProcessOrdering(&antiSigmaPlusMult_, idxAlongStep, 1);
+  pManager->SetProcessOrdering(&antiSigmaPlusMult_, idxPostStep, 1);
+
+  // XiMinus
+  pManager = G4XiMinus::XiMinus()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  leXiMinusModel_ = new G4LEXiMinusInelastic();
+  heXiMinusModel_ = new G4HEXiMinusInelastic();
+  xiMinusInelastic_.RegisterMe(leXiMinusModel_);
+  xiMinusInelastic_.RegisterMe(heXiMinusModel_);
+  pManager->AddDiscreteProcess(&xiMinusInelastic_);
+
+  pManager->AddProcess(&xiMinusIonisation_, ordInActive,2, 2);
+
+  pManager->AddProcess(&xiMinusMult_);
+  pManager->SetProcessOrdering(&xiMinusMult_, idxAlongStep, 1);
+  pManager->SetProcessOrdering(&xiMinusMult_, idxPostStep, 1);
+
+  // anti-XiMinus
+  pManager = G4AntiXiMinus::AntiXiMinus()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  leAntiXiMinusModel_ = new G4LEAntiXiMinusInelastic();
+  heAntiXiMinusModel_ = new G4HEAntiXiMinusInelastic();
+  antiXiMinusInelastic_.RegisterMe(leAntiXiMinusModel_);
+  antiXiMinusInelastic_.RegisterMe(heAntiXiMinusModel_);
+  pManager->AddDiscreteProcess(&antiXiMinusInelastic_);
+
+  pManager->AddProcess(&antiXiMinusIonisation_, ordInActive,2, 2);
+
+  pManager->AddProcess(&antiXiMinusMult_);
+  pManager->SetProcessOrdering(&antiXiMinusMult_, idxAlongStep, 1);
+  pManager->SetProcessOrdering(&antiXiMinusMult_, idxPostStep, 1);
+
+  // XiZero
+  pManager = G4XiZero::XiZero()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  leXiZeroModel_ = new G4LEXiZeroInelastic();
+  heXiZeroModel_ = new G4HEXiZeroInelastic();
+  xiZeroInelastic_.RegisterMe(leXiZeroModel_);
+  xiZeroInelastic_.RegisterMe(heXiZeroModel_);
+  pManager->AddDiscreteProcess(&xiZeroInelastic_);
+
+  // anti-XiZero
+  pManager = G4AntiXiZero::AntiXiZero()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  leAntiXiZeroModel_ = new G4LEAntiXiZeroInelastic();
+  heAntiXiZeroModel_ = new G4HEAntiXiZeroInelastic();
+  antiXiZeroInelastic_.RegisterMe(leAntiXiZeroModel_);
+  antiXiZeroInelastic_.RegisterMe(heAntiXiZeroModel_);
+  pManager->AddDiscreteProcess(&antiXiZeroInelastic_);
+
+  // OmegaMinus
+  pManager = G4OmegaMinus::OmegaMinus()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  leOmegaMinusModel_ = new G4LEOmegaMinusInelastic();
+  heOmegaMinusModel_ = new G4HEOmegaMinusInelastic();
+  omegaMinusInelastic_.RegisterMe(leOmegaMinusModel_);
+  omegaMinusInelastic_.RegisterMe(heOmegaMinusModel_);
+  pManager->AddDiscreteProcess(&omegaMinusInelastic_);
+
+  pManager->AddProcess(&omegaMinusIonisation_, ordInActive,2, 2);
+
+  pManager->AddProcess(&omegaMinusMult_);
+  pManager->SetProcessOrdering(&omegaMinusMult_, idxAlongStep, 1);
+  pManager->SetProcessOrdering(&omegaMinusMult_, idxPostStep, 1);
+
+  // anti-OmegaMinus
+  pManager = G4AntiOmegaMinus::AntiOmegaMinus()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  leAntiOmegaMinusModel_ = new G4LEAntiOmegaMinusInelastic();
+  heAntiOmegaMinusModel_ = new G4HEAntiOmegaMinusInelastic();
+  antiOmegaMinusInelastic_.RegisterMe(leAntiOmegaMinusModel_);
+  antiOmegaMinusInelastic_.RegisterMe(heAntiOmegaMinusModel_);
+  pManager->AddDiscreteProcess(&antiOmegaMinusInelastic_);
+
+  pManager->AddProcess(&antiOmegaMinusIonisation_, ordInActive,2, 2);
+
+  pManager->AddProcess(&antiOmegaMinusMult_);
+  pManager->SetProcessOrdering(&antiOmegaMinusMult_, idxAlongStep, 1);
+  pManager->SetProcessOrdering(&antiOmegaMinusMult_, idxPostStep, 1);
+
+}
diff --git a/src/G4BeamTestIonPhysics.cxx b/src/G4BeamTestIonPhysics.cxx
new file mode 100644
index 0000000..3c7e553
--- /dev/null
+++ b/src/G4BeamTestIonPhysics.cxx
@@ -0,0 +1,112 @@
+/*
+ * copyright  (C) 2010
+ * The Icecube Collaboration
+ *
+ * $Id: G4BeamTestIonPhysics.cxx 86420 2012-03-20 16:00:37Z kislat $
+ *
+ * @version $Revision: 86420 $
+ * @date $LastChangedDate: 2012-03-20 16:00:37 +0000 (Tue, 20 Mar 2012) $
+ * @author Fabian Kislat <fabian.kislat@desy.de>, Peter Nießen (tanktop)
+ * Last changed by: $LastChangedBy: kislat $
+ */
+
+#include <iomanip>   
+
+#include <globals.hh>
+#include <G4ios.hh>
+#include <G4ParticleDefinition.hh>
+#include <G4ParticleTable.hh>
+#include <G4IonConstructor.hh>
+#include <G4ProcessManager.hh>
+
+#include "G4BeamTestIonPhysics.h"
+
+G4BeamTestIonPhysics::G4BeamTestIonPhysics()
+  :  G4VPhysicsConstructor ("ion")
+{}
+
+/********************************************************************/
+
+G4BeamTestIonPhysics::~G4BeamTestIonPhysics()
+{}
+
+
+/********************************************************************/
+
+void G4BeamTestIonPhysics::ConstructParticle()
+{
+  //  Construct light ions
+  G4IonConstructor pConstructor;
+  pConstructor.ConstructParticle();
+}
+
+/********************************************************************/
+
+
+
+
+void G4BeamTestIonPhysics::ConstructProcess()
+{
+  G4ProcessManager *pManager = 0;
+  
+  // Elastic Process
+  elasticModel_ = new G4LElastic;
+  elasticProcess_.RegisterMe(elasticModel_);
+
+  // Generic Ion
+  pManager = G4GenericIon::GenericIon()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  pManager->AddProcess(&ionMultipleScattering_, -1, 1, 1);
+  pManager->AddProcess(&ionIonisation_,         -1, 2, 2);
+
+  // Deuteron 
+  pManager = G4Deuteron::Deuteron()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  deuteronModel_ = new G4LEDeuteronInelastic;
+  deuteronProcess_.RegisterMe(deuteronModel_);
+  pManager->AddDiscreteProcess(&deuteronProcess_);
+
+  pManager->AddProcess(&deuteronMultipleScattering_, -1, 1, 1);
+  pManager->AddProcess(&deuteronIonisation_,        -1, 2, 2);
+
+  // Triton 
+  pManager = G4Triton::Triton()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  tritonModel_ = new G4LETritonInelastic;
+  tritonProcess_.RegisterMe(tritonModel_);
+  pManager->AddDiscreteProcess(&tritonProcess_);
+
+  pManager->AddProcess(&tritonMultipleScattering_, -1, 1, 1);
+  pManager->AddProcess(&tritonIonisation_,         -1, 2, 2);
+ 
+  // Alpha 
+  pManager = G4Alpha::Alpha()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  alphaModel_ = new G4LEAlphaInelastic;
+  alphaProcess_.RegisterMe(alphaModel_);
+  pManager->AddDiscreteProcess(&alphaProcess_);
+
+  pManager->AddProcess(&alphaMultipleScattering_, -1, 1, 1);
+  pManager->AddProcess(&alphaIonisation_,         -1, 2, 2);
+
+  // He3
+  pManager = G4He3::He3()->GetProcessManager();
+  // add process
+  pManager->AddDiscreteProcess(&elasticProcess_);
+
+  pManager->AddProcess(&he3MultipleScattering_, -1, 1, 1);
+  pManager->AddProcess(&he3Ionisation_,         -1, 2, 2);
+
+   
+}
+
+
+
diff --git a/src/G4BeamTestMuonPhysics.cxx b/src/G4BeamTestMuonPhysics.cxx
new file mode 100644
index 0000000..92264a2
--- /dev/null
+++ b/src/G4BeamTestMuonPhysics.cxx
@@ -0,0 +1,94 @@
+/*
+ * copyright  (C) 2010
+ * The Icecube Collaboration
+ *
+ * $Id: G4BeamTestMuonPhysics.cxx 86420 2012-03-20 16:00:37Z kislat $
+ *
+ * @version $Revision: 86420 $
+ * @date $LastChangedDate: 2012-03-20 16:00:37 +0000 (Tue, 20 Mar 2012) $
+ * @author Fabian Kislat <fabian.kislat@desy.de>, Peter Nießen (tanktop)
+ * Last changed by: $LastChangedBy: kislat $
+ */
+
+
+#include <globals.hh>
+#include <G4ParticleDefinition.hh>
+#include <G4ParticleTable.hh>
+#include <G4MuonPlus.hh>
+#include <G4MuonMinus.hh>
+#include <G4TauMinus.hh>
+#include <G4TauPlus.hh>
+#include <G4NeutrinoTau.hh>
+#include <G4AntiNeutrinoTau.hh>
+#include <G4NeutrinoMu.hh>
+#include <G4AntiNeutrinoMu.hh>
+#include <G4ProcessManager.hh>
+
+#include "G4BeamTestMuonPhysics.h"
+
+
+
+G4BeamTestMuonPhysics::G4BeamTestMuonPhysics()
+  :  G4VPhysicsConstructor("muon")
+{}
+
+/********************************************************************/
+
+G4BeamTestMuonPhysics::~G4BeamTestMuonPhysics()
+{}
+
+/********************************************************************/
+
+void G4BeamTestMuonPhysics::ConstructParticle()
+{
+  // Mu
+  G4MuonPlus::MuonPlusDefinition();
+  G4MuonMinus::MuonMinusDefinition();
+  G4NeutrinoMu::NeutrinoMuDefinition();
+  G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
+
+  // Tau
+  G4TauMinus::TauMinusDefinition();
+  G4TauPlus::TauPlusDefinition();
+  G4NeutrinoTau::NeutrinoTauDefinition();
+  G4AntiNeutrinoTau::AntiNeutrinoTauDefinition();
+}
+
+/********************************************************************/
+
+
+void G4BeamTestMuonPhysics::ConstructProcess()
+{
+
+  G4ProcessManager *pManager = 0;
+
+  // Muon Plus Physics
+  pManager = G4MuonPlus::MuonPlus()->GetProcessManager();
+   
+  pManager->AddProcess(&muPlusMultipleScattering_, -1, 1, 1);
+  pManager->AddProcess(&muPlusIonisation_,         -1, 2, 2);
+  pManager->AddProcess(&muPlusBremsstrahlung_,     -1, 3, 3);
+  pManager->AddProcess(&muPlusPairProduction_,     -1, 4, 4);
+
+  // Muon Minus Physics
+  pManager = G4MuonMinus::MuonMinus()->GetProcessManager();
+   
+  pManager->AddProcess(&muMinusMultipleScattering_, -1,  1, 1);
+  pManager->AddProcess(&muMinusIonisation_,         -1,  2, 2);
+  pManager->AddProcess(&muMinusBremsstrahlung_,     -1,  3, 3);
+  pManager->AddProcess(&muMinusPairProduction_,     -1,  4, 4);
+
+  pManager->AddRestProcess(&muMinusCaptureAtRest_);
+
+  // Tau Plus Physics
+  pManager = G4TauPlus::TauPlus()->GetProcessManager();
+   
+  pManager->AddProcess(&tauPlusMultipleScattering_, -1, 1, 1);
+  pManager->AddProcess(&tauPlusIonisation_,         -1, 2, 2);
+
+  // Tau Minus Physics
+  pManager = G4TauMinus::TauMinus()->GetProcessManager();
+
+  pManager->AddProcess(&tauMinusMultipleScattering_, -1, 1, 1);
+  pManager->AddProcess(&tauMinusIonisation_,         -1, 2, 2);
+}
diff --git a/src/G4BeamTestPhysicsList.cxx b/src/G4BeamTestPhysicsList.cxx
new file mode 100644
index 0000000..5a5befd
--- /dev/null
+++ b/src/G4BeamTestPhysicsList.cxx
@@ -0,0 +1,81 @@
+/**
+ * Copyright (C) 2009
+ * The IceCube collaboration
+ * ID: $Id: G4BeamTestPhysicsList.cxx 152849 2017-01-20 21:44:25Z jgonzalez $
+ *
+ * @file G4BeamTestPhysicsList.cxx
+ * @version $Rev: 152849 $
+ * @date $Date: 2017-01-20 21:44:25 +0000 (Fri, 20 Jan 2017) $
+ * @author Tilo Waldenmaier, Thomas Melzig, Javier Gonzalez
+ */
+
+
+#include <globals.hh>
+#include <G4Version.hh>
+#include "G4BeamTestPhysicsList.h"
+#include "G4BeamTestGeneralPhysics.h"
+#if G4VERSION_NUMBER < 1000
+#include "G4BeamTestEMPhysics.h"
+#include "G4BeamTestMuonPhysics.h"
+#include "G4BeamTestHadronPhysics.h"
+#include "G4BeamTestIonPhysics.h"
+#else
+#include "G4DecayPhysics.hh"
+#include "G4EmStandardPhysics.hh"
+#include "G4EmExtraPhysics.hh"
+#include "G4IonPhysics.hh"
+#include "G4StoppingPhysics.hh"
+#include "G4HadronElasticPhysics.hh"
+//#include "G4NeutronTrackingCut.hh"
+#include "G4DataQuestionaire.hh"
+#include "G4HadronPhysicsFTFP_BERT.hh"
+#include <FTFP_BERT.hh>
+#endif
+
+#include <G4ProcessManager.hh>
+#include <G4ParticleTypes.hh>
+#include <G4UserSpecialCuts.hh>
+
+
+G4BeamTestPhysicsList::G4BeamTestPhysicsList()
+  : G4VUserPhysicsList()
+{
+  defaultCutValue = 0.7*CLHEP::mm;
+  SetVerboseLevel(1);
+
+  RegisterPhysics(new G4BeamTestGeneralPhysics);
+#if G4VERSION_NUMBER < 1000
+  RegisterPhysics(new G4BeamTestEMPhysics);
+  RegisterPhysics(new G4BeamTestMuonPhysics);
+  RegisterPhysics(new G4BeamTestHadronPhysics);
+  RegisterPhysics(new G4BeamTestIonPhysics);
+#else
+  // The following is basically Geant4's FTFP_BERT physics list
+  G4DataQuestionaire it(photon); // this checks that G4LEVELGAMMADATA is there
+
+  RegisterPhysics(new G4EmStandardPhysics());
+  RegisterPhysics(new G4EmExtraPhysics());
+  RegisterPhysics(new G4DecayPhysics());
+  RegisterPhysics(new G4HadronElasticPhysics());
+  RegisterPhysics(new G4HadronPhysicsFTFP_BERT());
+  RegisterPhysics(new G4StoppingPhysics());
+  RegisterPhysics(new G4IonPhysics());
+  //RegisterPhysics(new G4NeutronTrackingCut());
+#endif
+}
+
+
+G4BeamTestPhysicsList::~G4BeamTestPhysicsList()
+{
+}
+
+
+void G4BeamTestPhysicsList::SetCuts()
+{
+  //G4VUserPhysicsList::SetCutsWithDefault method sets
+  //the default cut value for all particle types
+  //
+  SetCutsWithDefault();
+  if (verboseLevel>0) DumpCutValuesTable();
+}
+
diff --git a/src/G4BeamTestRunManager.cxx b/src/G4BeamTestRunManager.cxx
new file mode 100644
index 0000000..753deb8
--- /dev/null
+++ b/src/G4BeamTestRunManager.cxx
@@ -0,0 +1,143 @@
+/**
+ * Copyright (C) 2009
+ * The IceCube collaboration
+ * ID: $Id: G4BeamTestRunManager.cxx 162348 2018-04-25 20:09:46Z nega $
+ *
+ * @file G4BeamTestRunManager.cxx
+ * @version $Rev: 162348 $
+ * @date $Date: 2018-04-25 21:09:46 +0100 (Wed, 25 Apr 2018) $
+ * @author Tilo Waldenmaier
+ */
+
+
+// On Sun, to prevent conflict with ObjectSpace, G4Timer.hh has to be
+// loaded *before* globals.hh...
+#include "G4Timer.hh"
+
+#include <g4-tankresponse/g4classes/G4BeamTestRunManager.h>
+#include <G4ParticleGun.hh>
+#include <G4Run.hh>
+
+
+G4BeamTestRunManager::G4BeamTestRunManager(): G4RunManager()
+{
+  
+}
+
+
+void G4BeamTestRunManager::BeamOn(G4int n_event,const char* macroFile,G4int n_select)
+{
+  G4String text = "BeamOn method is not supported in G4BeamTestRunManager!";
+  G4Exception("G4BeamTestRunManager::BeamOn()", "G4BeamTestRunManager001", FatalException, text);
+}
+
+
+void G4BeamTestRunManager::InitializeRun()
+{
+  G4bool cond = ConfirmBeamOnCondition();
+  if(cond)
+  {
+    // Reset the event counter 
+    numberOfEventToBeProcessed = 0;
+    ConstructScoringWorlds();
+    RunInitialization();
+    
+    if(verboseLevel>0) timer->Start();
+  }
+}
+
+
+void G4BeamTestRunManager::InjectParticle(G4ParticleGun* particleGun)
+{
+  if(!currentRun)
+  {
+    G4String text = "Run needs to be initialized before injecting a particle.";
+    G4Exception("G4BeamTestRunManager::InjectParticle()", "G4BeamTestRunManager002", FatalException, text);
+  }
+  assert(currentRun); // the G4Exception() above calls abort(). This assert() silences the clang static analyzer
+  
+  numberOfEventToBeProcessed++;
+  currentRun->SetNumberOfEventToBeProcessed(numberOfEventToBeProcessed);
+  
+  currentEvent = GenerateEvent(numberOfEventToBeProcessed);
+  particleGun->GeneratePrimaryVertex(currentEvent);
+  
+  eventManager->ProcessOneEvent(currentEvent);
+  AnalyzeEvent(currentEvent);
+  Update_Scoring();
+  StackPreviousEvent(currentEvent);
+  currentEvent = 0;
+  
+  if(runAborted) TerminateRun();
+}
+
+
+G4Event* G4BeamTestRunManager::GenerateEvent(G4int i_event)
+{
+  G4Event* anEvent = new G4Event(i_event);
+  
+  if(storeRandomNumberStatusToG4Event==1 || storeRandomNumberStatusToG4Event==3)
+  {
+    std::ostringstream oss;
+    CLHEP::HepRandom::saveFullState(oss);
+    randomNumberStatusForThisEvent = oss.str();
+    anEvent->SetRandomNumberStatus(randomNumberStatusForThisEvent);
+  }
+  
+  if(storeRandomNumberStatus)
+  {
+    G4String fileN = randomNumberStatusDir + "currentEvent.rndm"; 
+    CLHEP::HepRandom::saveEngineStatus(fileN);
+  }
+  
+  return anEvent;
+}
+
+
+void G4BeamTestRunManager::TerminateRun()
+{
+  if(verboseLevel>0)
+  {
+    timer->Stop();
+    G4cout << "Run terminated." << G4endl;
+    G4cout << "Run Summary" << G4endl;
+    if(runAborted)
+    {
+      G4cout << "  Run Aborted after " << numberOfEventToBeProcessed << " events processed." << G4endl;
+    }
+    else
+    {
+      G4cout << "  Number of events processed : " << numberOfEventToBeProcessed << G4endl;
+    }
+    G4cout << "  "  << *timer << G4endl;
+  }
+  
+  RunTermination();
+}
+
+
+//
+// The following method is an exact copy of 
+// UpdateScoring which is private in the G4RunManager
+//
+
+#include <G4ScoringManager.hh>
+#include <G4HCofThisEvent.hh>
+#include <G4VHitsCollection.hh>
+
+void G4BeamTestRunManager::Update_Scoring()
+{
+  G4ScoringManager* ScM = G4ScoringManager::GetScoringManagerIfExist();
+  if(!ScM) return;
+  G4int nPar = ScM->GetNumberOfMesh();
+  if(nPar<1) return;
+  
+  G4HCofThisEvent* HCE = currentEvent->GetHCofThisEvent();
+  if(!HCE) return;
+  G4int nColl = HCE->GetCapacity();
+  for(G4int i=0;i<nColl;i++)
+  {
+    G4VHitsCollection* HC = HCE->GetHC(i);
+    if(HC) ScM->Accumulate(HC);
+  }
+}
diff --git a/src/G4BeamTestTank.cxx b/src/G4BeamTestTank.cxx
new file mode 100644
index 0000000..4f12c7a
--- /dev/null
+++ b/src/G4BeamTestTank.cxx
@@ -0,0 +1,382 @@
+/**
+ * Copyright (C) 2009
+ * The IceCube collaboration
+ * ID: $Id: G4BeamTestTank.cxx 158930 2017-10-20 01:17:41Z cweaver $
+ *
+ * @file G4BeamTestTank.cxx
+ * @version $Rev: 158930 $
+ * @date $Date: 2017-10-20 02:17:41 +0100 (Fri, 20 Oct 2017) $
+ * @author Tilo Waldenmaier, Thomas Melzig
+ */
+
+
+#include <g4-tankresponse/g4classes/G4BeamTestTank.h>
+#include <g4-tankresponse/g4classes/G4TankIceSD.h>
+#include <icetray/I3Units.h>
+
+#include <icetray/OMKey.h>
+#include <dataclasses/TankKey.h>
+#include <dataclasses/geometry/I3TankGeo.h>
+#include <dataclasses/geometry/I3Geometry.h>
+#include <dataclasses/I3Position.h>
+
+#include <G4LogicalVolume.hh>
+#include <G4PVPlacement.hh>
+#include <G4Material.hh>
+#include <G4Tubs.hh>
+#include <G4Sphere.hh>
+#include <G4Box.hh>
+#include <G4SDManager.hh>
+#include <G4ThreeVector.hh>
+
+#include <G4VisAttributes.hh>
+
+#include <G4UserLimits.hh>
+
+#include <boost/foreach.hpp>
+
+//prevent gcc to make something stupid with pretended unused variables
+#ifdef __GNUC__
+#define ATTRIBUTE_UNUSED __attribute__((unused))
+#else
+#define ATTRIBUTE_UNUSED
+#endif
+
+
+G4BeamTestTank::G4BeamTestTank(const TankKey& tankKey, const I3Geometry& geometry):
+tankKey_(tankKey), geometry_(geometry)
+{
+  const I3StationGeoMap& stationMap = geometry.stationgeo;
+  unsigned int tankID = tankKey.tank==TankKey::TankA?0:1;
+  I3StationGeoMap::const_iterator station_iter = stationMap.find(tankKey.string);
+  
+  if(station_iter==stationMap.end())
+  {
+    log_fatal("The requested station %d in not in the geometry!", tankKey.string);
+    return;
+  }
+  
+  if(station_iter->second.size()<tankID)
+  {
+    log_fatal("The number of tanks in station %d is not correct!", tankKey.string);
+    return;
+  }
+  const I3TankGeo& tankGeo = station_iter->second.at(tankID);
+  
+  // Get tank dimensions
+  tankThickness_ = 0.5*CLHEP::cm;
+  tankHeight_  = (tankGeo.tankheight / I3Units::m) * CLHEP::m + tankThickness_;
+  innerRadius_ = (tankGeo.tankradius / I3Units::m) * CLHEP::m;
+  outerRadius_ = innerRadius_ + tankThickness_;
+  
+  // Get fill and snow heights
+  fillHeight_    = (tankGeo.fillheight / I3Units::m) * CLHEP::m;
+  snowHeight_    = (tankGeo.snowheight / I3Units::m) * CLHEP::m;
+  perliteHeight_ = tankHeight_ - tankThickness_ - fillHeight_;
+  
+  // Set DOM dimensions
+  glassOuterRadius_ = 6.5 * 2.54 * CLHEP::cm; // 6.5" outer glass sphere radius
+  glassThickness_   = 0.5 * 2.54 * CLHEP::cm; // 0.5" glass sphere thickness
+  
+  // Calculate tank position (tank center)
+  // tankGeo.position corresponds to the average position of the two DOMs in a tank
+  position_.set((tankGeo.position.GetX() / I3Units::m) * CLHEP::m,
+                (tankGeo.position.GetY() / I3Units::m) * CLHEP::m,
+                (tankGeo.position.GetZ() / I3Units::m) * CLHEP::m - fillHeight_ + 0.5*tankHeight_);
+  
+  // Get positions of the doms relativ to tank center
+  BOOST_FOREACH(const OMKey& omKey, tankGeo.omKeyList_)
+  {
+    I3OMGeoMap::const_iterator omGeo_iter = geometry_.omgeo.find(omKey);
+    if(omGeo_iter==geometry_.omgeo.end())
+    {
+      log_error_stream(omKey << " is missing in Tank " << tankKey_);
+      continue;
+    }
+    
+    G4ThreeVector relPos(omGeo_iter->second.position.GetX() - tankGeo.position.GetX(),
+                         omGeo_iter->second.position.GetY() - tankGeo.position.GetY(),
+                         omGeo_iter->second.position.GetZ() - tankGeo.position.GetZ());
+    
+    relDomPositions_[omKey] = (relPos / I3Units::m) * CLHEP::m;
+  }
+  
+  //
+  // Calculate Delaunay points
+  //
+  G4ThreeVector triangleDir(NAN, NAN, NAN);
+  switch(station_iter->second.size())
+  {
+    case 1: // Single tank 
+    {
+      // Vector orthogonal to DOM positions
+      triangleDir.set(relDomPositions_.begin()->second.y(),
+                      -relDomPositions_.begin()->second.x(),
+                      0.0);
+      break;
+    }
+    case 2: // Two tanks
+    {
+      const I3TankGeo& neighborGeo = station_iter->second.at(tankID==0?1:0);
+      G4ThreeVector neighborPos(neighborGeo.position.GetX(),
+                                neighborGeo.position.GetY(),
+                                neighborGeo.position.GetZ());
+      
+      // Convert to G4 units
+      neighborPos *= CLHEP::m/I3Units::m;
+      
+      // Same z position same as other tank
+      neighborPos.setZ(position_.z());
+      
+      triangleDir = position_ - neighborPos;
+      break;
+    }
+    default:
+    {
+      log_fatal("Invalid number of tanks (%zu) in station %d!",
+                station_iter->second.size(), tankKey_.string);
+      break;
+    }
+  }
+  
+  // side length
+  double triangleLength = 10.0 * CLHEP::m;
+  
+  triangleDir.setMag(0.5*triangleLength/cos(CLHEP::pi/6.0));
+  delaunayPoint1_ = position_ + triangleDir;
+  
+  // Rotate by 120 deg
+  triangleDir.rotateZ(CLHEP::pi/1.5);
+  delaunayPoint2_ = position_ + triangleDir;
+  
+  // Rotate by 120 deg
+  triangleDir.rotateZ(CLHEP::pi/1.5);
+  delaunayPoint3_ = position_ + triangleDir;
+}
+
+
+G4BeamTestTank::~G4BeamTestTank()
+{
+  
+}
+
+
+G4VPhysicalVolume* G4BeamTestTank::InstallTank(G4VPhysicalVolume* mother, const G4ThreeVector& origin)
+{
+  // User limits (energy cutoffs)
+  // Do not create photons or electrons below cherenkov threshold
+  // See also corresponding UserSpecialCuts in Physicslist !!!!
+  // Maybe do all of this as stepping action ??????
+  G4UserLimits* energyLimit = new G4UserLimits();
+  energyLimit->SetUserMinEkine(280.0 * CLHEP::keV);  // Cherenkov threshold of electrons in ice
+  
+  std::string tankName=boost::lexical_cast<std::string>(tankKey_);
+  
+  // Define plastic frame
+  G4Material* plastic = G4Material::GetMaterial("Plastic");
+  G4Tubs* solidTank = new G4Tubs(("solid_tank_" + tankName).c_str(),
+                                 0.0 * CLHEP::m, outerRadius_, 0.5 * tankHeight_,
+                                 0.0 * CLHEP::deg, 360.0 * CLHEP::deg);
+  tankLog_ = new G4LogicalVolume(solidTank, plastic,
+                                 ("log_tank_" + tankName).c_str(), 0, 0, 0);
+  
+  // Define ice volume
+  G4Material* ice = G4Material::GetMaterial("Ice");
+  G4Tubs* solidIce = new G4Tubs(("solid_ice_" + tankName).c_str(),
+                                0.0 * CLHEP::m, innerRadius_, 0.5 * fillHeight_,
+                                0.0 * CLHEP::deg, 360.0 * CLHEP::deg);
+  G4LogicalVolume* logIce =
+      new G4LogicalVolume(solidIce, ice, ("log_ice_" + tankName).c_str(), 0, 0, 0);
+  G4ThreeVector physIcePosition(0, 0, -0.5*tankHeight_ + tankThickness_ + 0.5*fillHeight_);
+  G4VPhysicalVolume* physIce ATTRIBUTE_UNUSED =
+      new G4PVPlacement(0, physIcePosition, logIce,
+                        ("ice_" + tankName).c_str(), tankLog_, false, 0);
+  
+  // Define perlite volume
+  G4Material* perlite = G4Material::GetMaterial("Perlite");
+  G4Tubs* solidPerlite = new G4Tubs(("solid_perlite_" + tankName).c_str(),
+                                    0.0 * CLHEP::m, innerRadius_, 0.5 * perliteHeight_,
+                                    0.0 * CLHEP::deg, 360.0 * CLHEP::deg);
+  G4LogicalVolume* logPerlite =
+      new G4LogicalVolume(solidPerlite, perlite, ("log_perlite_" + tankName).c_str(), 0, 0, 0);
+  G4ThreeVector physPerlitePosition(0, 0, -0.5 * tankHeight_ + 0.5 * CLHEP::cm + fillHeight_ +
+                                    0.5 * perliteHeight_);
+  G4VPhysicalVolume* physPerlite ATTRIBUTE_UNUSED =
+      new G4PVPlacement(0, physPerlitePosition, logPerlite,
+                        ("perlite_" + tankName).c_str(), tankLog_, false, 0);
+  
+  // Define glass sphere & effective DOM material splitted in upper and lower part
+  G4Material* glass        = G4Material::GetMaterial("Glass");
+  G4Material* effectiveDOM = G4Material::GetMaterial("effectiveDOM");
+  
+  std::map<OMKey, G4ThreeVector> domPosIce;
+  std::map<OMKey, G4ThreeVector>::const_iterator om_iter;
+  for(om_iter=relDomPositions_.begin(); om_iter!=relDomPositions_.end(); ++om_iter)
+  {
+    const OMKey& omKey = om_iter->first;
+    
+    G4ThreeVector upperDOMpos(om_iter->second.x(), om_iter->second.y(), -0.5 * perliteHeight_);
+    G4ThreeVector lowerDOMpos(om_iter->second.x(), om_iter->second.y(),  0.5 * fillHeight_);
+    
+    domPosIce[omKey] = lowerDOMpos;
+    
+    std::string omName=boost::lexical_cast<std::string>(omKey);
+    
+    G4Sphere *upperglasssphere = new G4Sphere (("solid_dom_up_" + omName).c_str(),
+                                               0.0 * CLHEP::m, glassOuterRadius_,
+                                               0.0 * CLHEP::deg, 360.0 * CLHEP::deg,
+                                               0.0 * CLHEP::deg,  90.0 * CLHEP::deg);
+    G4Sphere *lowerglasssphere = new G4Sphere (("solid_dom_lo_" + omName).c_str(),
+                                               0.0  * CLHEP::m, glassOuterRadius_,
+                                               0.0  * CLHEP::deg, 360.0 * CLHEP::deg,
+                                               90.0 * CLHEP::deg, 180.0 * CLHEP::deg);
+    
+    G4double domInnerRadius = glassOuterRadius_ - glassThickness_;
+    G4Sphere *upperdomsphere = new G4Sphere (("solid_inside_dom_up_" + omName).c_str(),
+                                             0.0 * CLHEP::m, domInnerRadius,
+                                             0.0 * CLHEP::deg, 360.0 * CLHEP::deg,
+                                             0.0 * CLHEP::deg,  90.0 * CLHEP::deg);
+    G4Sphere *lowerdomsphere = new G4Sphere (("solid_inside_dom_lo_" + omName).c_str(),
+                                             0.0  * CLHEP::m, domInnerRadius,
+                                             0.0  * CLHEP::deg, 360.0 * CLHEP::deg,
+                                             90.0 * CLHEP::deg, 180.0 * CLHEP::deg);
+    
+    G4LogicalVolume* logUpperGlass =
+        new G4LogicalVolume(upperglasssphere, glass,
+                            ("log_dom_up_" + omName).c_str(), 0, 0, 0);
+    G4LogicalVolume* logLowerGlass =
+        new G4LogicalVolume(lowerglasssphere, glass,
+                            ("log_dom_lo_" + omName).c_str(), 0, 0, 0);
+    G4LogicalVolume* logUpperDOM =
+        new G4LogicalVolume(upperdomsphere, effectiveDOM,
+                            ("log_inside_dom_up_" + omName).c_str(), 0, 0, 0);
+    G4LogicalVolume* logLowerDOM =
+        new G4LogicalVolume(lowerdomsphere, effectiveDOM,
+                            ("log_inside_dom_lo_" + omName).c_str(), 0, 0, 0);
+    G4VPhysicalVolume* physUpperGlass ATTRIBUTE_UNUSED =
+        new G4PVPlacement(0, upperDOMpos, logUpperGlass,
+                          ("dom_up_" + omName).c_str(), logPerlite, false, 0);
+    G4VPhysicalVolume* physLowerGlass ATTRIBUTE_UNUSED =
+        new G4PVPlacement(0, lowerDOMpos, logLowerGlass,
+                          ("dom_lo_" + omName).c_str(), logIce, false, 0);
+    G4VPhysicalVolume* physUpperDOM ATTRIBUTE_UNUSED =
+        new G4PVPlacement(0, G4ThreeVector(0,0,0), logUpperDOM,
+                          ("inside_dom_up_" + omName).c_str(), logUpperGlass, false, 0);
+    G4VPhysicalVolume* physLowerDOM ATTRIBUTE_UNUSED =
+        new G4PVPlacement(0, G4ThreeVector(0,0,0), logLowerDOM,
+                          ("inside_dom_lo_" + omName).c_str(), logLowerGlass, false, 0);
+    
+    // apply energy limits
+    logUpperGlass->SetUserLimits(energyLimit);
+    logLowerGlass->SetUserLimits(energyLimit);
+    logUpperDOM->SetUserLimits(energyLimit);
+    logLowerDOM->SetUserLimits(energyLimit);
+  }
+  
+  // Define sensitive detector
+  G4SDManager* sdManager = G4SDManager::GetSDMpointer();
+  iceSD_ = new G4TankIceSD(("ice_SD_" + tankName).c_str(), domPosIce);
+  sdManager->AddNewDetector(iceSD_);
+  logIce->SetSensitiveDetector(iceSD_);
+  
+  // Instantiation of a set of visualization attributes with red colour
+  G4VisAttributes * tankVisAtt = new G4VisAttributes(G4Colour(1,0,0));
+  // Set the forced wireframe style 
+  //snowVisAtt->SetForceWireFrame(true);
+  // Assignment of the visualization attributes to the logical volume
+  tankLog_->SetVisAttributes(tankVisAtt);
+  
+  G4ThreeVector tankPos = position_ - origin - mother->GetTranslation();
+  
+  G4VPhysicalVolume* tankPhys = new G4PVPlacement(0, tankPos, tankLog_,
+                                                  ("tank_" + tankName).c_str(),
+                                                  mother->GetLogicalVolume(), false, 0);
+  
+  // apply energy limits
+  tankLog_->SetUserLimits(energyLimit);
+  logPerlite->SetUserLimits(energyLimit);
+  logIce->SetUserLimits(energyLimit);
+  
+  return tankPhys;
+}
+
+
+double G4BeamTestTank::GetNumCherenkov(const OMKey& omKey)
+{
+  return std::max(iceSD_->GetNumCherenkov(omKey), 0.);
+}
+
+
+double G4BeamTestTank::GetNumCherenkovWeight(const OMKey& omKey)
+{
+  return std::max(iceSD_->GetNumCherenkovWeight(omKey), 0.);
+}
+
+
+double G4BeamTestTank::GetEDep_G4(const OMKey& omKey)
+{
+  return std::max(iceSD_->GetEDep(omKey), 0.);
+}
+
+
+double G4BeamTestTank::GetTime_G4(const OMKey& omKey)
+{
+  return iceSD_->GetTime(omKey);
+}
+
+
+double G4BeamTestTank::GetEDep_I3(const OMKey& omKey)
+{
+  return std::max(iceSD_->GetEDep(omKey), 0.) / CLHEP::keV * I3Units::keV;
+}
+
+
+double G4BeamTestTank::GetTime_I3(const OMKey& omKey)
+{
+  return (iceSD_->GetTime(omKey) / CLHEP::s) * I3Units::s;
+}
+
+
+I3Position G4BeamTestTank::GetPos_I3()
+{
+  I3Position pos((position_.x() / CLHEP::m) * I3Units::m,
+                 (position_.y() / CLHEP::m) * I3Units::m,
+                 (position_.z() / CLHEP::m) * I3Units::m);
+  return pos;
+}
+
+
+double G4BeamTestTank::GetX_I3()
+{
+  return (position_.x() / CLHEP::m) * I3Units::m;
+}
+
+
+double G4BeamTestTank::GetY_I3()
+{
+  return (position_.y() / CLHEP::m) * I3Units::m;
+}
+
+
+double G4BeamTestTank::GetZ_I3()
+{
+  return (position_.z() / CLHEP::m) * I3Units::m;
+}
+
+
+double G4BeamTestTank::GetTankHeight_I3()
+{
+  return (tankHeight_ / CLHEP::m) * I3Units::m;
+}
+
+
+double G4BeamTestTank::GetTankRadius_I3()
+{
+  return (outerRadius_ / CLHEP::m) * I3Units::m;
+}
+
+
+double G4BeamTestTank::GetSnowHeight_I3()
+{
+  return (snowHeight_ / CLHEP::m) * I3Units::m;
+}
diff --git a/src/G4BeamTestUserSteppingAction.cxx b/src/G4BeamTestUserSteppingAction.cxx
new file mode 100644
index 0000000..b16098c
--- /dev/null
+++ b/src/G4BeamTestUserSteppingAction.cxx
@@ -0,0 +1,42 @@
+/**
+ * Copyright (C) 2011
+ * The IceCube collaboration
+ * ID: $Id$
+ *
+ * @file G4BeamTestUserSteppingAction.cxx
+ * @version $Revision$
+ * @date $Date$
+ * @author Thomas Melzig
+ *
+ * $LastChangedBy$
+ */
+
+
+#include <g4-tankresponse/g4classes/G4BeamTestUserSteppingAction.h>
+
+#include "G4Step.hh"
+#include "G4Track.hh"
+#include "G4LogicalVolume.hh"
+#include "G4UserLimits.hh"
+#include "G4SteppingManager.hh"
+
+G4BeamTestUserSteppingAction::G4BeamTestUserSteppingAction(){}
+
+void G4BeamTestUserSteppingAction::UserSteppingAction(const G4Step*)
+{
+  G4Track* track = fpSteppingManager->GetTrack();
+  if(track)
+  {
+    const G4LogicalVolume *volume = track->GetLogicalVolumeAtVertex();
+    G4UserLimits *limit = volume->GetUserLimits();
+    if(!limit) G4cout << "----> G4LogicalVolume: " << volume->GetName() << " has no defined G4UserLimit" << G4endl;
+    G4double threshold = limit->GetUserMinEkine(*track);
+    //check if particle is a gamma
+    if(track->GetDefinition()->GetParticleName() == "gamma")
+    {
+      //check if particle energy is below threshold; if true, kill the particle
+      G4double energy = track->GetTotalEnergy();
+      if(energy < threshold) track->SetTrackStatus(fStopAndKill);
+    }
+  }
+}
diff --git a/src/G4BeamTestUserTrackingAction.cxx b/src/G4BeamTestUserTrackingAction.cxx
new file mode 100644
index 0000000..d554152
--- /dev/null
+++ b/src/G4BeamTestUserTrackingAction.cxx
@@ -0,0 +1,50 @@
+/**
+ * Copyright (C) 2011
+ * The IceCube collaboration
+ * ID: $Id$
+ *
+ * @file G4BeamTestUserTrackingAction.cxx
+ * @version $Revision$
+ * @date $Date$
+ * @author Thomas Melzig
+ *
+ * $LastChangedBy$
+ */
+
+
+#include <g4-tankresponse/g4classes/G4BeamTestUserTrackingAction.h>
+
+#include "G4Track.hh"
+#include "G4UserLimits.hh"
+#include "G4TrackVector.hh"
+#include "G4TrackingManager.hh"
+
+G4BeamTestUserTrackingAction::G4BeamTestUserTrackingAction(){}
+
+void G4BeamTestUserTrackingAction::PreUserTrackingAction(const G4Track*){}
+
+void G4BeamTestUserTrackingAction::PostUserTrackingAction(const G4Track* track)
+{
+  const G4LogicalVolume *volume = track->GetLogicalVolumeAtVertex();
+  G4UserLimits *limit = volume->GetUserLimits();
+  if(!limit) G4cout << "----> G4LogicalVolume: " << volume->GetName() << " has no defined G4UserLimit" << G4endl;
+  G4double threshold = limit->GetUserMinEkine(*track);
+  G4TrackVector* secondaries = fpTrackingManager->GimmeSecondaries();
+  if(secondaries)
+  {
+    size_t nSeco = secondaries->size();
+    if(nSeco>0)
+    {
+      for(size_t i=0;i<nSeco;i++)
+      {
+        //check if secondary particle is a gamma
+        if((*secondaries)[i]->GetDefinition()->GetParticleName() == "gamma")
+        {
+          //check if particle energy is below threshold; if true, kill the particle
+          G4double energy = (*secondaries)[i]->GetTotalEnergy();
+          if(energy < threshold) (*secondaries)[i]->SetTrackStatus(fStopAndKill);
+        }
+      }
+    }
+  }
+}
diff --git a/src/G4Interface.cxx b/src/G4Interface.cxx
new file mode 100644
index 0000000..36ba3d9
--- /dev/null
+++ b/src/G4Interface.cxx
@@ -0,0 +1,292 @@
+/**
+ * Copyright (C) 2009
+ * The IceCube collaboration
+ * ID: $Id: G4Interface.cxx 149388 2016-08-18 21:50:04Z jgonzalez $
+ *
+ * @file G4Interface.cxx
+ * @version $Rev: 149388 $
+ * @date $Date: 2016-08-18 22:50:04 +0100 (Thu, 18 Aug 2016) $
+ * @author Tilo Waldenmaier
+ */
+
+#include <dataclasses/physics/I3Particle.h>
+
+#include <g4-tankresponse/g4classes/G4Interface.h>
+#include <g4-tankresponse/g4classes/G4BeamTestDetectorConstruction.h>
+#include <g4-tankresponse/g4classes/G4BeamTestPhysicsList.h>
+#include <g4-tankresponse/g4classes/G4BeamTestUserTrackingAction.h>
+#include <g4-tankresponse/g4classes/G4BeamTestUserSteppingAction.h>
+
+#include <icetray/I3Logging.h>
+
+#ifdef G4VIS_USE
+#include <G4VisExecutive.hh>
+#endif
+
+#include <G4ParticleGun.hh>
+#include <G4ParticleTable.hh>
+#include <G4ParticleDefinition.hh>
+#include <G4UImanager.hh>
+
+
+G4Interface* G4Interface::g4Interface_ = NULL;
+
+G4Interface::G4Interface(const std::string& visMacro):
+  detector_(NULL), initialized_(false),
+  eventInitialized_(false), visMacro_(visMacro)
+{
+  g4Interface_ = this;
+
+  // Visualization manager
+#ifdef G4VIS_USE
+  visManager_ = NULL;
+  if(!visMacro_.empty())
+  {
+    visManager_ = new G4VisExecutive();
+    visManager_->Initialize();
+  }
+#endif
+}
+
+
+G4Interface::~G4Interface()
+{
+  g4Interface_ = NULL;
+
+#ifdef G4VIS_USE
+  if(visManager_) delete visManager_;
+#endif
+}
+
+
+void G4Interface::InstallTank(G4BeamTestTank* tank)
+{
+  if(initialized_)
+  {
+    log_fatal("G4Interface aleady initialized. Cannot install tank!");
+    return;
+  }
+
+  if(!detector_) detector_ = new G4BeamTestDetectorConstruction();
+  detector_->InstallTank(tank);
+}
+
+
+void G4Interface::InitializeEvent()
+{
+  ///
+  /// An IceTray EVENT corresponds to a G4RUN
+  /// where each injected particle initiates an G4EVENT
+  ///
+
+  if(!initialized_)
+  {
+    Initialize();
+  }
+
+  if(!eventInitialized_)
+  {
+    runManager_.InitializeRun();
+    eventInitialized_ = true;
+  }
+}
+
+
+void G4Interface::InjectParticle(const I3Particle& particle)
+{
+  if(!eventInitialized_)
+  {
+    log_fatal("No event initialized. Cannot inject particle!");
+    return;
+  }
+
+  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
+  G4ParticleDefinition* particleDef = NULL;
+  switch(particle.GetType())
+  {
+    case I3Particle::Gamma:
+       particleDef = particleTable->FindParticle("gamma");
+       break;
+    case I3Particle::EMinus:
+       particleDef = particleTable->FindParticle("e-");
+       break;
+    case I3Particle::EPlus:
+       particleDef = particleTable->FindParticle("e+");
+       break;
+    case I3Particle::MuMinus:
+       particleDef = particleTable->FindParticle("mu-");
+       break;
+    case I3Particle::MuPlus:
+       particleDef = particleTable->FindParticle("mu+");
+       break;
+    case I3Particle::PPlus:
+       particleDef = particleTable->FindParticle("proton");
+       break;
+    case I3Particle::PMinus:
+       particleDef = particleTable->FindParticle("anti_proton");
+       break;
+    case I3Particle::Neutron:
+       particleDef = particleTable->FindParticle("neutron");
+       break;
+#ifdef I3PARTICLE_SUPPORTS_PDG_ENCODINGS
+    case I3Particle::NeutronBar:
+#else
+    case 25:
+#endif
+       particleDef = particleTable->FindParticle("anti_neutron");
+       break;
+    case I3Particle::PiPlus:
+       particleDef = particleTable->FindParticle("pi+");
+       break;
+    case I3Particle::PiMinus:
+       particleDef = particleTable->FindParticle("pi-");
+       break;
+    case I3Particle::Pi0:
+       particleDef = particleTable->FindParticle("pi0");
+       break;
+    case I3Particle::KPlus:
+       particleDef = particleTable->FindParticle("kaon+");
+       break;
+    case I3Particle::KMinus:
+       particleDef = particleTable->FindParticle("kaon-");
+       break;
+    case I3Particle::K0_Long:
+       particleDef = particleTable->FindParticle("kaon0L");
+       break;
+    case I3Particle::K0_Short:
+       particleDef = particleTable->FindParticle("kaon0S");
+       break;
+    case I3Particle::NuE:
+       particleDef = particleTable->FindParticle("nu_e");
+       break;
+    case I3Particle::NuEBar:
+       particleDef = particleTable->FindParticle("anti_nu_e");
+       break;
+    case I3Particle::NuMu:
+       particleDef = particleTable->FindParticle("nu_mu");
+       break;
+    case I3Particle::NuMuBar:
+       particleDef = particleTable->FindParticle("anti_nu_mu");
+       break;
+    case I3Particle::NuTau:
+       particleDef = particleTable->FindParticle("nu_tau");
+       break;
+    case I3Particle::NuTauBar:
+       particleDef = particleTable->FindParticle("anti_nu_tau");
+       break;
+    default:
+      log_warn("Man, check out that strange particle \"%s\" ?!", particle.GetTypeString().c_str());
+      return;
+  }
+  
+  // Particle position in G4 units
+  G4ThreeVector position((particle.GetX() / I3Units::m) * CLHEP::m,
+                         (particle.GetY() / I3Units::m) * CLHEP::m,
+                         (particle.GetZ() / I3Units::m) * CLHEP::m);
+
+  // Transform I3 coorinates to world system
+  position -= detector_->GetWorldOrigin();
+
+  G4ThreeVector direction(particle.GetDir().GetX(),
+                          particle.GetDir().GetY(),
+                          particle.GetDir().GetZ());
+
+  G4ParticleGun gun(1);
+  gun.SetParticleDefinition(particleDef);
+  gun.SetParticleEnergy((particle.GetEnergy() / I3Units::GeV) * CLHEP::GeV);
+  gun.SetParticlePosition(position);
+  gun.SetParticleMomentumDirection(direction);
+
+  log_trace("Injecting %s: x=%.2f m, y=%.2f m, z=%.2f m, E=%.3f MeV",
+            particle.GetTypeString().c_str(),
+            position.x() / CLHEP::m,
+            position.y() / CLHEP::m,
+            position.z() / CLHEP::m,
+            gun.GetParticleEnergy() / CLHEP::MeV);
+
+  runManager_.InjectParticle(&gun);
+}
+
+
+void G4Interface::TerminateEvent()
+{
+  ///
+  /// An IceTray EVENT corresponds to a G4RUN
+  /// where each injected particle initiates an G4EVENT
+  ///
+
+  if(eventInitialized_)
+  {
+    runManager_.TerminateRun();
+    eventInitialized_ = false;
+  }
+}
+
+
+void G4Interface::Initialize()
+{
+  if(initialized_)
+  {
+    log_error("G4Interface has already been initialized. Ignoring this call!");
+    return;
+  }
+
+  log_debug("Init geometry ...");
+  runManager_.SetUserInitialization(detector_);
+
+  log_debug("Init physics list ...");
+  runManager_.SetUserInitialization(new G4BeamTestPhysicsList());
+
+  log_debug("Init UserTrackingAction ...");
+  runManager_.SetUserAction(new G4BeamTestUserTrackingAction());
+
+  log_debug("Init UserSteppingAction ...");
+  runManager_.SetUserAction(new G4BeamTestUserSteppingAction());
+
+  // Initialize G4 kernel
+  log_debug("Init run manager ...");
+  runManager_.Initialize();
+
+  // Set verbosity
+  int verboseLevel = 0;
+  switch (GetIcetrayLogger()->LogLevelForUnit("G4Interface"))
+  {
+    case I3LOG_FATAL:
+    case I3LOG_ERROR:
+    case I3LOG_WARN:
+    case I3LOG_INFO:
+    case I3LOG_NOTICE:
+    default:
+      verboseLevel = 0;
+      break;
+    case I3LOG_DEBUG:
+      verboseLevel = 1;
+      break;
+    case I3LOG_TRACE:
+      verboseLevel = 2;
+      break;
+  }
+
+  runManager_.SetVerboseLevel(verboseLevel);
+  G4EventManager::GetEventManager()->SetVerboseLevel(verboseLevel);
+  G4EventManager::GetEventManager()->GetStackManager()->SetVerboseLevel(verboseLevel);
+  G4EventManager::GetEventManager()->GetTrackingManager()->SetVerboseLevel(verboseLevel);
+#ifdef G4VIS_USE
+  if(visManager_) visManager_->SetVerboseLevel(verboseLevel);
+#endif
+
+  // Execute visualization macro (if specified)
+  if(!visMacro_.empty())
+  {
+    G4UImanager* uim = G4UImanager::GetUIpointer();
+
+    // Checking geometry
+    uim->ApplyCommand("/geometry/test/grid_test");
+
+    // Execute visualization macro
+    std::string visCmd = "/control/execute " + visMacro_;
+    uim->ApplyCommand(visCmd.c_str());
+  }
+
+  initialized_ = true;
+}
diff --git a/src/G4TankIceSD.cxx b/src/G4TankIceSD.cxx
new file mode 100644
index 0000000..4a1be0b
--- /dev/null
+++ b/src/G4TankIceSD.cxx
@@ -0,0 +1,176 @@
+/**
+ * Copyright (C) 2009
+ * The IceCube collaboration
+ * ID: $Id: G4TankIceSD.cxx 152814 2017-01-19 21:34:52Z jgonzalez $
+ *
+ * @file G4TankIceSD.cxx
+ * @version $Rev: 152814 $
+ * @date $Date: 2017-01-19 21:34:52 +0000 (Thu, 19 Jan 2017) $
+ * @author Tilo Waldenmaier, Thomas Melzig
+ */
+
+
+#include <g4-tankresponse/g4classes/G4TankIceSD.h>
+#include <g4-tankresponse/g4classes/G4BeamTestTank.h>
+
+#include <G4Step.hh>
+#include <G4HCofThisEvent.hh>
+#include <G4TouchableHistory.hh>
+#include <G4ios.hh>
+#include <iostream>
+#include <ios>
+#include <fstream>
+
+#include <G4VProcess.hh>
+#include <G4OpticalPhoton.hh>
+#include "G4Poisson.hh"
+
+#include <icetray/I3Units.h>
+
+#include "G4Version.hh"
+
+#if G4VERSION_NUMBER >= 950
+// The G4MaterialPropertyVector is gone since 4.9.5.
+// It has been typedef'd to a G4UnorderedPhysicsVector
+// with a different interface. Try to support both
+// versions with an ifdef.
+#define MATERIAL_PROPERTY_VECTOR_IS_PHYSICS_VECTOR
+#endif
+
+
+G4TankIceSD::G4TankIceSD(G4String name, const std::map<OMKey, G4ThreeVector>& domPositions)
+  : G4VSensitiveDetector(name), domPositions_(domPositions)
+{}
+
+
+G4TankIceSD::~G4TankIceSD() {}
+
+
+void G4TankIceSD::Initialize(G4HCofThisEvent* HCE)
+{
+  std::map<OMKey, G4ThreeVector>::const_iterator om_iter;
+  for(om_iter=domPositions_.begin(); om_iter!=domPositions_.end(); ++om_iter)
+  {
+    sumEdep_[om_iter->first] = 0;
+    cogTime_[om_iter->first] = 0;
+    cherenkovCounter_[om_iter->first] = 0;
+    cherenkovCounterWeight_[om_iter->first] = 0;
+  }
+}
+
+
+G4bool G4TankIceSD::ProcessHits(G4Step* aStep, G4TouchableHistory*)
+{
+  // Get energy deposit/time --> CHECK WHAT THESE QUANTITIES ACTUALLY MEAN !!!!!!!!!!!!!!!!!!!!!!!!
+  G4double edep = aStep->GetTotalEnergyDeposit();
+  G4double time = aStep->GetPreStepPoint()->GetGlobalTime();
+  G4double cherenkovNumber = GetCerenkovNumber(aStep);
+
+  if(edep<=0 && cherenkovNumber<=0) return false;
+
+  // Get Position relative to ice center
+  G4ThreeVector preStepPoint = aStep->GetPreStepPoint()->GetPosition();
+  G4ThreeVector postStepPoint = aStep->GetPostStepPoint()->GetPosition();
+  G4TouchableHandle theTouchable = aStep->GetPreStepPoint()->GetTouchableHandle();
+  G4ThreeVector worldPosition = (preStepPoint + postStepPoint) / 2.0;
+  G4ThreeVector localPosition = theTouchable->GetHistory()->
+                                GetTopTransform().TransformPoint(worldPosition);
+
+  std::map<OMKey, G4ThreeVector>::const_iterator om_iter;
+  for(om_iter=domPositions_.begin(); om_iter!=domPositions_.end(); ++om_iter)
+  {
+    G4double radius = sqrt(pow(om_iter->second.x() - localPosition.x(), 2) +
+                           pow(om_iter->second.y() - localPosition.y(), 2));
+    G4double height = (om_iter->second.z() - localPosition.z());
+
+    sumEdep_[om_iter->first] += edep;
+    cogTime_[om_iter->first] += edep*time;
+    cherenkovCounterWeight_[om_iter->first] += GetProbability(radius, height) * cherenkovNumber;
+    cherenkovCounter_[om_iter->first] += cherenkovNumber;
+  }
+  return true;
+}
+
+
+void G4TankIceSD::EndOfEvent(G4HCofThisEvent*)
+{
+  std::map<OMKey, G4ThreeVector>::const_iterator om_iter;
+  for(om_iter=domPositions_.begin(); om_iter!=domPositions_.end(); ++om_iter)
+  {
+    if(sumEdep_[om_iter->first]>0)
+    {
+      cogTime_[om_iter->first] /= sumEdep_[om_iter->first];
+    }
+  }
+}
+
+
+G4double G4TankIceSD::GetCerenkovNumber(G4Step* aStep)
+{
+  // same function as in "source/processes/electromagnetic/xrays/src/G4Cerenkov.cc" but without
+  // cerenkov angle and only for an energy independent refraction index
+  G4Track* aTrack = aStep->GetTrack();
+  const G4DynamicParticle* aParticle = aTrack->GetDynamicParticle();
+  const G4Material* aMaterial = aTrack->GetMaterial();
+
+  G4StepPoint* pPreStepPoint  = aStep->GetPreStepPoint();
+  G4StepPoint* pPostStepPoint = aStep->GetPostStepPoint();
+
+  G4MaterialPropertiesTable* aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable();
+  if (!aMaterialPropertiesTable) return 0.0;
+
+  G4MaterialPropertyVector* Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
+  if (!Rindex) return 0.0;
+
+  const G4double Rfact = 369.81 / (CLHEP::eV * CLHEP::cm);
+  const G4double charge = aParticle->GetDefinition()->GetPDGCharge();
+  const G4double beta = (pPreStepPoint->GetBeta() + pPostStepPoint->GetBeta()) / 2.;
+  if (beta <= 0.0) return 0.0;
+  G4double BetaInverse = 1. / beta;
+
+  // Min and Max photon energies
+#ifdef MATERIAL_PROPERTY_VECTOR_IS_PHYSICS_VECTOR
+  G4double Pmin = Rindex->GetMinLowEdgeEnergy();
+  G4double Pmax = Rindex->GetMaxLowEdgeEnergy();
+#else
+  G4double Pmin = Rindex->GetMinPhotonEnergy();
+  G4double Pmax = Rindex->GetMaxPhotonEnergy();
+#endif
+
+  // Min and Max Refraction Indices
+#ifdef MATERIAL_PROPERTY_VECTOR_IS_PHYSICS_VECTOR
+  G4double nMin = Rindex->GetMinValue();
+  G4double nMax = Rindex->GetMaxValue();
+#else
+  G4double nMin = Rindex->GetMinProperty();
+  G4double nMax = Rindex->GetMaxProperty();
+#endif
+  if (nMin!=nMax)
+  {
+    log_error("Support for energy dependent refraction index not yet implemented!");
+    return 0.0;
+  }
+
+  // If n(Pmax) < 1/Beta -- no photons generated
+  if (nMax < BetaInverse) return 0.0;
+
+  G4double MeanNumberOfPhotons = Rfact * charge/CLHEP::eplus * charge/CLHEP::eplus * (Pmax - Pmin) *
+                                 (1. - BetaInverse * BetaInverse / nMin / nMin);
+
+  G4double step_length = aStep->GetStepLength();
+
+  return MeanNumberOfPhotons * step_length;
+}
+
+
+G4double G4TankIceSD::GetProbability(G4double radius, G4double height)
+{
+  height = 0.90 - height / CLHEP::m;
+  radius = radius / CLHEP::m;
+  G4double p0 = 0.0340 +  0.0120 * height;
+  G4double p1 = 0.0400 + 0.000622 * exp(-1.13 + 8.03 * height);
+  G4double p2 = 0.917 + 2.34 * exp(-1.82 + 3.12 * height);
+  G4double p3 = -0.00325 - 0.00199 * height - 0.0121 * height*height;
+
+  return (p0 + p3 * radius) + (p1 - p0) * exp(-p2*p2 * radius*radius);
+}