Thu 22 Aug 01:37:19 BST 2019

32 files changed, 456 insertions, 2315 deletions

diff --git a/#G4BeamTest.cxx# b/#G4BeamTest.cxx#
deleted file mode 100644
index 287d689..0000000
--- a/#G4BeamTest.cxx#
+++ /dev/null
@@ -1,100 +0,0 @@
-#include "G4UImanager.hh"
-#include "G4UIExecutive.hh"
-
-#include "G4Interface.h"
-#include "G4BeamTestTank.h"
-
-//....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);
-
-  // Seed the random number generator manually
-  G4Random::setTheSeed(myseed);
-
-  // Initialize G4 kernel
-  //
-  G4Interface *g4Interface_ = G4Interface::GetInstance();
-  if (!g4Interface_) g4Interface_ = new G4Interface(macro);
-  G4BeamTestTank *g4Tank_ = new G4BeamTestTank(/* tankKey_, geometry */);
-  g4Interface_->InstallTank(g4Tank_);
-  g4Interface_->InitializeEvent();
-
-  // 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 G4BeamTest.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 !
-  delete g4Interface_;
-
-  return 0;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/#run1.mac# b/#run1.mac#
deleted file mode 100644
index db94b6e..0000000
--- a/#run1.mac#
+++ /dev/null
@@ -1,14 +0,0 @@
-/vis/ogl/set/displayListLimit 1
-
-/gps/particle e-
-/gps/energy 8 GeV
-
-/gps/pos/type Plane
-/gps/pos/shape Circle
-/gps/pos/centre 60 60 -25 cm
-/gps/pos/radius 0.5 cm
-/gps/ang/type iso
-
-/gps/direction 1 0 0
-/gps/position -2 0 0 m
-/run/beamOn 1
diff --git a/G4BeamTest.cxx b/G4BeamTest.cxx
index 287d689..6426b65 100644
--- a/G4BeamTest.cxx
+++ b/G4BeamTest.cxx
@@ -3,12 +3,15 @@
 
 #include "G4Interface.h"
 #include "G4BeamTestTank.h"
+#include "G4BeamTestSiHit.h"
+
+std::fstream testnew;
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 namespace {
   void PrintUsage() {
     G4cerr << " Usage: " << G4endl;
-    G4cerr << " OpNovice [-m macro ] [-u UIsession] [-t nThreads] [-r seed] "
+    G4cerr << " OpNovice [-m macro ] [-u UIsession] [-t nThreads] [-r seed] [-n outName]"
            << G4endl;
     G4cerr << "   note: -t option is available only for multi-threaded mode."
            << G4endl;
@@ -33,10 +36,14 @@ int main(int argc,char** argv)
 #endif
 
   G4long myseed = 345354;
+  std::string outName = "./testnew.txt";
   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]);
+     else if ( G4String(argv[i]) == "-n" ) {
+         outName = "./"+std::string(argv[i+1]);
+     }
 #ifdef G4MULTITHREADED
      else if ( G4String(argv[i]) == "-t" ) {
                     nThreads = G4UIcommand::ConvertToInt(argv[i+1]);
@@ -48,6 +55,8 @@ int main(int argc,char** argv)
     }
   }
 
+  testnew.open(outName, std::ofstream::out);
+
   // Choose the Random engine
   //
   G4Random::setTheEngine(new CLHEP::RanecuEngine);
diff --git a/I3G4TankResponse.cxx.backup b/I3G4TankResponse.cxx.backup
deleted file mode 100644
index 53cee3d..0000000
--- a/I3G4TankResponse.cxx.backup
+++ /dev/null
@@ -1,360 +0,0 @@
-#include <topsimulator/interface/I3IceTopResponseFactory.h>
-#include <topsimulator/GeoFunctions.h>
-#include "G4IceTopTank.h"
-#include "G4Interface.h"
-#include "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/#G4BeamTestEventAction.h# b/include/#G4BeamTestEventAction.h#
deleted file mode 100644
index 2c8b933..0000000
--- a/include/#G4BeamTestEventAction.h#
+++ /dev/null
@@ -1,44 +0,0 @@
-#ifndef G4BeamTestEventAction_h
-#define G4BeamTestEventAction_h 1
-
-#include "G4UserEventAction.hh"
-#include "globals.hh"
-
-/// Event action class
-///
-
-class G4BeamTestEventAction : public G4UserEventAction
-{
-  public:
-    G4BeamTestEventAction();
-    virtual ~G4BeamTestEventAction();
-    
-    virtual void BeginOfEventAction(const G4Event* );
-    virtual void EndOfEventAction(const G4Event* );
-
-    void AddEdep(G4double edep) { fEdep += edep; }
-    void AddTime(G4double time) { ftime += ->GetTime()}
-    void AddPath(G4double path) { fIntegralZ +=path; }
-    G4double GetPath(){return fIntegralZ;}
-    void SetXY (G4double xhit, G4double yhit) {fXIn=xhit;fYIn=yhit;}
-	G4double GetX()const {return fXIn;}
-
-
-	G4double GetY()const {return fYIn;}
-  private:
-  G4double  fEdep;
-  G4double  fIntegralZ;
-  G4double  fXIn;
-  G4double ftime;
-  G4int SiCollID;
-  G4int hcID;
-
-	G4double  fYIn;
-};
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-#endif
-
-    
-
diff --git a/include/G4BeamTestDetectorConstruction.h b/include/G4BeamTestDetectorConstruction.h
index 7e109b2..c2720ce 100644
--- a/include/G4BeamTestDetectorConstruction.h
+++ b/include/G4BeamTestDetectorConstruction.h
@@ -5,6 +5,7 @@
 #include <G4ThreeVector.hh>
 
 #include "G4BeamTestTank.h"
+#include "G4BeamTestSC4SD.h"
 
 class G4BeamTestDetectorConstruction: public G4VUserDetectorConstruction
 {
@@ -26,12 +27,8 @@ class G4BeamTestDetectorConstruction: public G4VUserDetectorConstruction
   void CreateMaterials();
 
   void CreateAir();
-  /* void CreateIce(); */
-  /* void CreateSnow(); */
   void CreateWater();
   void CreatePlastic();
-  /* void CreateTyvek(); */
-  /* void CreatePerlite(); */
   void CreateGlassSphere();
   void CreateEffectiveDOMMaterial();
   void CreateSC4();
@@ -41,6 +38,7 @@ class G4BeamTestDetectorConstruction: public G4VUserDetectorConstruction
   G4int verboseLevel_;
 
   G4BeamTestTank* tank_;
+  G4BeamTestSC4SD* sc4SD_;
 };
 
 #endif
diff --git a/include/G4BeamTestEMPhysics.h b/include/G4BeamTestEMPhysics.h
index bf20b62..24c5006 100644
--- a/include/G4BeamTestEMPhysics.h
+++ b/include/G4BeamTestEMPhysics.h
@@ -41,11 +41,11 @@ private:
   G4eMultipleScattering electronMultipleScattering;
   G4eIonisation electronIonisation;
   G4eBremsstrahlung electronBremsStrahlung;
-  
+
   //Positron physics
   G4eMultipleScattering positronMultipleScattering;
-  G4eIonisation positronIonisation; 
-  G4eBremsstrahlung positronBremsStrahlung;  
+  G4eIonisation positronIonisation;
+  G4eBremsstrahlung positronBremsStrahlung;
   G4eplusAnnihilation annihilation;
 };
 
diff --git a/include/G4BeamTestEventAction.h b/include/G4BeamTestEventAction.h
index 2d6de49..69f0af7 100644
--- a/include/G4BeamTestEventAction.h
+++ b/include/G4BeamTestEventAction.h
@@ -12,7 +12,7 @@ class G4BeamTestEventAction : public G4UserEventAction
   public:
     G4BeamTestEventAction();
     virtual ~G4BeamTestEventAction();
-    
+
     virtual void BeginOfEventAction(const G4Event* );
     virtual void EndOfEventAction(const G4Event* );
 
@@ -29,6 +29,7 @@ class G4BeamTestEventAction : public G4UserEventAction
   G4double  fIntegralZ;
   G4double  fXIn;
   G4int SiCollID;
+  G4int SC4CollID;
   G4int hcID;
 
 	G4double  fYIn;
