Tue 20 Nov 17:26:02 GMT 2018

23 files changed, 922 insertions, 133 deletions

diff --git a/include/G4BeamTestDetectorConstruction.h b/include/G4BeamTestDetectorConstruction.h
index a3c0c9f..8b45fde 100644
--- a/include/G4BeamTestDetectorConstruction.h
+++ b/include/G4BeamTestDetectorConstruction.h
@@ -4,7 +4,7 @@
 #include <G4VUserDetectorConstruction.hh>
 #include <G4ThreeVector.hh>
 
-#include "G4BeamTestTank.h";
+#include "G4BeamTestTank.h"
 
 class G4BeamTestDetectorConstruction: public G4VUserDetectorConstruction
 {
diff --git a/include/G4BeamTestEMPhysics.h b/include/G4BeamTestEMPhysics.h
index e7a88dd..bf20b62 100644
--- a/include/G4BeamTestEMPhysics.h
+++ b/include/G4BeamTestEMPhysics.h
@@ -11,17 +11,6 @@
 #include <G4eBremsstrahlung.hh>
 #include <G4eplusAnnihilation.hh>
 
-#include <G4Cerenkov.hh>
-#include <G4Scintillation.hh>
-#include <G4OpAbsorption.hh>
-#include <G4OpRayleigh.hh>
-#include <G4OpMieHG.hh>
-#include <G4OpBoundaryProcess.hh>
-
-#include <G4LossTableManager.hh>
-#include <G4EmSaturation.hh>
-
-
 /**
    @class G4BeamTestEMPhysics
    @brief Electromagnetic physics. Used only if Geant4 version is earlier than 4.10.
@@ -41,7 +30,6 @@ public:
 
   void ConstructParticle();
   void ConstructProcess();
-  void ConstructOp();
 
 private:
   // Gamma physics
@@ -59,9 +47,6 @@ private:
   G4eIonisation positronIonisation; 
   G4eBremsstrahlung positronBremsStrahlung;  
   G4eplusAnnihilation annihilation;
-
-  // Cerenkov physics
-  G4Cerenkov cerenkov;
 };
 
 #endif // G4TANKRESPONSE_G4BEAMTESTEMPHYSICS_H_INCLUDED
diff --git a/include/G4BeamTestEventAction.h b/include/G4BeamTestEventAction.h
new file mode 100644
index 0000000..2d6de49
--- /dev/null
+++ b/include/G4BeamTestEventAction.h
@@ -0,0 +1,42 @@
+#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 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;
+  G4int SiCollID;
+  G4int hcID;
+
+	G4double  fYIn;
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif
+
+    
+
diff --git a/include/G4BeamTestRunManager.h b/include/G4BeamTestRunManager.h
new file mode 100644
index 0000000..a7eaff0
--- /dev/null
+++ b/include/G4BeamTestRunManager.h
@@ -0,0 +1,34 @@
+#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/G4BeamTestSiHit.h b/include/G4BeamTestSiHit.h
new file mode 100644
index 0000000..44dd69a
--- /dev/null
+++ b/include/G4BeamTestSiHit.h
@@ -0,0 +1,104 @@
+#ifndef G4BeamTestSiHit_h
+#define G4BeamTestSiHit_h 1
+#include <iostream>
+#include <fstream>
+
+#include "G4VHit.hh"
+#include "G4THitsCollection.hh"
+#include "G4Allocator.hh"
+#include "G4ThreeVector.hh"
+#include "tls.hh"
+
+/// Tracker hit class
+///
+/// It defines data members to store the trackID, chamberNb, energy deposit,
+/// and position of charged particles in a selected volume:
+/// - fTrackID, fChamberNB, fEdep, fPos
+
+class G4BeamTestSiHit : public G4VHit
+{
+public:
+  G4BeamTestSiHit();
+  G4BeamTestSiHit(const G4BeamTestSiHit&);
+  virtual ~G4BeamTestSiHit();
+
+  // operators
+  const G4BeamTestSiHit& operator=(const G4BeamTestSiHit&);
+  G4int operator==(const G4BeamTestSiHit&) const;
+
+  inline void* operator new(size_t);
+  inline void  operator delete(void*);
+
+  // methods from base class
+  virtual void Draw();
+  virtual void Print();
+  virtual void Dataout();
+
+  // Set methods
+  void SetTrackID  (G4int track)      { fTrackID = track; };
+  void SetEdep     (G4double de)      { fEdep = de; };
+  void SetPos      (G4ThreeVector xyz){ fPos = xyz; };
+  void SetTime(G4double val) { fTime = val; };
+
+
+  // Get methods
+  G4int GetTrackID() const     { return fTrackID; };
+  G4double GetEdep() const     { return fEdep; };
+  G4ThreeVector GetPos() const { return fPos; };
+  G4double GetTime() const { return fTime; }
+
+  /* const std::string testnew_out = "./testnew.txt"; */
+
+private:
+
+  G4ThreeVector position;
+  G4ThreeVector momentum;
+  G4double      energy;
+  G4double time;
+  // std::ofstream testnew;
+
+
+public:
+  inline void SetPosition(G4ThreeVector pos)
+  { position=pos; }
+  inline G4ThreeVector GetPosition()
+  { return position; }
+  inline void SetEnergy(G4double ene)
+  { energy = ene; }
+  inline G4double GetEnergy()
+  { return energy; }
+ 
+
+private:
+  G4double    fTime;
+  G4int         fTrackID;
+  G4double      fEdep;
+  G4ThreeVector fPos;
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+typedef G4THitsCollection<G4BeamTestSiHit> G4BeamTestSiHitsCollection;
+
+extern G4ThreadLocal G4Allocator<G4BeamTestSiHit>* G4BeamTestSiHitAllocator;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+inline void* G4BeamTestSiHit::operator new(size_t)
+{
+  if(!G4BeamTestSiHitAllocator)
+    G4BeamTestSiHitAllocator = new G4Allocator<G4BeamTestSiHit>;
+  return (void *) G4BeamTestSiHitAllocator->MallocSingle();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+inline void G4BeamTestSiHit::operator delete(void *hit)
+{
+  G4BeamTestSiHitAllocator->FreeSingle((G4BeamTestSiHit*) hit);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif
+
diff --git a/include/G4BeamTestSiSD.h b/include/G4BeamTestSiSD.h
new file mode 100644
index 0000000..20992c0
--- /dev/null
+++ b/include/G4BeamTestSiSD.h
@@ -0,0 +1,42 @@
+#ifndef G4BeamTestSiSD_h
+#define G4BeamTestSiSD_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 G4BeamTestSiSD : public G4VSensitiveDetector
+{
+public:
+  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);
+  virtual void   EndOfEvent(G4HCofThisEvent* hitCollection);
+
+
+private:
+  G4BeamTestSiHitsCollection* fHitsCollection;
+
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif
+
diff --git a/include/G4BeamTestTank.h b/include/G4BeamTestTank.h
index 3f2fd8c..c4d3400 100644
--- a/include/G4BeamTestTank.h
+++ b/include/G4BeamTestTank.h
@@ -8,12 +8,13 @@
 
 class G4LogicalVolume;
 class G4VPhysicalVolume;
-class G4TankIceSD;
 
 #include <G4Types.hh>
 #include <G4ThreeVector.hh>
 #include <map>
 
+#include "G4BeamTestSiSD.h"
+
 /**
  * This class constructs the physical volume of single tanks, serves as bridge to the sensitive detector in the tank, and computes the location of the point that define the snow surface. It also stores the basic tank dimensions (height, radius, snow).
  *
@@ -53,9 +54,9 @@ class G4BeamTestTank
   /* double   GetSnowHeight_I3(); */
 
   /// Energy deposit for a given OM, in Geant4 units (from G4TankIceSD)
-  G4double GetEDep_G4(/* const OMKey& omKey */);
+  /* G4double GetEDep_G4(#<{(| const OMKey& omKey |)}>#); */
   /// Emission time for a given OM, in Geant4 units (from G4TankIceSD)
-  G4double GetTime_G4(/* const OMKey& omKey */);
+  /* G4double GetTime_G4(#<{(| const OMKey& omKey |)}>#); */
 
   /* /// Energy deposit for a given OM, in Geant4 units (from G4TankIceSD) */
   /* double   GetEDep_I3(const OMKey& omKey); */
@@ -63,9 +64,9 @@ class G4BeamTestTank
   /* double   GetTime_I3(const OMKey& omKey); */
 
   /// Number of Cherenkovs for a given OM (from G4TankIceSD)
-  double GetNumCherenkov(/* const OMKey& omKey */);
+  /* double GetNumCherenkov(#<{(| const OMKey& omKey |)}>#); */
   /// Number of Cherenkovs weighted by emission point (from G4TankIceSD)
-  double GetNumCherenkovWeight(/* const OMKey& omKey */);
+  /* double GetNumCherenkovWeight(#<{(| const OMKey& omKey |)}>#); */
 
   /* const TankKey& GetTankKey() const {return tankKey_;} */
 
@@ -94,7 +95,7 @@ class G4BeamTestTank
   // std::map<OMKey, G4ThreeVector> relDomPositions_;
 
   G4LogicalVolume* tankLog_;
-  G4TankIceSD* iceSD_;
+  G4BeamTestSiSD* iceSD_;
 
   /* const TankKey& tankKey_; */
   /* const I3Geometry& geometry_; */
diff --git a/include/G4BeamTestUserStackingAction.h b/include/G4BeamTestUserStackingAction.h
new file mode 100644
index 0000000..556f3d1
--- /dev/null
+++ b/include/G4BeamTestUserStackingAction.h
@@ -0,0 +1,28 @@
+#ifndef G4BeamTestUserStackingAction_H
+#define G4BeamTestUserStackingAction_H 1
+
+#include "globals.hh"
+#include "G4UserStackingAction.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+class G4BeamTestUserStackingAction : public G4UserStackingAction
+{
+  public:
+    G4BeamTestUserStackingAction();
+    virtual ~G4BeamTestUserStackingAction();
+
+  public:
+    virtual G4ClassificationOfNewTrack ClassifyNewTrack(const G4Track* aTrack);
+    virtual void NewStage();
+    virtual void PrepareNewEvent();
+
+  private:
+    G4int fScintillationCounter;
+    G4int fCerenkovCounter;
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif
+
diff --git a/include/G4BeamTestUserSteppingAction.h b/include/G4BeamTestUserSteppingAction.h
index 2d09799..c4cc863 100644
--- a/include/G4BeamTestUserSteppingAction.h
+++ b/include/G4BeamTestUserSteppingAction.h
@@ -1,6 +1,7 @@
 #ifndef G4BEAMTESTUSERSTEPPINGACTION_H_INCLUDED
 #define G4BEAMTESTUSERSTEPPINGACTION_H_INCLUDED
 #include "G4UserSteppingAction.hh"
+#include "globals.hh"
 
 /**
  * Implementation of G4UserSteppingAction. This class kills gammas below threshold (set by G4BeamTestTank).
@@ -12,6 +13,11 @@ class G4BeamTestUserSteppingAction : public G4UserSteppingAction {
    ~G4BeamTestUserSteppingAction() {}
 
     void UserSteppingAction(const G4Step*);
+
+  private:
+    G4int fScintillationCounter;
+    G4int fCerenkovCounter;
+    G4int fEventNumber;
 };
 
 #endif  // G4BEAMTESTUSERSTEPPINGACTION_H_INCLUDED
diff --git a/run1.mac b/run1.mac
index 7f9c02a..b6fe271 100644
--- a/run1.mac
+++ b/run1.mac
@@ -1,3 +1,5 @@
+/vis/ogl/set/displayListLimit 1000000
+
 /gps/particle mu-
 /gps/energy 8 GeV
 
diff --git a/src/#G4TestEMPhysics.cxx# b/src/#G4TestEMPhysics.cxx#
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/src/#G4TestEMPhysics.cxx#
diff --git a/src/G4BeamTestDetectorConstruction.cxx b/src/G4BeamTestDetectorConstruction.cxx
index 95a2071..ebbd9b8 100644
--- a/src/G4BeamTestDetectorConstruction.cxx
+++ b/src/G4BeamTestDetectorConstruction.cxx
@@ -60,7 +60,8 @@ G4VPhysicalVolume* G4BeamTestDetectorConstruction::Construct()
   // Do not create photons or electrons below cherenkov threshold
   // See also corresponding UserSpecialCuts in Physicslist !!!!
   G4UserLimits* energyLimit = new G4UserLimits();
-  energyLimit->SetUserMinEkine(280.0 * CLHEP::keV);  // Cherenkov threshold of electrons in ice TODO(shivesh)
+  /* energyLimit->SetUserMinEkine(280.0 * CLHEP::keV);  // Cherenkov threshold of electrons in ice */
+  energyLimit->SetUserMinEkine(1.907 * CLHEP::eV);  // Lower threshold of PMT - 600nm
   worldLog->SetUserLimits(energyLimit);
   /* snowLog->SetUserLimits(energyLimit); */
   
diff --git a/src/G4BeamTestEMPhysics.cxx b/src/G4BeamTestEMPhysics.cxx
index 1aa9a7c..55ecf99 100644
--- a/src/G4BeamTestEMPhysics.cxx
+++ b/src/G4BeamTestEMPhysics.cxx
@@ -55,66 +55,4 @@ void G4BeamTestEMPhysics::ConstructProcess()
   pManager->AddProcess(&positronIonisation,         -1, 2, 2);
   pManager->AddProcess(&positronBremsStrahlung,     -1, 3, 3);  
   pManager->AddProcess(&annihilation,                0,-1, 4);  
-
-  // Cerenkov Physics
-  ConstructOp();
-}
-
-void G4BeamTestEMPhysics::ConstructOp()
-{
-  G4Cerenkov* cerenkovProcess = &cerenkov;
-  G4int fMaxNumPhotonStep = 20;
-
-  cerenkovProcess->SetMaxNumPhotonsPerStep(fMaxNumPhotonStep);
-  G4cout << "MaxNumPhoton" << fMaxNumPhotonStep << G4endl; 
-
-  cerenkovProcess->SetMaxBetaChangePerStep(10.0);
-  cerenkovProcess->SetTrackSecondariesFirst(true);
- G4Scintillation* scintillationProcess = new G4Scintillation("Scintillation");
- scintillationProcess->SetScintillationYieldFactor(1.);
- scintillationProcess->SetTrackSecondariesFirst(true);
-  G4OpAbsorption* absorptionProcess = new G4OpAbsorption();
-  G4OpRayleigh* rayleighScatteringProcess = new G4OpRayleigh();
-  G4OpMieHG* mieHGScatteringProcess = new G4OpMieHG();
-  G4OpBoundaryProcess* boundaryProcess = new G4OpBoundaryProcess();
-
-/*   cerenkovProcess->SetVerboseLevel(fVerboseLebel); */
-/*  // scintillationProcess->SetVerboseLevel(fVerboseLebel); */
-/*   absorptionProcess->SetVerboseLevel(fVerboseLebel); */
-/*   rayleighScatteringProcess->SetVerboseLevel(fVerboseLebel); */
-/*   mieHGScatteringProcess->SetVerboseLevel(fVerboseLebel); */
-/*   boundaryProcess->SetVerboseLevel(fVerboseLebel); */
-  
-  // Use Birks Correction in the Scintillation process
-
- if(!G4Threading::IsWorkerThread())
- {
-   G4EmSaturation* emSaturation =
-             G4LossTableManager::Instance()->EmSaturation();
-     scintillationProcess->AddSaturation(emSaturation);
- }
-
-  auto theParticleIterator = GetParticleIterator();
-  theParticleIterator->reset();
-  while( (*theParticleIterator)() ){
-    G4ParticleDefinition* particle = theParticleIterator->value();
-    G4ProcessManager* pmanager = particle->GetProcessManager();
-    G4String particleName = particle->GetParticleName();
-    if (cerenkovProcess->IsApplicable(*particle)) {
-      pmanager->AddProcess(cerenkovProcess);
-      pmanager->SetProcessOrdering(cerenkovProcess,idxPostStep);
-    }
-   if (scintillationProcess->IsApplicable(*particle)) {
-     pmanager->AddProcess(scintillationProcess);
-     pmanager->SetProcessOrderingToLast(scintillationProcess, idxAtRest);
-     pmanager->SetProcessOrderingToLast(scintillationProcess, idxPostStep);
-   }
-    if (particleName == "opticalphoton") {
-      G4cout << " AddDiscreteProcess to OpticalPhoton " << G4endl;
-      pmanager->AddDiscreteProcess(absorptionProcess);
-      pmanager->AddDiscreteProcess(rayleighScatteringProcess);
-      pmanager->AddDiscreteProcess(mieHGScatteringProcess);
-      pmanager->AddDiscreteProcess(boundaryProcess);
-    }
-  }
 }
diff --git a/src/G4BeamTestEventAction.cxx b/src/G4BeamTestEventAction.cxx
new file mode 100644
index 0000000..176991d
--- /dev/null
+++ b/src/G4BeamTestEventAction.cxx
@@ -0,0 +1,91 @@
+#include "G4Event.hh"
+
+#include "G4HCofThisEvent.hh"
+#include "G4VHitsCollection.hh"
+#include "G4SDManager.hh"
+#include "G4ios.hh"
+#include "G4VDigitizerModule.hh"
+#include "G4DigiManager.hh"
+#include "G4Event.hh"
+#include "G4RunManager.hh"
+
+#include "G4BeamTestSiSD.h"
+#include "G4BeamTestSiHit.h"
+#include "G4BeamTestEventAction.h"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4BeamTestEventAction::G4BeamTestEventAction()
+  : G4UserEventAction(),
+ 
+    SiCollID(0)
+{
+} 
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4BeamTestEventAction::~G4BeamTestEventAction()
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4BeamTestEventAction::BeginOfEventAction(const G4Event* event )
+{    
+ 
+  
+  G4SDManager * SDman = G4SDManager::GetSDMpointer();
+  if(SiCollID<0)
+    {
+      G4String colNam;
+      SiCollID = SDman->GetCollectionID(colNam="G4BeamTestSiSDCollection");
+     
+    }
+  
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4BeamTestEventAction::EndOfEventAction(const G4Event* event)
+{   
+
+  G4cout << ">>> Summary of Event " << event->GetEventID() << G4endl;
+
+  if(SiCollID<0) return;
+
+  G4HCofThisEvent* HCE = event->GetHCofThisEvent();
+  G4BeamTestSiHitsCollection* SiHC = 0;
+
+  if(HCE)
+    {
+      SiHC = (G4BeamTestSiHitsCollection*)(HCE->GetHC(SiCollID));
+
+    }
+
+
+
+
+  if(SiHC)
+    {
+      int n_hit = SiHC->entries();
+      G4cout << G4endl;
+      G4cout << "Si hits " <<
+      "--------------------------------------------------------------"
+	     << G4endl;
+      G4cout << n_hit << " hits are stored in G4BeamTestSiHitsCollection."
+	     << G4endl;
+      /* G4cout << "List of hits in tracker" << G4endl; */
+      for(int i=0;i<n_hit;i++)
+      {
+          /* (*SiHC)[i]->Print(); */
+          (*SiHC)[i]->Dataout();
+      }
+
+      //      G4cout << "sid + " << SiCollID << G4endl;
+   
+    } 
+
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
diff --git a/src/G4BeamTestPhysicsList.cxx b/src/G4BeamTestPhysicsList.cxx
index 0c27b0b..65c575c 100644
--- a/src/G4BeamTestPhysicsList.cxx
+++ b/src/G4BeamTestPhysicsList.cxx
@@ -2,12 +2,12 @@
 #include <G4Version.hh>
 #include "G4BeamTestPhysicsList.h"
 #include "G4BeamTestGeneralPhysics.h"
-#if G4VERSION_NUMBER < 1000
-#include "G4BeamTestEMPhysics.h"
-#include "G4BeamTestMuonPhysics.h"
-#include "G4BeamTestHadronPhysics.h"
-#include "G4BeamTestIonPhysics.h"
-#else
+/* #if G4VERSION_NUMBER < 1000 */
+/* #include "G4BeamTestEMPhysics.h" */
+/* #include "G4BeamTestMuonPhysics.h" */
+/* #include "G4BeamTestHadronPhysics.h" */
+/* #include "G4BeamTestIonPhysics.h" */
+/* #else */
 #include "G4DecayPhysics.hh"
 #include "G4EmStandardPhysics.hh"
 #include "G4EmExtraPhysics.hh"
@@ -18,7 +18,9 @@
 #include "G4DataQuestionaire.hh"
 #include "G4HadronPhysicsFTFP_BERT.hh"
 #include <FTFP_BERT.hh>
-#endif
+/* #endif */
+
+#include <G4OpticalPhysics.hh>
 
 #include <G4ProcessManager.hh>
 #include <G4ParticleTypes.hh>
@@ -32,15 +34,17 @@ G4BeamTestPhysicsList::G4BeamTestPhysicsList()
   SetVerboseLevel(1);
 
   RegisterPhysics(new G4BeamTestGeneralPhysics);
-#if G4VERSION_NUMBER < 1000
-  RegisterPhysics(new G4BeamTestEMPhysics);
-  RegisterPhysics(new G4BeamTestMuonPhysics);
-  RegisterPhysics(new G4BeamTestHadronPhysics);
-  RegisterPhysics(new G4BeamTestIonPhysics);
-#else
+/* #if G4VERSION_NUMBER < 1000 */
+/*   RegisterPhysics(new G4BeamTestEMPhysics); */
+/*   RegisterPhysics(new G4BeamTestMuonPhysics); */
+/*   RegisterPhysics(new G4BeamTestHadronPhysics); */
+/*   RegisterPhysics(new G4BeamTestIonPhysics); */
+/* #else */
   // The following is basically Geant4's FTFP_BERT physics list
   G4DataQuestionaire it(photon); // this checks that G4LEVELGAMMADATA is there
 
+  RegisterPhysics(new G4OpticalPhysics());
+
   RegisterPhysics(new G4EmStandardPhysics());
   RegisterPhysics(new G4EmExtraPhysics());
   RegisterPhysics(new G4DecayPhysics());
@@ -49,7 +53,7 @@ G4BeamTestPhysicsList::G4BeamTestPhysicsList()
   RegisterPhysics(new G4StoppingPhysics());
   RegisterPhysics(new G4IonPhysics());
   //RegisterPhysics(new G4NeutronTrackingCut());
-#endif
+/* #endif */
 }
 
 
diff --git a/src/G4BeamTestRunManager.cxx b/src/G4BeamTestRunManager.cxx
new file mode 100644
index 0000000..5ed2e05
--- /dev/null
+++ b/src/G4BeamTestRunManager.cxx
@@ -0,0 +1,131 @@
+// 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/G4BeamTestSiHit.cxx b/src/G4BeamTestSiHit.cxx
new file mode 100644
index 0000000..591274c
--- /dev/null
+++ b/src/G4BeamTestSiHit.cxx
@@ -0,0 +1,116 @@
+#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/G4BeamTestSiSD.cxx b/src/G4BeamTestSiSD.cxx
new file mode 100644
index 0000000..b67092b
--- /dev/null
+++ b/src/G4BeamTestSiSD.cxx
@@ -0,0 +1,137 @@
+
+#include <sstream>
+
+#include "G4UnitsTable.hh"
+
+#include "G4BeamTestSiSD.h"
+#include "G4HCofThisEvent.hh"
+#include "G4Step.hh"
+#include "G4ThreeVector.hh"
+#include "G4SDManager.hh"
+#include "G4ios.hh"
+
+std::ofstream testdata("EVENTDATA.txt");
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4BeamTestSiSD::G4BeamTestSiSD(const G4String& name,
+	   const G4String& hitsCollectionName) 
+  : G4VSensitiveDetector(name)
+      , fHitsCollection(NULL)
+{
+   collectionName.insert(hitsCollectionName);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4BeamTestSiSD::~G4BeamTestSiSD() 
+{ // testdata  << "Arrival time"   << "  "  << "Energy "  <<  "  "  << "Distance"  <<  std::endl;
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4BeamTestSiSD::Initialize(G4HCofThisEvent* hce)
+{
+  // Create hits collection
+
+  testdata  <<  "*"    <<  std::endl;
+
+  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 G4BeamTestSiSD::ProcessHits(G4Step* aStep, 
+				G4TouchableHistory*)
+{  
+  //  G4double stepl = 10;
+  // aStep->SetStepLength(stepl);
+
+  G4String name=aStep->GetTrack()->GetDefinition()->GetParticleName();
+
+  /* G4cout << " Particle_name = " <<  name  << G4endl;   */
+
+  if (name == "opticalphoton" || name == "gamma") {
+/*  G4cout << " Particle_name = " <<  name  << G4endl; */
+
+// energy deposit
+  G4double edep = aStep->GetTotalEnergyDeposit();
+
+// if (edep==0.) return false;
+
+  /* G4cout << " Particle_name_after_edep = " <<  name  << G4endl;   */
+
+  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 zfix =  (zdelta + (zbef- 4.7))*1000;
+
+  //  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 << "TIME == " << time << G4endl; */
+ 
+ G4double vel =  aStep->GetPreStepPoint()->GetVelocity();
+ 
+ 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; 
+
+
+*/
+  return true;
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4BeamTestSiSD::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/G4BeamTestTank.cxx b/src/G4BeamTestTank.cxx
index 8a5193d..56570a4 100644
--- a/src/G4BeamTestTank.cxx
+++ b/src/G4BeamTestTank.cxx
@@ -14,7 +14,8 @@
 /* #include <boost/foreach.hpp> */
 
 #include "G4BeamTestTank.h"
-#include "G4TankIceSD.h"
+#include "G4BeamTestSiSD.h"
+/* #include "G4TankIceSD.h" */
 
 //prevent gcc to make something stupid with pretended unused variables
 #ifdef __GNUC__
@@ -55,7 +56,7 @@ G4VPhysicalVolume* G4BeamTestTank::InstallTank(G4VPhysicalVolume* mother, const
   // See also corresponding UserSpecialCuts in Physicslist !!!!
   // TODO(shivesh): Maybe do all of this as stepping action ??????
   G4UserLimits* energyLimit = new G4UserLimits();
-  energyLimit->SetUserMinEkine(264.1 * CLHEP::keV);  // Cherenkov threshold of electrons in water
+  energyLimit->SetUserMinEkine(1.907 * CLHEP::eV);  // Lower threshold of PMT - 600nm
 
   // std::string tankName=boost::lexical_cast<std::string>(tankKey_);
   std::string tankName = "BTT";
@@ -139,7 +140,7 @@ G4VPhysicalVolume* G4BeamTestTank::InstallTank(G4VPhysicalVolume* mother, const
       new G4LogicalVolume(upperglasssphere, glass,
                           ("log_dom_up_" + omName).c_str(), 0, 0, 0);
   G4LogicalVolume* logLowerGlass =
-      new G4LogicalVolume(lowerglasssphere, glass,
+      new G4LogicalVolume(lowerglasssphere, water,
                           ("log_dom_lo_" + omName).c_str(), 0, 0, 0);
   G4LogicalVolume* logUpperDOM =
       new G4LogicalVolume(upperdomsphere, effectiveDOM,
@@ -167,11 +168,12 @@ G4VPhysicalVolume* G4BeamTestTank::InstallTank(G4VPhysicalVolume* mother, const
   logLowerDOM->SetUserLimits(energyLimit);
   // }
 
-  // Define sensitive detector TODO(shivesh): make the PMT the SD
+  // Define sensitive detector
   G4SDManager* sdManager = G4SDManager::GetSDMpointer();
-  iceSD_ = new G4TankIceSD(("ice_SD_" + tankName).c_str());
+  iceSD_ = new G4BeamTestSiSD(("ice_SD_" + tankName).c_str(), "HitsCollection");
   sdManager->AddNewDetector(iceSD_);
-  logWater->SetSensitiveDetector(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));
@@ -196,29 +198,29 @@ G4VPhysicalVolume* G4BeamTestTank::InstallTank(G4VPhysicalVolume* mother, const
 }
 
 
-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::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) */
 /* { */
diff --git a/src/G4BeamTestUserStackingAction.cxx b/src/G4BeamTestUserStackingAction.cxx
new file mode 100644
index 0000000..73a354c
--- /dev/null
+++ b/src/G4BeamTestUserStackingAction.cxx
@@ -0,0 +1,76 @@
+#include "G4BeamTestUserStackingAction.h"
+
+#include "G4VProcess.hh"
+
+#include "G4ParticleDefinition.hh"
+#include "G4ParticleTypes.hh"
+#include "G4Track.hh"
+#include "G4ios.hh"
+
+#include "G4Event.hh"
+
+#include "G4HCofThisEvent.hh"
+#include "G4VHitsCollection.hh"
+#include "G4SDManager.hh"
+#include "G4ios.hh"
+#include "G4VDigitizerModule.hh"
+#include "G4DigiManager.hh"
+#include "G4Event.hh"
+#include "G4RunManager.hh"
+//#include "TH1F.h"
+//#include "TH2F.h"
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4BeamTestUserStackingAction::G4BeamTestUserStackingAction()
+  : G4UserStackingAction(),
+    fScintillationCounter(0), 
+    fCerenkovCounter(0)
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4BeamTestUserStackingAction::~G4BeamTestUserStackingAction()
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4ClassificationOfNewTrack
+G4BeamTestUserStackingAction::ClassifyNewTrack(const G4Track * aTrack)
+{
+  if(aTrack->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition())
+  { // particle is optical photon
+    if(aTrack->GetParentID()>0)
+    { // particle is secondary
+      if(aTrack->GetCreatorProcess()->GetProcessName() == "Scintillation")
+       fScintillationCounter++;
+      if(aTrack->GetCreatorProcess()->GetProcessName() == "Cerenkov")
+        fCerenkovCounter++;
+    }
+  }
+  return fUrgent;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4BeamTestUserStackingAction::NewStage()
+{
+  G4cout << "Number of Scintillation photons produced in this event : "
+        << fScintillationCounter << G4endl;
+  G4cout << "Number of Cerenkov photons produced in this event : "
+         << fCerenkovCounter << G4endl;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4BeamTestUserStackingAction::PrepareNewEvent()
+{
+ fScintillationCounter = 0;
+  fCerenkovCounter = 0;
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
diff --git a/src/G4BeamTestUserSteppingAction.cxx b/src/G4BeamTestUserSteppingAction.cxx
index e064628..f81ab62 100644
--- a/src/G4BeamTestUserSteppingAction.cxx
+++ b/src/G4BeamTestUserSteppingAction.cxx
@@ -5,11 +5,26 @@
 #include "G4LogicalVolume.hh"
 #include "G4UserLimits.hh"
 #include "G4SteppingManager.hh"
+#include "G4OpticalPhoton.hh"
+#include "G4RunManager.hh"
 
 G4BeamTestUserSteppingAction::G4BeamTestUserSteppingAction(){}
 
-void G4BeamTestUserSteppingAction::UserSteppingAction(const G4Step*)
+void G4BeamTestUserSteppingAction::UserSteppingAction(const G4Step* step)
 {
+  G4int eventNumber = G4RunManager::GetRunManager()->
+                                              GetCurrentEvent()->GetEventID();
+
+  if (eventNumber != fEventNumber) {
+     /* G4cout << " Number of Scintillation Photons in previous event: " */
+     /*       << fScintillationCounter << G4endl; */
+     /* G4cout << " Number of Cerenkov Photons in previous event: " */
+     /*        << fCerenkovCounter << G4endl; */
+     fEventNumber = eventNumber;
+     fScintillationCounter = 0;
+     fCerenkovCounter = 0;
+  }
+
   G4Track* track = fpSteppingManager->GetTrack();
   if(track)
   {
@@ -18,11 +33,33 @@ void G4BeamTestUserSteppingAction::UserSteppingAction(const G4Step*)
     if(!limit) G4cout << "----> G4LogicalVolume: " << volume->GetName() << " has no defined G4UserLimit" << G4endl;
     G4double threshold = limit->GetUserMinEkine(*track);
     //check if particle is a gamma
-    if(track->GetDefinition()->GetParticleName() == "gamma")
+    G4String particle = track->GetDefinition()->GetParticleName();
+    /* G4cout << "particle = " << particle << G4endl; */
+    if(particle == "gamma")
     {
       //check if particle energy is below threshold; if true, kill the particle
       G4double energy = track->GetTotalEnergy();
-      if(energy < threshold) track->SetTrackStatus(fStopAndKill);
+      if(energy < threshold){
+          G4cout << "SteppingAction: killing particle " << particle << " with energy " << energy << " < " << threshold << G4endl;
+          track->SetTrackStatus(fStopAndKill);
+      }
     }
   }
+
+  const std::vector<const G4Track*>* secondaries =
+                                            step->GetSecondaryInCurrentStep();
+
+  if (secondaries->size()>0) {
+     for(unsigned int i=0; i<secondaries->size(); ++i) {
+        if (secondaries->at(i)->GetParentID()>0) {
+           if(secondaries->at(i)->GetDynamicParticle()->GetParticleDefinition()
+               == G4OpticalPhoton::OpticalPhotonDefinition()){
+             if (secondaries->at(i)->GetCreatorProcess()->GetProcessName()
+              == "Scintillation")fScintillationCounter++;
+              if (secondaries->at(i)->GetCreatorProcess()->GetProcessName()
+               == "Cerenkov")fCerenkovCounter++;
+           }
+        }
+     }
+  }
 }
diff --git a/src/G4BeamTestUserTrackingAction.cxx b/src/G4BeamTestUserTrackingAction.cxx
index c6317df..f7c053b 100644
--- a/src/G4BeamTestUserTrackingAction.cxx
+++ b/src/G4BeamTestUserTrackingAction.cxx
@@ -24,11 +24,15 @@ void G4BeamTestUserTrackingAction::PostUserTrackingAction(const G4Track* track)
       for(size_t i=0;i<nSeco;i++)
       {
         //check if secondary particle is a gamma
-        if((*secondaries)[i]->GetDefinition()->GetParticleName() == "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) (*secondaries)[i]->SetTrackStatus(fStopAndKill);
+          if(energy < threshold){
+              G4cout << "TrackingAction: killing particle " << particle << " with energy " << energy << " < " << threshold << G4endl;
+              (*secondaries)[i]->SetTrackStatus(fStopAndKill);
+          }
         }
       }
     }
diff --git a/src/G4Interface.cxx b/src/G4Interface.cxx
index 75d6f78..e552173 100644
--- a/src/G4Interface.cxx
+++ b/src/G4Interface.cxx
@@ -5,7 +5,9 @@
 #include "G4BeamTestPhysicsList.h"
 #include "G4BeamTestUserTrackingAction.h"
 #include "G4BeamTestUserSteppingAction.h"
+#include "G4BeamTestUserStackingAction.h"
 #include "G4BeamTestPrimaryGeneratorAction.h"
+#include "G4BeamTestEventAction.h"
 #include "G4BeamTestRunAction.h"
 
 /* #include <icetray/I3Logging.h> */
@@ -165,6 +167,12 @@ void G4Interface::Initialize()
   G4cout << "Init UserSteppingAction ..." << G4endl;
   runManager_.SetUserAction(new G4BeamTestUserSteppingAction());
 
+  G4cout << "Init UserStackingAction ..." << G4endl;
+  runManager_.SetUserAction(new G4BeamTestUserStackingAction());
+
+  G4cout << "Init G4BeamTestEventAction ..." << G4endl;
+  runManager_.SetUserAction(new G4BeamTestEventAction());
+
   // Initialize G4 kernel
   /* log_debug("Init run manager ..."); */
   G4cout << "Init run manager ..." << G4endl;