initial commit

13 files changed, 2123 insertions, 0 deletions

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);
+}