@@ -37,6 +38,3 @@ class G4BeamTestEventAction : public G4UserEventAction
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 #endif
-
-    
-
diff --git a/include/G4BeamTestHadronPhysics.h.backup b/include/G4BeamTestHadronPhysics.h.backup
deleted file mode 100644
index dc2bf67..0000000
--- a/include/G4BeamTestHadronPhysics.h.backup
+++ /dev/null
@@ -1,327 +0,0 @@
-#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.backup b/include/G4BeamTestIonPhysics.h.backup
deleted file mode 100644
index 0a9d5f1..0000000
--- a/include/G4BeamTestIonPhysics.h.backup
+++ /dev/null
@@ -1,88 +0,0 @@
-#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
index 0c53595..8ebead1 100644
--- a/include/G4BeamTestMuonPhysics.h
+++ b/include/G4BeamTestMuonPhysics.h
@@ -27,10 +27,10 @@
 */
 class G4BeamTestMuonPhysics : public G4VPhysicsConstructor
 {
-public: 
+public:
   G4BeamTestMuonPhysics();
   ~G4BeamTestMuonPhysics();
-  
+
   void ConstructParticle();
   void ConstructProcess();
 
diff --git a/include/G4BeamTestPrimaryGeneratorMessenger.h b/include/G4BeamTestPrimaryGeneratorMessenger.h
index 0063f95..62181d0 100644
--- a/include/G4BeamTestPrimaryGeneratorMessenger.h
+++ b/include/G4BeamTestPrimaryGeneratorMessenger.h
@@ -15,9 +15,9 @@ class G4BeamTestPrimaryGeneratorMessenger: public G4UImessenger
   public:
     G4BeamTestPrimaryGeneratorMessenger(G4BeamTestPrimaryGeneratorAction* );
     virtual ~G4BeamTestPrimaryGeneratorMessenger();
- 
+
     virtual void SetNewValue(G4UIcommand*, G4String);
- 
+
   private:
     G4BeamTestPrimaryGeneratorAction* fG4BeamTestAction;
     G4UIdirectory*                  fGunDir;
diff --git a/include/G4BeamTestRunManager.h.backup b/include/G4BeamTestRunManager.h.backup
deleted file mode 100644
index a7eaff0..0000000
--- a/include/G4BeamTestRunManager.h.backup
+++ /dev/null
@@ -1,34 +0,0 @@
-#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 (G4BeamTestRunManager*)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/G4BeamTestSC4SD.h b/include/G4BeamTestSC4SD.h
new file mode 100644
index 0000000..d27c7e0
--- /dev/null
+++ b/include/G4BeamTestSC4SD.h
@@ -0,0 +1,42 @@
+#ifndef G4BeamTestSC4SD_h
+#define G4BeamTestSC4SD_h 1
+
+#include "G4VSensitiveDetector.hh"
+
+#include "G4BeamTestSiHit.h"
+
+#include <vector>
+
+class G4Step;
+class G4HCofThisEvent;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+/// B2Tracker sensitive detector class
+///
+/// The hits are accounted in hits in ProcessHits() function which is called
+/// by Geant4 kernel at each step. A hit is created with each step with non zero
+/// energy deposit.
+
+class G4BeamTestSC4SD : public G4VSensitiveDetector
+{
+public:
+  G4BeamTestSC4SD(const G4String& name, const G4String& hitsCollectionName);
+  virtual ~G4BeamTestSC4SD();
+
+  // methods from base class
+  virtual void   Initialize(G4HCofThisEvent* hitCollection);
+  virtual G4bool ProcessHits(G4Step* step, G4TouchableHistory* history);
+  virtual void   EndOfEvent(G4HCofThisEvent* hitCollection);
+
+
+private:
+  G4BeamTestSiHitsCollection* fHitsCollection;
+
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif
+
+
diff --git a/include/G4BeamTestSiHit.h b/include/G4BeamTestSiHit.h
index eaed902..ce093a6 100644
--- a/include/G4BeamTestSiHit.h
+++ b/include/G4BeamTestSiHit.h
@@ -9,7 +9,7 @@
 #include "G4ThreeVector.hh"
 #include "tls.hh"
 
-static std::fstream testnew("./testnew.txt", std::ofstream::out);
+#include "G4Interface.h"
 
 /// Tracker hit class
 ///
diff --git a/include/G4BeamTestSiSD.h b/include/G4BeamTestSiSD.h
index 20992c0..029175d 100644
--- a/include/G4BeamTestSiSD.h
+++ b/include/G4BeamTestSiSD.h
@@ -15,16 +15,15 @@ class G4HCofThisEvent;
 /// B2Tracker sensitive detector class
 ///
 /// The hits are accounted in hits in ProcessHits() function which is called
-/// by Geant4 kernel at each step. A hit is created with each step with non zero 
+/// by Geant4 kernel at each step. A hit is created with each step with non zero
 /// energy deposit.
 
 class G4BeamTestSiSD : public G4VSensitiveDetector
 {
 public:
-  G4BeamTestSiSD(const G4String& name,
-                const G4String& hitsCollectionName);
+  G4BeamTestSiSD(const G4String& name, const G4String& hitsCollectionName);
   virtual ~G4BeamTestSiSD();
-  
+
   // methods from base class
   virtual void   Initialize(G4HCofThisEvent* hitCollection);
   virtual G4bool ProcessHits(G4Step* step, G4TouchableHistory* history);
diff --git a/include/G4BeamTestTank.h b/include/G4BeamTestTank.h
index c4d3400..422d49a 100644
--- a/include/G4BeamTestTank.h
+++ b/include/G4BeamTestTank.h
@@ -30,7 +30,7 @@ class G4BeamTestTank
 
   ~G4BeamTestTank();
 
-  const G4ThreeVector& GetPos() const {return position_;}
+  // const G4ThreeVector& GetPos() const {return position_;}
 
   /* /// Position of center of the tank */
   /* I3Position GetPos_I3(); */
@@ -86,7 +86,7 @@ class G4BeamTestTank
   G4double glassOuterRadius_;
   G4double glassThickness_;
 
-  G4ThreeVector position_;
+  // G4ThreeVector position_;
 
   // G4ThreeVector delaunayPoint1_;
   // G4ThreeVector delaunayPoint2_;
diff --git a/include/G4Interface.h b/include/G4Interface.h
index e2956af..6fdce4b 100644
--- a/include/G4Interface.h
+++ b/include/G4Interface.h
@@ -2,8 +2,6 @@
 #define _TOPSIM_G4INTERFACE_H_
 
 #include "G4RunManager.hh"
-/* #include "G4BeamTestRunManager.h" */
-/* #include <icetray/I3Logging.h> */
 
 #ifdef G4VIS_USE
 class G4VisManager;
@@ -13,6 +11,8 @@ class I3Particle;
 class G4BeamTestTank;
 class G4BeamTestDetectorConstruction;
 
+extern std::fstream testnew;
+
 /**
  * Top-level class to handle Geant4. All global things are initialized here (run manager, visualization manager, detector construction, physics list and user actions).
  */
diff --git a/include/G4TankIceSD.h b/include/G4TankIceSD.h
deleted file mode 100644
index 25ce7ca..0000000
--- a/include/G4TankIceSD.h
+++ /dev/null
@@ -1,57 +0,0 @@
-#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. TODO(shivesh): make the PMT the SD
- *
- * 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 G4ThreeVector domPositions_;
-
-  /// Cherenkov production. See technical note
-  G4double GetCerenkovNumber(G4Step* aStep);
-  G4double GetProbability(G4double radius, G4double height);
-
-  G4double sumEdep_;
-  G4double cogTime_;
-  G4double cherenkovCounter_;
-  G4double cherenkovCounterWeight_;
-};
-
-#endif
diff --git a/src/#G4BeamTestSiHit.cxx# b/src/#G4BeamTestSiHit.cxx#
deleted file mode 100644
index 663a03f..0000000
--- a/src/#G4BeamTestSiHit.cxx#
+++ /dev/null
@@ -1,115 +0,0 @@
-#include <iostream>
-#include <iomanip>
-
-#include "G4BeamTestSiHit.h"
-#include "G4UnitsTable.hh"
-#include "G4VVisManager.hh"
-#include "G4Circle.hh"
-#include "G4Colour.hh"
-#include "G4VisAttributes.hh"
-
-G4ThreadLocal G4Allocator<G4BeamTestSiHit>* G4BeamTestSiHitAllocator=0;
-static std::fstream testnew("./testnew.txt", std::ofstream::out);
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-G4BeamTestSiHit::G4BeamTestSiHit()
-  : G4VHit(),
-    fTrackID(-1),
-  
-    fEdep(0.),
-              fTime(0.),
-    fPos(G4ThreeVector())
-{
-  /* G4cout << "opening file" << G4endl; */
-  /* testnew.open(testnew_out); */
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-G4BeamTestSiHit ::~G4BeamTestSiHit() {
-  /* G4cout << "closing file" << G4endl; */
-  /* testnew.close(); */
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-G4BeamTestSiHit::G4BeamTestSiHit(const G4BeamTestSiHit& right)
-  : G4VHit()
-{
-  fTrackID   = right.fTrackID;
-  // fChamberNb = right.fChamberNb;
-  fEdep      = right.fEdep;
-  fPos       = right.fPos;
-  fTime = right.fTime;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-const G4BeamTestSiHit& G4BeamTestSiHit::operator=(const G4BeamTestSiHit& right)
-{
-  fTrackID   = right.fTrackID;
-  //  fChamberNb = right.fChamberNb;
-  fEdep      = right.fEdep;
-  fPos       = right.fPos;
-  fTime = right.fTime;
-
-  return *this;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-G4int G4BeamTestSiHit::operator==(const G4BeamTestSiHit& right) const
-{
-  return ( this == &right ) ? 1 : 0;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-void G4BeamTestSiHit::Draw()
-{
-  /* G4VVisManager* pVVisManager = G4VVisManager::GetConcreteInstance(); */
-  /* if(pVVisManager) */
-  /*   { */
-  /*     G4Circle circle(fPos); */
-  /*     circle.SetScreenSize(4.); */
-  /*     circle.SetFillStyle(G4Circle::filled); */
-  /*     G4Colour colour(1.,0.,0.); */
-  /*     G4VisAttributes attribs(colour); */
-  /*     circle.SetVisAttributes(attribs); */
-  /*     pVVisManager->Draw(circle); */
-  /*   } */
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-
-void G4BeamTestSiHit::Print()
-{
- G4cout
-    << " trackID: " << fTrackID
-    << "Edep: "
-    << std::setw(7) << G4BestUnit(fEdep,"Energy")
-    << " Position: "
-    << std::setw(7) << G4BestUnit(fPos ,"Length")
-    << "Time: "
-    << std::setw(7) << G4BestUnit(fTime,"Time")
-    << G4endl;
-}
-
-void G4BeamTestSiHit::Dataout()
-{
-testnew
-  << " trackID: " << fTrackID
-  << "Edep: "
-  << std::setw(7) << G4BestUnit(fEdep,"Energy")
-  << " Position: "
-  << std::setw(7) << G4BestUnit( fPos ,"Length")
-  << "Time: "
-  << std::setw(7) << G4BestUnit( fTime,"Time")
-  << G4endl;
-
- }
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/src/#G4TestEMPhysics.cxx# b/src/#G4TestEMPhysics.cxx#
deleted file mode 100644
index e69de29..0000000
--- a/src/#G4TestEMPhysics.cxx#
+++ /dev/null
diff --git a/src/G4BeamTestDetectorConstruction.cxx b/src/G4BeamTestDetectorConstruction.cxx
index 495330a..2942ec6 100644
--- a/src/G4BeamTestDetectorConstruction.cxx
+++ b/src/G4BeamTestDetectorConstruction.cxx
@@ -13,6 +13,7 @@
 #include "G4BeamTestDetectorConstruction.h"
 #include "G4BeamTestTank.h"
 
+#include "G4BeamTestSC4SD.h"
 
 G4BeamTestDetectorConstruction::G4BeamTestDetectorConstruction():
   origin_(0, 0, 0), verboseLevel_(0)/* , tankList_(0) */
@@ -27,22 +28,18 @@ G4BeamTestDetectorConstruction::~G4BeamTestDetectorConstruction()
 
 G4VPhysicalVolume* G4BeamTestDetectorConstruction::Construct()
 {
-  /* if(tankList_.empty()) return NULL; */
-
   CreateMaterials();
 
-  /* // World origin in IceCube coordinates */
-  /* origin_.set(delaunay.GetOrigin().x(), delaunay.GetOrigin().y(), zSnowBottom + zHalfLength); */
-
   // Determine World dimensions
+  // G4double xWorld = 20.0 * CLHEP::m;
   G4double xWorld = 4.0 * CLHEP::m;
   G4double yWorld = 4.0 * CLHEP::m;
   G4double zWorld = 4.0 * CLHEP::m;
 
   // SC4 dimensions
-  G4double scinHeight_ = 1 * 2.54 * CLHEP::cm;
-  G4double scinWidth_ = 1 * 2.54 * CLHEP::cm;
-  G4double scinThickness_ = 1 * 2.54 * CLHEP::cm;
+  G4double scinX_ = 0.25 * 2.54 * CLHEP::cm;
+  G4double scinY_ = 1 * 2.54 * CLHEP::cm;
+  G4double scinZ_ = 1 * 2.54 * CLHEP::cm;
 
   // Create world volume
   G4Box* world_box = new G4Box("solid_world", xWorld*0.5, yWorld*0.5, zWorld*0.5);
@@ -51,39 +48,61 @@ G4VPhysicalVolume* G4BeamTestDetectorConstruction::Construct()
   G4VPhysicalVolume* worldPhys =
       new G4PVPlacement(0, G4ThreeVector(), worldLog, "world", 0, 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 tank
-  /* BOOST_FOREACH(G4BeamTestTank* tank, tankList_) */
-  /* { */
-  tank_->InstallTank(worldPhys, origin_);
-  // }
+  // Location calculations
+  G4ThreeVector mwpc1Vec = G4ThreeVector(198.8*CLHEP::cm-xWorld*0.5, 0, 0);
+  G4ThreeVector sc1Vec   = mwpc1Vec + G4ThreeVector(-50.0*CLHEP::cm, 0, 0);
+  G4ThreeVector sc2Vec   = mwpc1Vec + G4ThreeVector( 80.0*CLHEP::cm, 0, 0);
+  G4ThreeVector sc3Vec   = mwpc1Vec + G4ThreeVector(180.0*CLHEP::cm, 0, 0);
+  G4ThreeVector mwpc2Vec = mwpc1Vec + G4ThreeVector(275.3*CLHEP::cm, 0, 0);
+  G4ThreeVector mwpc3Vec = mwpc2Vec + G4ThreeVector(284.6*CLHEP::cm, 0, 0);
+  G4ThreeVector mwpc4Vec = mwpc3Vec + G4ThreeVector(771.8*CLHEP::cm, 0, 0);
+  // G4ThreeVector tankVec  = mwpc4Vec + G4ThreeVector(300.0*CLHEP::cm, 0, 0);
+  G4ThreeVector tankVec  = G4ThreeVector(0, 0, 0);
+  G4ThreeVector sc4Vec   = tankVec +  G4ThreeVector(120.0*CLHEP::cm, 0, 0);
+
+  tank_->InstallTank(worldPhys, tankVec);
+
+  // // Define SC1
+  // G4Box* sc1_box = new G4Box("sc1",scinX_*0.5, scinY_*2.0, scinZ_*2.0);
+  // G4LogicalVolume* sc1Log =
+  //     new G4LogicalVolume(sc1_box, G4Material::GetMaterial("SC4"), "log_sc1", 0, 0, 0);
+  // G4VPhysicalVolume* sc1Phys =
+  //     new G4PVPlacement(0, sc1Vec, sc1Log, "sc1", worldLog, false, 0);
+
+  // // Define SC2
+  // G4Box* sc2_box = new G4Box("sc2",scinX_*0.5, scinY_*2.0, scinZ_*2.0);
+  // G4LogicalVolume* sc2Log =
+  //     new G4LogicalVolume(sc2_box, G4Material::GetMaterial("SC4"), "log_sc2", 0, 0, 0);
+  // G4VPhysicalVolume* sc2Phys =
+  //     new G4PVPlacement(0, sc2Vec, sc2Log, "sc2", worldLog, false, 0);
+
+  // // Define SC3
+  // G4Box* sc3_box = new G4Box("sc3",scinX_*0.5, scinY_*2.0, scinZ_*2.0);
+  // G4LogicalVolume* sc3Log =
+  //     new G4LogicalVolume(sc3_box, G4Material::GetMaterial("SC4"), "log_sc3", 0, 0, 0);
+  // G4VPhysicalVolume* sc3Phys =
+  //     new G4PVPlacement(0, sc3Vec, sc3Log, "sc3", worldLog, false, 0);
 
   // Define SC4
-  G4Box* sc4_box = new G4Box("sc4",scinHeight_*0.5, scinWidth_*0.5, scinThickness_*0.5);
+  G4Box* sc4_box = new G4Box("sc4",scinX_, scinY_*0.5, scinZ_*0.5);
   G4LogicalVolume* sc4Log =
       new G4LogicalVolume(sc4_box, G4Material::GetMaterial("SC4"), "log_sc4", 0, 0, 0);
   G4VPhysicalVolume* sc4Phys =
-      new G4PVPlacement(0, G4ThreeVector(1.2*CLHEP::m,0,0), sc4Log, "sc4", worldLog, false, 0);
+      new G4PVPlacement(0, sc4Vec, sc4Log, "sc4", worldLog, false, 0);
 
   // 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 */
   energyLimit->SetUserMinEkine(2.26 * CLHEP::eV);  // Lower threshold of PMT - 550nm
-  // energyLimit->SetUserMaxEkine(3.55 * CLHEP::eV); //upper threshold of PMT - 350nm
   worldLog->SetUserLimits(energyLimit);
-  /* snowLog->SetUserLimits(energyLimit); */
+  // sc1Log->SetUserLimits(energyLimit);
+  // sc2Log->SetUserLimits(energyLimit);
+  // sc3Log->SetUserLimits(energyLimit);
   sc4Log->SetUserLimits(energyLimit);
-  
-  // G4SDManager* sdManager = G4SDManager::GetSDMpointer();
-  // sc4SD_ = new G4BeamTestSC4SD("sc4_SD_", "HitsCollection");
-  // sdManager->AddNewDetector(sc4SD_);
-  // sc4Log->SetSensitiveDetector(sc4SD_);
+
+  G4SDManager* sdManager = G4SDManager::GetSDMpointer();
+  sc4SD_ = new G4BeamTestSC4SD("sc4_SD_", "SDHitsCollection");
+  sdManager->AddNewDetector(sc4SD_);
+  sc4Log->SetSensitiveDetector(sc4SD_);
 
   return worldPhys;
 }
@@ -93,16 +112,11 @@ G4VPhysicalVolume* G4BeamTestDetectorConstruction::Construct()
 void G4BeamTestDetectorConstruction::CreateMaterials()
 {
   CreateAir();
-  /* CreateIce(); */
-  /* CreateSnow(); */
   CreateWater();
   CreatePlastic();
-  /* CreatePerlite(); */
   CreateGlassSphere();
   CreateEffectiveDOMMaterial();
   CreateSC4();
-
-  //if(verboseLevel_>0) G4cout << *G4Material::GetMaterialTable() << G4endl;
 }
 
 /*****************************************************************/
@@ -122,16 +136,86 @@ void G4BeamTestDetectorConstruction::CreateWater()
   water->AddElement(nistManager->FindOrBuildElement("H"), 2);
   water->AddElement(nistManager->FindOrBuildElement("O"), 1);
 
-  // pmt spectral response 300-650nm
-  // https://docushare.icecube.wisc.edu/dsweb/Get/Document-6637/R7081-02%20data%20sheet.pdf<Paste>
-  // TODO(shivesh): add more properties?
-  const G4int water_bins = 2;
-  // G4double water_ephot[water_bins] = {1.91 * CLHEP::eV, 4.13 * CLHEP::eV};
-  G4double water_ephot[water_bins] = {0.1 * CLHEP::eV, 10 * CLHEP::eV};
-  G4double water_refr[water_bins] = {1.33, 1.33};
+  // const G4int water_bins = 2;
+  // G4double water_ephot[water_bins] = {0.1 * CLHEP::eV, 10 * CLHEP::eV};
+  // G4double water_refr[water_bins] = {1.33, 1.33};
+
+   //From SFDETSIM water absorption
+   const G4int NUMENTRIES_water=60;
+
+   G4double ENERGY_water[NUMENTRIES_water] =
+     { 1.56962e-09*GeV, 1.58974e-09*GeV, 1.61039e-09*GeV, 1.63157e-09*GeV, 
+       1.65333e-09*GeV, 1.67567e-09*GeV, 1.69863e-09*GeV, 1.72222e-09*GeV, 
+       1.74647e-09*GeV, 1.77142e-09*GeV,1.7971e-09*GeV, 1.82352e-09*GeV, 
+       1.85074e-09*GeV, 1.87878e-09*GeV, 1.90769e-09*GeV, 1.93749e-09*GeV, 
+       1.96825e-09*GeV, 1.99999e-09*GeV, 2.03278e-09*GeV, 2.06666e-09*GeV,
+       2.10169e-09*GeV, 2.13793e-09*GeV, 2.17543e-09*GeV, 2.21428e-09*GeV, 
+       2.25454e-09*GeV, 2.29629e-09*GeV, 2.33962e-09*GeV, 2.38461e-09*GeV, 
+       2.43137e-09*GeV, 2.47999e-09*GeV, 2.53061e-09*GeV, 2.58333e-09*GeV, 
+       2.63829e-09*GeV, 2.69565e-09*GeV, 2.75555e-09*GeV, 2.81817e-09*GeV, 
+       2.88371e-09*GeV, 2.95237e-09*GeV, 3.02438e-09*GeV, 3.09999e-09*GeV,
+       3.17948e-09*GeV, 3.26315e-09*GeV, 3.35134e-09*GeV, 3.44444e-09*GeV, 
+       3.54285e-09*GeV, 3.64705e-09*GeV, 3.75757e-09*GeV, 3.87499e-09*GeV, 
+       3.99999e-09*GeV, 4.13332e-09*GeV, 4.27585e-09*GeV, 4.42856e-09*GeV, 
+       4.59258e-09*GeV, 4.76922e-09*GeV, 4.95999e-09*GeV, 5.16665e-09*GeV, 
+       5.39129e-09*GeV, 5.63635e-09*GeV, 5.90475e-09*GeV, 6.19998e-09*GeV }; 
+
+   // M Fechner : new ; define the water refraction index using refsg.F 
+   //from skdetsim using the whole range.   
+    G4double RINDEX1[NUMENTRIES_water] = 
+     {1.32885, 1.32906, 1.32927, 1.32948, 1.3297, 1.32992, 1.33014, 
+      1.33037, 1.3306, 1.33084, 1.33109, 1.33134, 1.3316, 1.33186, 1.33213,
+      1.33241, 1.3327, 1.33299, 1.33329, 1.33361, 1.33393, 1.33427, 1.33462,
+      1.33498, 1.33536, 1.33576, 1.33617, 1.3366, 1.33705, 1.33753, 1.33803,
+      1.33855, 1.33911, 1.3397, 1.34033, 1.341, 1.34172, 1.34248, 1.34331,
+      1.34419, 1.34515, 1.3462, 1.34733, 1.34858, 1.34994, 1.35145, 1.35312,
+      1.35498, 1.35707, 1.35943, 1.36211, 1.36518, 1.36872, 1.37287, 1.37776,
+      1.38362, 1.39074, 1.39956, 1.41075, 1.42535};
+
+    G4double ABWFF = 1.0;
+    //T. Akiri: Values from Skdetsim 
+    G4double ABSORPTION_water[NUMENTRIES_water] =
+      {
+        16.1419*cm*ABWFF,  18.278*cm*ABWFF, 21.0657*cm*ABWFF, 24.8568*cm*ABWFF, 30.3117*cm*ABWFF, 
+        38.8341*cm*ABWFF, 54.0231*cm*ABWFF, 81.2306*cm*ABWFF, 120.909*cm*ABWFF, 160.238*cm*ABWFF, 
+        193.771*cm*ABWFF, 215.017*cm*ABWFF, 227.747*cm*ABWFF,  243.85*cm*ABWFF, 294.036*cm*ABWFF, 
+        321.647*cm*ABWFF,  342.81*cm*ABWFF, 362.827*cm*ABWFF, 378.041*cm*ABWFF, 449.378*cm*ABWFF,
+        739.434*cm*ABWFF, 1114.23*cm*ABWFF, 1435.56*cm*ABWFF, 1611.06*cm*ABWFF, 1764.18*cm*ABWFF, 
+        2100.95*cm*ABWFF,  2292.9*cm*ABWFF, 2431.33*cm*ABWFF,  3053.6*cm*ABWFF, 4838.23*cm*ABWFF, 
+        6539.65*cm*ABWFF, 7682.63*cm*ABWFF, 9137.28*cm*ABWFF, 12220.9*cm*ABWFF, 15270.7*cm*ABWFF, 
+        19051.5*cm*ABWFF, 23671.3*cm*ABWFF, 29191.1*cm*ABWFF, 35567.9*cm*ABWFF,   42583*cm*ABWFF,
+        49779.6*cm*ABWFF, 56465.3*cm*ABWFF,   61830*cm*ABWFF, 65174.6*cm*ABWFF, 66143.7*cm*ABWFF,   
+          64820*cm*ABWFF,   61635*cm*ABWFF, 57176.2*cm*ABWFF, 52012.1*cm*ABWFF, 46595.7*cm*ABWFF, 
+        41242.1*cm*ABWFF, 36146.3*cm*ABWFF, 31415.4*cm*ABWFF, 27097.8*cm*ABWFF, 23205.7*cm*ABWFF, 
+        19730.3*cm*ABWFF, 16651.6*cm*ABWFF, 13943.6*cm*ABWFF, 11578.1*cm*ABWFF, 9526.13*cm*ABWFF
+      };
+
+   G4double RAYFF = 0.625;
+   //T. Akiri: Values from Skdetsim 
+   G4double RAYLEIGH_water[NUMENTRIES_water] = {
+      386929*cm*RAYFF,  366249*cm*RAYFF,  346398*cm*RAYFF,  327355*cm*RAYFF,  309097*cm*RAYFF,  
+      291603*cm*RAYFF,  274853*cm*RAYFF,  258825*cm*RAYFF,  243500*cm*RAYFF,  228856*cm*RAYFF,  
+      214873*cm*RAYFF,  201533*cm*RAYFF,  188816*cm*RAYFF,  176702*cm*RAYFF,  165173*cm*RAYFF,
+      154210*cm*RAYFF,  143795*cm*RAYFF,  133910*cm*RAYFF,  124537*cm*RAYFF,  115659*cm*RAYFF,  
+      107258*cm*RAYFF, 99318.2*cm*RAYFF, 91822.2*cm*RAYFF,   84754*cm*RAYFF, 78097.3*cm*RAYFF, 
+     71836.5*cm*RAYFF,   65956*cm*RAYFF, 60440.6*cm*RAYFF, 55275.4*cm*RAYFF, 50445.6*cm*RAYFF,
+       45937*cm*RAYFF, 41735.2*cm*RAYFF, 37826.6*cm*RAYFF, 34197.6*cm*RAYFF, 30834.9*cm*RAYFF, 
+     27725.4*cm*RAYFF, 24856.6*cm*RAYFF, 22215.9*cm*RAYFF, 19791.3*cm*RAYFF, 17570.9*cm*RAYFF,   
+       15543*cm*RAYFF, 13696.6*cm*RAYFF, 12020.5*cm*RAYFF, 10504.1*cm*RAYFF, 9137.15*cm*RAYFF,
+     7909.45*cm*RAYFF,  6811.3*cm*RAYFF, 5833.25*cm*RAYFF,  4966.2*cm*RAYFF, 4201.36*cm*RAYFF, 
+     3530.28*cm*RAYFF, 2944.84*cm*RAYFF, 2437.28*cm*RAYFF, 2000.18*cm*RAYFF,  1626.5*cm*RAYFF, 
+     1309.55*cm*RAYFF, 1043.03*cm*RAYFF, 821.016*cm*RAYFF,  637.97*cm*RAYFF, 488.754*cm*RAYFF
+   };
+
+  // G4MaterialPropertiesTable *mpt_water = new G4MaterialPropertiesTable ();
+  // mpt_water->AddProperty ("RINDEX", water_ephot, water_refr, water_bins);
+  // water->SetMaterialPropertiesTable(mpt_water);
   G4MaterialPropertiesTable *mpt_water = new G4MaterialPropertiesTable ();
-  mpt_water->AddProperty ("RINDEX", water_ephot, water_refr, water_bins);
+  mpt_water->AddProperty ("RINDEX", ENERGY_water, RINDEX1, NUMENTRIES_water);
+  mpt_water->AddProperty("ABSLENGTH",ENERGY_water, ABSORPTION_water, NUMENTRIES_water);
+  mpt_water->AddProperty("RAYLEIGH",ENERGY_water,RAYLEIGH_water,NUMENTRIES_water);
   water->SetMaterialPropertiesTable(mpt_water);
+
 }
 
 /*****************************************************************/
@@ -139,52 +223,32 @@ void G4BeamTestDetectorConstruction::CreateWater()
 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);
-
   // http://hypernews.slac.stanford.edu/HyperNews/geant4/get/AUX/2013/03/11/12.39-85121-chDetectorConstruction.cc
+  G4NistManager* nistManager = G4NistManager::Instance();
   // Define elements for all materials not found in the NIST database
   G4Element* Si = nistManager->FindOrBuildElement("Si");
   G4Element* B = nistManager->FindOrBuildElement("B");
   G4Element* O = nistManager->FindOrBuildElement("O");
   G4Element* Na = nistManager->FindOrBuildElement("Na");
   G4Element* Al = nistManager->FindOrBuildElement("Al");
-  G4Element* K = nistManager->FindOrBuildElement("K");  
+  G4Element* K = nistManager->FindOrBuildElement("K");
 
   G4double density;
   G4int ncomponents;
   G4double fractionmass;
   G4Material* glass = new G4Material("Glass", density= 2.23*CLHEP::g/CLHEP::cm3, ncomponents=6);
   glass->AddElement(B, fractionmass=0.040064);
-  glass->AddElement(O, fractionmass=0.539562); 
+  glass->AddElement(O, fractionmass=0.539562);
   glass->AddElement(Na, fractionmass=0.028191);
   glass->AddElement(Al, fractionmass=0.011644);
   glass->AddElement(Si, fractionmass=0.377220);
@@ -192,7 +256,6 @@ void G4BeamTestDetectorConstruction::CreateGlassSphere()
 
   // pmt spectral response 300-650nm
   // https://docushare.icecube.wisc.edu/dsweb/Get/Document-6637/R7081-02%20data%20sheet.pdf<Paste>
-  // TODO(shivesh): add more properties?
   const G4int glass_bins = 2;
   // G4double glass_ephot[glass_bins] = {1.91 * CLHEP::eV, 4.13 * CLHEP::eV};
   G4double glass_ephot[glass_bins] = {0.1 * CLHEP::eV, 10 * CLHEP::eV};
@@ -216,7 +279,7 @@ void G4BeamTestDetectorConstruction::CreateEffectiveDOMMaterial()
 
   // const G4int glass_bins = 2;
   // G4double glass_ephot[glass_bins] = {1.91 * CLHEP::eV, 4.13 * CLHEP::eV};
-  // glass 
+  // glass
   // G4double glass_ephot[glass_bins] = {0.1 * CLHEP::eV, 10 * CLHEP::eV};
   // G4double glass_refr[glass_bins] = {1.47, 1.47};
   // G4MaterialPropertiesTable *mpt_glass = new G4MaterialPropertiesTable ();
@@ -227,7 +290,6 @@ void G4BeamTestDetectorConstruction::CreateEffectiveDOMMaterial()
 void G4BeamTestDetectorConstruction::CreateSC4()
 {
   G4NistManager* nistManager = G4NistManager::Instance();
-  // POM
   G4Material* plastic = new G4Material("SC4", 1.425 * CLHEP::g / CLHEP::cm3, 3, kStateSolid);
   plastic->AddElement(nistManager->FindOrBuildElement("H"), 2);
   plastic->AddElement(nistManager->FindOrBuildElement("C"), 1);
diff --git a/src/G4BeamTestEventAction.cxx b/src/G4BeamTestEventAction.cxx
index 53402c9..bf6b09b 100644
--- a/src/G4BeamTestEventAction.cxx
+++ b/src/G4BeamTestEventAction.cxx
@@ -9,6 +9,7 @@
 #include "G4Event.hh"
 #include "G4RunManager.hh"
 
+#include "G4Interface.h"
 #include "G4BeamTestSiSD.h"
 #include "G4BeamTestSiHit.h"
 #include "G4BeamTestEventAction.h"
@@ -16,10 +17,11 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 G4BeamTestEventAction::G4BeamTestEventAction()
   : G4UserEventAction(),
- 
-    SiCollID(0)
+
+    SiCollID(0),
+    SC4CollID(0)
 {
-} 
+}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
@@ -29,75 +31,67 @@ G4BeamTestEventAction::~G4BeamTestEventAction()
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 void G4BeamTestEventAction::BeginOfEventAction(const G4Event* event )
-{    
- 
-  
+{
   G4SDManager * SDman = G4SDManager::GetSDMpointer();
   SDman->ListTree();
-  if(SiCollID<0)
-    {
-      G4String colNam;
-      SiCollID = SDman->GetCollectionID(colNam="G4BeamTestSiSDCollection");
-     
-    }
-
-
+  if(SiCollID<0) SiCollID = SDman->GetCollectionID("ice_SD_");
+  if(SC4CollID<0) SC4CollID = SDman->GetCollectionID("sc4_SD_");
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 void G4BeamTestEventAction::EndOfEventAction(const G4Event* event)
-{   
+{
 
   G4cout << ">>> Summary of Event " << event->GetEventID() << G4endl;
   // testnew << ">>> Summary of Event " << event->GetEventID() << G4endl;
 
   G4cout << SiCollID << G4endl;
+  G4cout << SC4CollID << G4endl;
 
   if(SiCollID<0) return;
+  if(SC4CollID<0) return;
 
   G4HCofThisEvent* HCE = event->GetHCofThisEvent();
   G4BeamTestSiHitsCollection* SiHC = 0;
+  G4BeamTestSiHitsCollection* SC4HC = 0;
 
-  if(HCE)
-    {
-      SiHC = (G4BeamTestSiHitsCollection*)(HCE->GetHC(SiCollID));
-
-    }
-
+  G4cout << "# collections = " << HCE->GetNumberOfCollections() << G4endl;
 
+  if(HCE) {
+      SiHC = (G4BeamTestSiHitsCollection*)(HCE->GetHC(SiCollID));
+      SC4HC = (G4BeamTestSiHitsCollection*)(HCE->GetHC(1));
+  }
 
-  if(SiHC)
+  if(SiHC && SC4HC)
     {
+      // SC4 veto
+      int sc4_hit = SC4HC->entries();
+      if (sc4_hit < 1) return;
+
       int n_hit = SiHC->entries();
-      // testnew << std::flush;
       G4cout << G4endl;
       // G4cout << "Si hits " <<
       // "--------------------------------------------------------------"
       //    << G4endl;
-      G4cout << n_hit << " hits are stored in G4BeamTestSiHitsCollection."
+      G4cout << n_hit << " hits are stored in ice_SD_"
 	     << G4endl;
       /* G4cout << "List of hits in tracker" << G4endl; */
-      // testnew << G4endl;
-      // testnew << "Si hits " <<
       // "--------------------------------------------------------------"
-      testnew << n_hit << " hits are stored in G4BeamTestSiHitsCollection."
+      testnew << n_hit << " hits are stored in ice_SD_"
 	     << G4endl;
-      /* testnew << "List of hits in tracker" << G4endl; */
       for(int i=0;i<n_hit;i++)
       {
           /* (*SiHC)[i]->Print(); */
           (*SiHC)[i]->Dataout();
       }
-
       //      G4cout << "sid + " << SiCollID << G4endl;
       //      testnew << "sid + " << SiCollID << G4endl;
       testnew << std::flush;
 
-    } 
+    }
 
 }
 
-
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
diff --git a/src/G4BeamTestHadronPhysics.cxx.backup b/src/G4BeamTestHadronPhysics.cxx.backup
deleted file mode 100644
index 1e37ed5..0000000
--- a/src/G4BeamTestHadronPhysics.cxx.backup
+++ /dev/null
@@ -1,412 +0,0 @@
-#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.backup b/src/G4BeamTestIonPhysics.cxx.backup
deleted file mode 100644
index a8cc573..0000000
--- a/src/G4BeamTestIonPhysics.cxx.backup
+++ /dev/null
@@ -1,100 +0,0 @@
-#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/G4BeamTestRunManager.cxx.backup b/src/G4BeamTestRunManager.cxx.backup
deleted file mode 100644
index 5ed2e05..0000000
--- a/src/G4BeamTestRunManager.cxx.backup
+++ /dev/null
@@ -1,131 +0,0 @@
-// On Sun, to prevent conflict with ObjectSpace, G4Timer.hh has to be
-// loaded *before* globals.hh...
-#include "G4Timer.hh"
-
-#include "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/G4BeamTestSC4SD.cxx b/src/G4BeamTestSC4SD.cxx
new file mode 100644
index 0000000..f9ef766
--- /dev/null
+++ b/src/G4BeamTestSC4SD.cxx
@@ -0,0 +1,90 @@
+
+#include <sstream>
+
+#include "G4UnitsTable.hh"
+
+#include "G4BeamTestSC4SD.h"
+#include "G4HCofThisEvent.hh"
+#include "G4Step.hh"
+#include "G4ThreeVector.hh"
+#include "G4SDManager.hh"
+#include "G4ios.hh"
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4BeamTestSC4SD::G4BeamTestSC4SD(const G4String& name, const G4String& hitsCollectionName)
+  : G4VSensitiveDetector(name)
+      , fHitsCollection(NULL)
+{
+   collectionName.insert(hitsCollectionName);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4BeamTestSC4SD::~G4BeamTestSC4SD()
+{
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4BeamTestSC4SD::Initialize(G4HCofThisEvent* hce)
+{
+  // Create hits collection
+
+  fHitsCollection
+    = new G4BeamTestSiHitsCollection(SensitiveDetectorName, collectionName[0]);
+
+  // Add this collection in hce
+
+  G4int hcID
+    = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
+  hce->AddHitsCollection( hcID, fHitsCollection );
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4bool G4BeamTestSC4SD::ProcessHits(G4Step* aStep, G4TouchableHistory*)
+{
+  //  G4double stepl = 10;
+  // aStep->SetStepLength(stepl);
+
+  G4String name=aStep->GetTrack()->GetDefinition()->GetParticleName();
+
+
+  if (name == "pi-" || name == "e-") {
+    G4cout << " Particle_name hit SC4 = " <<  name  << G4endl;
+    // total energy
+    G4double etot = aStep->GetTrack()->GetTotalEnergy();
+    // energy deposit
+    G4double edep = aStep->GetTotalEnergyDeposit();
+
+    G4BeamTestSiHit* newHit = new G4BeamTestSiHit();
+    // G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
+
+    newHit->SetTrackID  (aStep->GetTrack()->GetTrackID());
+    newHit->SetEdep(edep);
+    newHit->SetPos (aStep->GetPostStepPoint()->GetPosition());
+    newHit-> SetTime (aStep->GetPreStepPoint()->GetProperTime());
+    fHitsCollection->insert( newHit );
+
+    return true;
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4BeamTestSC4SD::EndOfEvent(G4HCofThisEvent*)
+{
+  if ( verboseLevel>1 ) {
+    G4int nofHits = fHitsCollection->entries();
+    G4cout << G4endl
+	   << "-------->Hits Collection: in this event they are " << nofHits
+	   << " hits in the tracker chambers: " << G4endl;
+      for ( G4int i=0; i<nofHits; i++ ) (*fHitsCollection)[i]->Print();
+
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
diff --git a/src/G4BeamTestSiSD.cxx b/src/G4BeamTestSiSD.cxx
index afa6bf8..194a0fc 100644
--- a/src/G4BeamTestSiSD.cxx
+++ b/src/G4BeamTestSiSD.cxx
@@ -15,8 +15,7 @@ std::ofstream testdata("EVENTDATA.txt");
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-G4BeamTestSiSD::G4BeamTestSiSD(const G4String& name,
-	   const G4String& hitsCollectionName) 
+G4BeamTestSiSD::G4BeamTestSiSD(const G4String& name, const G4String& hitsCollectionName)
   : G4VSensitiveDetector(name)
       , fHitsCollection(NULL)
 {
@@ -25,7 +24,7 @@ G4BeamTestSiSD::G4BeamTestSiSD(const G4String& name,
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-G4BeamTestSiSD::~G4BeamTestSiSD() 
+G4BeamTestSiSD::~G4BeamTestSiSD()
 { // testdata  << "Arrival time"   << "  "  << "Energy "  <<  "  "  << "Distance"  <<  std::endl;
 
 }
@@ -38,100 +37,92 @@ void G4BeamTestSiSD::Initialize(G4HCofThisEvent* hce)
 
   testdata  <<  "*"    <<  std::endl;
 
-  fHitsCollection 
-    = new G4BeamTestSiHitsCollection(SensitiveDetectorName, collectionName[0]); 
+  fHitsCollection
+    = new G4BeamTestSiHitsCollection(SensitiveDetectorName, collectionName[0]);
 
   // Add this collection in hce
 
-  G4int hcID 
+  G4int hcID
     = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
-  G4cout << "hcID " << hcID << G4endl;
-  hce->AddHitsCollection( hcID, fHitsCollection ); 
+  hce->AddHitsCollection( hcID, fHitsCollection );
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-G4bool G4BeamTestSiSD::ProcessHits(G4Step* aStep, 
-				G4TouchableHistory*)
-{  
+G4bool G4BeamTestSiSD::ProcessHits(G4Step* aStep, G4TouchableHistory*)
+{
   //  G4double stepl = 10;
   // aStep->SetStepLength(stepl);
 
   G4String name=aStep->GetTrack()->GetDefinition()->GetParticleName();
 
-  /* G4cout << " Particle_name = " <<  name  << G4endl;   */
+  // G4cout << " Particle_name = " <<  name  << G4endl;
 
   if (name == "opticalphoton" || name == "gamma") {
-  // G4cout << " Particle_name = " <<  name  << G4endl;
+    // G4cout << " Particle_name = " <<  name  << G4endl;
 
-// total energy
-  G4double etot = aStep->GetTrack()->GetTotalEnergy();
-// energy deposit
-  G4double edep = aStep->GetTotalEnergyDeposit();
+    // total energy
+    G4double etot = aStep->GetTrack()->GetTotalEnergy();
+    // energy deposit
+    G4double edep = aStep->GetTotalEnergyDeposit();
 
-  if (etot < 2.26 * CLHEP::eV) { // Lower threshold of PMT - 550nm
-  // if (etot < 2.48 * CLHEP::eV) { // Lower threshold of PMT - 500nm
-      // G4cout << "particle " << name << " under threshold with energy " << etot << G4endl;
-      return true;
-  }
-  if (etot > 3.55 * CLHEP::eV) { // Upper threshold of PMT - 350nm
-  // if (etot > 3.10 * CLHEP::eV) { // Upper threshold of PMT - 400nm
-      // G4cout << "particle " << name << " over threshold with energy " << etot << G4endl;
-      return true;
-  }
-  // G4cout << "inserting particle " << name << " with energy " << etot << " into record" << G4endl;
+    if (etot < 2.26 * CLHEP::eV) { // Lower threshold of PMT - 550nm
+    // if (etot < 2.48 * CLHEP::eV) { // Lower threshold of PMT - 500nm
+        // G4cout << "particle " << name << " under threshold with energy " << etot << G4endl;
+        return true;
+    }
+    if (etot > 3.55 * CLHEP::eV) { // Upper threshold of PMT - 350nm
+    // if (etot > 3.10 * CLHEP::eV) { // Upper threshold of PMT - 400nm
+        // G4cout << "particle " << name << " over threshold with energy " << etot << G4endl;
+        return true;
+    }
+    // G4cout << "inserting particle " << name << " with energy " << etot << " into record" << G4endl;
 
-// if (edep==0.) return false;
+    G4BeamTestSiHit* newHit = new G4BeamTestSiHit();
+    // G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
 
-  /* G4cout << " Particle_name_after_edep = " <<  name  << G4endl;   */
+    newHit->SetTrackID  (aStep->GetTrack()->GetTrackID());
+    newHit->SetEdep(edep);
+    newHit->SetPos (aStep->GetPostStepPoint()->GetPosition());
+    newHit-> SetTime (aStep->GetPreStepPoint()->GetProperTime());
+    fHitsCollection->insert( newHit );
 
-  G4BeamTestSiHit* newHit = new G4BeamTestSiHit();
-  // G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
 
-  newHit->SetTrackID  (aStep->GetTrack()->GetTrackID());
-  newHit->SetEdep(edep);
-  newHit->SetPos (aStep->GetPostStepPoint()->GetPosition());
-  newHit-> SetTime (aStep->GetPreStepPoint()->GetProperTime());
-  fHitsCollection->insert( newHit );
-  
+    G4double zbef =aStep->GetPreStepPoint()->GetPosition().z();
+    G4double zaft =aStep->GetPostStepPoint()->GetPosition().z();
+    G4double zdelta = zaft-zbef;
+    //G4double zintial = aStep
 
-  G4double zbef =aStep->GetPreStepPoint()->GetPosition().z();
-  G4double zaft =aStep->GetPostStepPoint()->GetPosition().z();
-  G4double zdelta = zaft-zbef;
-  //G4double zintial = aStep
+    G4double zfix =  (zdelta + (zbef- 4.7))*1000;
 
-  G4double zfix =  (zdelta + (zbef- 4.7))*1000;
+    //  G4double zfix =  (zdelta + (zbef- 375));
 
-  //  G4double zfix =  (zdelta + (zbef- 375));
+    //  G4cout << " Z VALUE = " <<  zfix  << G4endl;
+    G4double time = aStep->GetPreStepPoint()->GetProperTime();
+    G4double gtime = aStep->GetPreStepPoint()->GetGlobalTime();
+    //  G4cout << "Global Time " << gtime << G4endl;
+    G4double timefix = time - 0.288861775183071;
 
-  //  G4cout << " Z VALUE = " <<  zfix  << G4endl;
-  G4double time = aStep->GetPreStepPoint()->GetProperTime();
-  G4double gtime = aStep->GetPreStepPoint()->GetGlobalTime();
-  //  G4cout << "Global Time " << gtime << G4endl;
- G4double timefix = time - 0.288861775183071;
+    /*  G4cout << "TIME == " << time << G4endl; */
 
-/*  G4cout << "TIME == " << time << G4endl; */
- 
- G4double vel =  aStep->GetPreStepPoint()->GetVelocity();
- 
- G4double vel2 = zfix/time;
+     G4double vel =  aStep->GetPreStepPoint()->GetVelocity();
 
- // G4cout << G4BestUnit(vel, "Speed")<< G4endl; 
+     G4double vel2 = zfix/time;
 
-/*
+     // G4cout << G4BestUnit(vel, "Speed")<< G4endl;
 
- testdata << timefix << "  "   << edep  << "  "  <<  zfix  <<  std::endl;
-   
- h_bragg9->Fill(zfix,edep);
- enetime->Fill(timefix, edep);
- // newHit->Print();
- sally->Fill(zfix,vel);
+    /*
 
- G4cout << "DISTANCE = " << zfix << " Velocity = " << vel2 <<  " Velocityold = " << vel <<G4endl; 
+     testdata << timefix << "  "   << edep  << "  "  <<  zfix  <<  std::endl;
 
+     h_bragg9->Fill(zfix,edep);
+     enetime->Fill(timefix, edep);
+     // newHit->Print();
+     sally->Fill(zfix,vel);
 
-*/
-  return true;
+     G4cout << "DISTANCE = " << zfix << " Velocity = " << vel2 <<  " Velocityold = " << vel <<G4endl;
+    */
+    return true;
   }
 }
 
@@ -139,10 +130,10 @@ G4bool G4BeamTestSiSD::ProcessHits(G4Step* aStep,
 
 void G4BeamTestSiSD::EndOfEvent(G4HCofThisEvent*)
 {
-  if ( verboseLevel>1 ) { 
+  if ( verboseLevel>1 ) {
     G4int nofHits = fHitsCollection->entries();
     G4cout << G4endl
-	   << "-------->Hits Collection: in this event they are " << nofHits 
+	   << "-------->Hits Collection: in this event they are " << nofHits
 	   << " hits in the tracker chambers: " << G4endl;
       for ( G4int i=0; i<nofHits; i++ ) (*fHitsCollection)[i]->Print();
 
diff --git a/src/G4BeamTestTank.cxx b/src/G4BeamTestTank.cxx
index 02f6751..c4bc6b8 100644
--- a/src/G4BeamTestTank.cxx
+++ b/src/G4BeamTestTank.cxx
@@ -28,8 +28,8 @@
 G4BeamTestTank::G4BeamTestTank()
 {
   // Get tank dimensions
-  // tankThickness_ = 0.0*CLHEP::cm; // TODO(shivesh) : check thickness
-  tankThickness_ = 0.44 * 2.54 *CLHEP::cm; // TODO(shivesh) : check thickness
+  // tankThickness_ = 0.0*CLHEP::cm;
+  tankThickness_ = 0.44 * 2.54 *CLHEP::cm;
   tankHeight_  = 76.83 * 2.54 * CLHEP::cm;
   innerRadius_ = 32 * 2.54 * CLHEP::cm;
   outerRadius_ = innerRadius_ + tankThickness_;
@@ -55,112 +55,103 @@ G4VPhysicalVolume* G4BeamTestTank::InstallTank(G4VPhysicalVolume* mother, const
   // User limits (energy cutoffs)
   // Do not create photons or electrons below cherenkov threshold
   // See also corresponding UserSpecialCuts in Physicslist !!!!
-  // TODO(shivesh): Maybe do all of this as stepping action ??????
   G4UserLimits* energyLimit = new G4UserLimits();
   energyLimit->SetUserMinEkine(2.26 * CLHEP::eV);  // Lower threshold of PMT - 550nm
 
-  // std::string tankName=boost::lexical_cast<std::string>(tankKey_);
-  std::string tankName = "BTT";
-
-  // Define plastic frame TODO(shivesh): plastic or polyethelene?
   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);
+  G4Tubs* solidTank = new G4Tubs(
+          "solid_tank_",
+          0.0 * CLHEP::m, outerRadius_, 0.5 * tankHeight_,
+          0.0 * CLHEP::deg, 360.0 * CLHEP::deg
+          );
+  tankLog_ = new G4LogicalVolume(
+          solidTank, plastic, "log_tank_", 0, 0, 0
+          );
 
   // Define water volume
   G4Material* water = G4Material::GetMaterial("Water");
-  G4Tubs* solidWater = new G4Tubs(("solid_water_" + tankName).c_str(),
-                                0.0 * CLHEP::m, innerRadius_, 0.5 * fillHeight_,
-                                0.0 * CLHEP::deg, 360.0 * CLHEP::deg);
+  G4Tubs* solidWater = new G4Tubs(
+          "solid_water_",
+          0.0 * CLHEP::m, innerRadius_, 0.5 * fillHeight_,
+          0.0 * CLHEP::deg, 360.0 * CLHEP::deg
+          );
   G4LogicalVolume* logWater =
-      new G4LogicalVolume(solidWater, water, ("log_water_" + tankName).c_str(), 0, 0, 0);
-  G4ThreeVector physWaterPosition(0, 0, -0.5*tankHeight_ + tankThickness_ + 0.5*fillHeight_);
+      new G4LogicalVolume(solidWater, water, "log_water_", 0, 0, 0);
+  G4ThreeVector physWaterPosition(
+          0, 0, -0.5*tankHeight_ + tankThickness_ + 0.5*fillHeight_
+          );
   G4VPhysicalVolume* physWater ATTRIBUTE_UNUSED =
-      new G4PVPlacement(0, physWaterPosition, logWater,
-                        ("water_" + tankName).c_str(), tankLog_, false, 0);
+      new G4PVPlacement(
+              0, physWaterPosition, logWater, "water_", tankLog_, false, 0
+              );
 
   // Define air volume
   G4Material* air = G4Material::GetMaterial("Air");
-  G4cout << "airHeight_ = " << airHeight_ << G4endl;
-  G4cout << "fillHeight_ = " << fillHeight_ << G4endl;
-  G4Tubs* solidAir = new G4Tubs(("solid_air_" + tankName).c_str(),
-                                    0.0 * CLHEP::m, innerRadius_, 0.5 * airHeight_,
-                                    0.0 * CLHEP::deg, 360.0 * CLHEP::deg);
+  G4Tubs* solidAir = new G4Tubs(
+          "solid_air_",
+          0.0 * CLHEP::m, innerRadius_, 0.5 * airHeight_,
+          0.0 * CLHEP::deg, 360.0 * CLHEP::deg
+          );
   G4LogicalVolume* logAir =
-      new G4LogicalVolume(solidAir, air, ("log_air_" + tankName).c_str(), 0, 0, 0);
-  G4ThreeVector physAirPosition(0, 0, -0.5 * tankHeight_ + 0.5 * CLHEP::cm + fillHeight_ +
-                                    0.5 * airHeight_);
-  G4cout << "physAirPosition = " << physAirPosition << G4endl;
+      new G4LogicalVolume(solidAir, air, "log_air_", 0, 0, 0);
+  G4ThreeVector physAirPosition(0, 0,
+          -0.5 * tankHeight_ + 0.5 * CLHEP::cm + fillHeight_ + 0.5 * airHeight_
+          );
   G4VPhysicalVolume* physAir_UNUSED =
-      new G4PVPlacement(0, physAirPosition, logAir,
-                        ("air_" + tankName).c_str(), tankLog_, false, 0);
+      new G4PVPlacement(
+              0, physAirPosition, logAir, "air_", 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(0, 0, -0.5 * airHeight_);
   G4ThreeVector lowerDOMpos(0, 0,  0.5 * fillHeight_);
-  G4cout << "upperDOMpos = " << upperDOMpos << G4endl;
-  G4cout << "lowerDOMpos = " << lowerDOMpos << G4endl;
-
-  // domPosIce[omKey] = lowerDOMpos;
-
-  // std::string omName=boost::lexical_cast<std::string>(omKey);
-  std::string omName="BTOM";
 
-  G4Sphere *upperglasssphere = new G4Sphere (("solid_dom_up_" + omName).c_str(),
+  G4Sphere *upperglasssphere = new G4Sphere ("solid_dom_up_",
                                              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(),
+  G4Sphere *lowerglasssphere = new G4Sphere ("solid_dom_lo_",
                                              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(),
+  G4Sphere *upperdomsphere = new G4Sphere ("solid_inside_dom_up_",
                                            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(),
+  G4Sphere *lowerdomsphere = new G4Sphere ("solid_inside_dom_lo_",
                                            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);
+                          "log_dom_up_", 0, 0, 0);
   G4LogicalVolume* logLowerGlass =
       new G4LogicalVolume(lowerglasssphere, glass,
-                          ("log_dom_lo_" + omName).c_str(), 0, 0, 0);
+                          "log_dom_lo_", 0, 0, 0);
   G4LogicalVolume* logUpperDOM =
       new G4LogicalVolume(upperdomsphere, effectiveDOM,
-                          ("log_inside_dom_up_" + omName).c_str(), 0, 0, 0);
+                          "log_inside_dom_up_", 0, 0, 0);
   G4LogicalVolume* logLowerDOM =
       new G4LogicalVolume(lowerdomsphere, effectiveDOM,
-                          ("log_inside_dom_lo_" + omName).c_str(), 0, 0, 0);
+                          "log_inside_dom_lo_", 0, 0, 0);
   G4VPhysicalVolume* physUpperGlass ATTRIBUTE_UNUSED =
       new G4PVPlacement(0, upperDOMpos, logUpperGlass,
-                        ("dom_up_" + omName).c_str(), logAir, false, 0);
+                        "dom_up_", logAir, false, 0);
   G4VPhysicalVolume* physLowerGlass ATTRIBUTE_UNUSED =
       new G4PVPlacement(0, lowerDOMpos, logLowerGlass,
-                        ("dom_lo_" + omName).c_str(), logWater, false, 0);
+                        "dom_lo_", logWater, false, 0);
   G4VPhysicalVolume* physUpperDOM ATTRIBUTE_UNUSED =
       new G4PVPlacement(0, G4ThreeVector(0,0,0), logUpperDOM,
-                        ("inside_dom_up_" + omName).c_str(), logUpperGlass, false, 0);
+                        "inside_dom_up_", logUpperGlass, false, 0);
   G4VPhysicalVolume* physLowerDOM ATTRIBUTE_UNUSED =
       new G4PVPlacement(0, G4ThreeVector(0,0,0), logLowerDOM,
-                        ("inside_dom_lo_" + omName).c_str(), logLowerGlass, false, 0);
+                        "inside_dom_lo_", logLowerGlass, false, 0);
 
   // apply energy limits
   logUpperGlass->SetUserLimits(energyLimit);
@@ -171,24 +162,21 @@ G4VPhysicalVolume* G4BeamTestTank::InstallTank(G4VPhysicalVolume* mother, const
 
   // Define sensitive detector
   G4SDManager* sdManager = G4SDManager::GetSDMpointer();
-  iceSD_ = new G4BeamTestSiSD(("ice_SD_" + tankName).c_str(), "HitsCollection");
+  iceSD_ = new G4BeamTestSiSD("ice_SD_", "HitsCollection");
   sdManager->AddNewDetector(iceSD_);
   // logLowerDOM->SetSensitiveDetector(iceSD_);
   logLowerGlass->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();
-  G4cout << "tankPos = " << tankPos << G4endl;
-
-  G4VPhysicalVolume* tankPhys = new G4PVPlacement(0, tankPos, tankLog_,
-                                                  ("tank_" + tankName).c_str(),
-                                                  mother->GetLogicalVolume(), false, 0);
+  // G4ThreeVector tankPos = position_ - origin - mother->GetTranslation();
+  G4ThreeVector tankPos = origin;
+  G4VPhysicalVolume* tankPhys = new G4PVPlacement(
+          0, tankPos, tankLog_, "tank_", mother->GetLogicalVolume(), false, 0
+          );
 
   // apply energy limits
   tankLog_->SetUserLimits(energyLimit);
@@ -197,84 +185,3 @@ G4VPhysicalVolume* G4BeamTestTank::InstallTank(G4VPhysicalVolume* mother, const
 
   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) */
-/* } */
-/*  */
-/*  */
-/* double G4BeamTestTank::GetTankRadius_I3() */
-/* { */
-/*   return (outerRadius_ / CLHEP::m) */
-/* } */
-/*  */
-/*  */
-/* double G4BeamTestTank::GetFillHeight_I3() */
-/* { */
-/*   return (fillHeight_ / CLHEP::m) */
-/* } */
diff --git a/src/G4BeamTestUserSteppingAction.cxx b/src/G4BeamTestUserSteppingAction.cxx
index f81ab62..816dc4a 100644
--- a/src/G4BeamTestUserSteppingAction.cxx
+++ b/src/G4BeamTestUserSteppingAction.cxx
@@ -40,7 +40,7 @@ void G4BeamTestUserSteppingAction::UserSteppingAction(const G4Step* step)
       //check if particle energy is below threshold; if true, kill the particle
       G4double energy = track->GetTotalEnergy();
       if(energy < threshold){
-          G4cout << "SteppingAction: killing particle " << particle << " with energy " << energy << " < " << threshold << G4endl;
+          if (energy > 0) G4cout << "SteppingAction: killing particle " << particle << " with energy " << energy << " < " << threshold << G4endl;
           track->SetTrackStatus(fStopAndKill);
       }
     }
diff --git a/src/G4BeamTestUserTrackingAction.cxx b/src/G4BeamTestUserTrackingAction.cxx
index e41ccc3..09a7e59 100644
--- a/src/G4BeamTestUserTrackingAction.cxx
+++ b/src/G4BeamTestUserTrackingAction.cxx
@@ -11,36 +11,35 @@ 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);
-  G4double max_threshold = 3.54;
-  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
-        G4String particle = (*secondaries)[i]->GetDefinition()->GetParticleName();
-        if(particle == "gamma" || particle == "opticalphoton")
-        {
-          //check if particle energy is below threshold; if true, kill the particle
-          G4double energy = (*secondaries)[i]->GetTotalEnergy();
-          // if(energy < threshold){
-          //     G4cout << "TrackingAction: killing particle " << particle << " with energy " << energy << " < " << threshold << G4endl;
-          //     (*secondaries)[i]->SetTrackStatus(fStopAndKill);
-          // }
+  // 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);
+  // G4double max_threshold = 3.54;
+  // 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
+  //       G4String particle = (*secondaries)[i]->GetDefinition()->GetParticleName();
+  //       if(particle == "gamma" || particle == "opticalphoton")
+  //       {
+  //         //check if particle energy is below threshold; if true, kill the particle
+  //         G4double energy = (*secondaries)[i]->GetTotalEnergy();
+  //         if(energy < threshold){
+  //             G4cout << "TrackingAction: killing particle " << particle << " with energy " << energy << " < " << threshold << G4endl;
+  //             (*secondaries)[i]->SetTrackStatus(fStopAndKill);
+  //         }
 	  // if (energy >  max_threshold * CLHEP::eV){
 	      // G4cout << "TrackingAction: killing particle " << particle << " with energy " << energy << " > " << max_threshold << G4endl;
-	      // (*secondaries)[i]->SetTrackStatus(fStopAndKill); 
-	    
+	      // (*secondaries)[i]->SetTrackStatus(fStopAndKill);
 	  // }
-        }
-      }
-    }
-  }
+  //       }
+  //     }
+  //   }
+  // }
 }
diff --git a/src/G4TankIceSD.cxx b/src/G4TankIceSD.cxx
deleted file mode 100644
index 2974837..0000000
--- a/src/G4TankIceSD.cxx
+++ /dev/null
@@ -1,166 +0,0 @@
-#include "G4TankIceSD.h"
-#include "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_ = 0;
-  cogTime_ = 0;
-  cherenkovCounter_ = 0;
-  cherenkovCounterWeight_ = 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(domPositions_.x() - localPosition.x(), 2) +
-                         pow(domPositions_.y() - localPosition.y(), 2));
-  G4double height = (domPositions_.z() - localPosition.z());
-
-  sumEdep_ += edep;
-  cogTime_ += edep*time;
-  cherenkovCounterWeight_ += GetProbability(radius, height) * cherenkovNumber;
-  cherenkovCounter_ += 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_>0)
-  {
-    cogTime_ /= sumEdep_;
-  }
-  // }
-}
-
-
-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); // TODO(shivesh): check this
-  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!"); */
-    G4cout << "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;
-}
-
-
-// TODO(shivesh): what is this
-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);
-}