Tue 20 Nov 17:26:02 GMT 2018

13 files changed, 656 insertions, 111 deletions

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;