Add Cerenkov Processes

16 files changed, 503 insertions, 48 deletions

diff --git a/include/G4BeamTestEMPhysics.h b/include/G4BeamTestEMPhysics.h
index bf20b62..e7a88dd 100644
--- a/include/G4BeamTestEMPhysics.h
+++ b/include/G4BeamTestEMPhysics.h
@@ -11,6 +11,17 @@
 #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.
@@ -30,6 +41,7 @@ public:
 
   void ConstructParticle();
   void ConstructProcess();
+  void ConstructOp();
 
 private:
   // Gamma physics
@@ -47,6 +59,9 @@ private:
   G4eIonisation positronIonisation; 
   G4eBremsstrahlung positronBremsStrahlung;  
   G4eplusAnnihilation annihilation;
+
+  // Cerenkov physics
+  G4Cerenkov cerenkov;
 };
 
 #endif // G4TANKRESPONSE_G4BEAMTESTEMPHYSICS_H_INCLUDED
diff --git a/include/G4BeamTestPrimaryGeneratorAction.h b/include/G4BeamTestPrimaryGeneratorAction.h
new file mode 100644
index 0000000..83fbacc
--- /dev/null
+++ b/include/G4BeamTestPrimaryGeneratorAction.h
@@ -0,0 +1,33 @@
+#ifndef G4BeamTestPrimaryGeneratorAction_h
+#define G4BeamTestPrimaryGeneratorAction_h 1
+
+#include "G4VUserPrimaryGeneratorAction.hh"
+#include "globals.hh"
+
+class G4GeneralParticleSource;
+class G4Event;
+class G4BeamTestPrimaryGeneratorMessenger;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+class G4BeamTestPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
+{
+  public:
+    G4BeamTestPrimaryGeneratorAction();
+    virtual ~G4BeamTestPrimaryGeneratorAction();
+
+  public:
+    virtual void GeneratePrimaries(G4Event*);
+
+    void SetOptPhotonPolar();
+    void SetOptPhotonPolar(G4double);
+
+  private:
+    G4GeneralParticleSource* fParticleGun;
+    G4BeamTestPrimaryGeneratorMessenger* fGunMessenger;
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif /*G4BeamTestPrimaryGeneratorAction_h*/
+
diff --git a/include/G4BeamTestPrimaryGeneratorMessenger.h b/include/G4BeamTestPrimaryGeneratorMessenger.h
new file mode 100644
index 0000000..0063f95
--- /dev/null
+++ b/include/G4BeamTestPrimaryGeneratorMessenger.h
@@ -0,0 +1,30 @@
+#ifndef G4BeamTestPrimaryGeneratorMessenger_h
+#define G4BeamTestPrimaryGeneratorMessenger_h 1
+
+#include "G4UImessenger.hh"
+#include "globals.hh"
+
+class G4BeamTestPrimaryGeneratorAction;
+class G4UIdirectory;
+class G4UIcmdWithADoubleAndUnit;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+class G4BeamTestPrimaryGeneratorMessenger: public G4UImessenger
+{
+  public:
+    G4BeamTestPrimaryGeneratorMessenger(G4BeamTestPrimaryGeneratorAction* );
+    virtual ~G4BeamTestPrimaryGeneratorMessenger();
+ 
+    virtual void SetNewValue(G4UIcommand*, G4String);
+ 
+  private:
+    G4BeamTestPrimaryGeneratorAction* fG4BeamTestAction;
+    G4UIdirectory*                  fGunDir;
+    G4UIcmdWithADoubleAndUnit*      fPolarCmd;
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif
+
diff --git a/include/G4BeamTestRunAction.h b/include/G4BeamTestRunAction.h
new file mode 100644
index 0000000..b25ad7e
--- /dev/null
+++ b/include/G4BeamTestRunAction.h
@@ -0,0 +1,30 @@
+#ifndef G4BeamTestRunAction_h
+#define G4BeamTestRunAction_h 1
+
+#include "globals.hh"
+#include "G4UserRunAction.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+class G4Timer;
+class G4Run;
+
+class G4BeamTestRunAction : public G4UserRunAction
+{
+  public:
+    G4BeamTestRunAction();
+    virtual ~G4BeamTestRunAction();
+
+  public:
+    virtual G4Run* GenerateRun();
+    virtual void BeginOfRunAction(const G4Run* aRun);
+    virtual void EndOfRunAction(const G4Run* aRun);
+
+  private:
+    G4Timer* fTimer;
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif /*G4BeamTestRunAction_h*/
+
diff --git a/include/G4BeamTestRunManager.h b/include/G4BeamTestRunManager.h.backup
index a7eaff0..a7eaff0 100644
--- a/include/G4BeamTestRunManager.h
+++ b/include/G4BeamTestRunManager.h.backup
diff --git a/include/G4BeamTestTank.h b/include/G4BeamTestTank.h
index c7a1713..3f2fd8c 100644
--- a/include/G4BeamTestTank.h
+++ b/include/G4BeamTestTank.h
@@ -80,7 +80,6 @@ class G4BeamTestTank
   G4double innerRadius_;
   G4double outerRadius_;
 
-  G4double waterHeight_;
   G4double airHeight_;
 
   G4double glassOuterRadius_;
diff --git a/include/G4Interface.h b/include/G4Interface.h
index 4b88f00..e2956af 100644
--- a/include/G4Interface.h
+++ b/include/G4Interface.h
@@ -1,7 +1,8 @@
 #ifndef _TOPSIM_G4INTERFACE_H_
 #define _TOPSIM_G4INTERFACE_H_
 
-#include "G4BeamTestRunManager.h"
+#include "G4RunManager.hh"
+/* #include "G4BeamTestRunManager.h" */
 /* #include <icetray/I3Logging.h> */
 
 #ifdef G4VIS_USE
@@ -33,17 +34,17 @@ class G4Interface
   /// To be called after simulating each IceTray event.
   void TerminateEvent();
   /// Simulate a single particle (InitializeEvent must be called first)
-  void InjectParticle(const std::string& particleName,
-          const G4ThreeVector& particlePosition,
-          const G4ThreeVector& particleDirection,
-          const G4double particleEnergy);
+  /* void InjectParticle(const std::string& particleName, */
+  /*         const G4ThreeVector& particlePosition, */
+  /*         const G4ThreeVector& particleDirection, */
+  /*         const G4double particleEnergy); */
 
  private:
   void Initialize();
 
   static G4Interface* g4Interface_;
 
-  G4BeamTestRunManager runManager_;
+  G4RunManager runManager_;
 
 #ifdef G4VIS_USE
   G4VisManager* visManager_;
diff --git a/run1.mac b/run1.mac
index fceb1c4..7f9c02a 100644
--- a/run1.mac
+++ b/run1.mac
@@ -8,4 +8,5 @@
 /gps/ang/type iso
 
 /gps/direction 1 0 0
+/gps/position -2 0 0 m
 /run/beamOn 1
diff --git a/src/G4BeamTestDetectorConstruction.cxx b/src/G4BeamTestDetectorConstruction.cxx
index 2bf76a2..95a2071 100644
--- a/src/G4BeamTestDetectorConstruction.cxx
+++ b/src/G4BeamTestDetectorConstruction.cxx
@@ -34,14 +34,14 @@ G4VPhysicalVolume* G4BeamTestDetectorConstruction::Construct()
   /* origin_.set(delaunay.GetOrigin().x(), delaunay.GetOrigin().y(), zSnowBottom + zHalfLength); */
 
   // Determine World dimensions
-  G4double xWorld = 20.0 * CLHEP::m;
-  G4double yWorld = 20.0 * CLHEP::m;
-  G4double zWorld = 20.0 * CLHEP::m;
+  G4double xWorld = 2.0 * CLHEP::m;
+  G4double yWorld = 2.0 * CLHEP::m;
+  G4double zWorld = 2.0 * CLHEP::m;
 
   // Create world volume
   G4Box* world_box = new G4Box("solid_world", xWorld, yWorld, zWorld);
   G4LogicalVolume* worldLog =
-      new G4LogicalVolume(world_box, G4Material::GetMaterial("G4_AIR"), "log_world", 0, 0, 0);
+      new G4LogicalVolume(world_box, G4Material::GetMaterial("Air"), "log_world", 0, 0, 0);
   G4VPhysicalVolume* worldPhys =
       new G4PVPlacement(0, G4ThreeVector(), worldLog, "world", 0, false, 0);
 
diff --git a/src/G4BeamTestEMPhysics.cxx b/src/G4BeamTestEMPhysics.cxx
index 55ecf99..1aa9a7c 100644
--- a/src/G4BeamTestEMPhysics.cxx
+++ b/src/G4BeamTestEMPhysics.cxx
@@ -55,4 +55,66 @@ 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/G4BeamTestPrimaryGeneratorAction.cxx b/src/G4BeamTestPrimaryGeneratorAction.cxx
new file mode 100644
index 0000000..192ab12
--- /dev/null
+++ b/src/G4BeamTestPrimaryGeneratorAction.cxx
@@ -0,0 +1,84 @@
+#include "G4BeamTestPrimaryGeneratorAction.h"
+#include "G4BeamTestPrimaryGeneratorMessenger.h"
+
+#include "Randomize.hh"
+
+#include "G4Event.hh"
+#include "G4GeneralParticleSource.hh"
+#include "G4ParticleTable.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4SystemOfUnits.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4BeamTestPrimaryGeneratorAction::G4BeamTestPrimaryGeneratorAction()
+ : G4VUserPrimaryGeneratorAction(), 
+fParticleGun(0)
+{
+  G4int n_particle = 1;
+//  fParticleGun = new G4GeneralParticleSource(n_particle);
+  fParticleGun = new G4GeneralParticleSource();  
+
+  //create a messenger for this class
+  fGunMessenger = new G4BeamTestPrimaryGeneratorMessenger(this);
+
+  //default kinematic
+  //
+  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
+  G4ParticleDefinition* particle = particleTable->FindParticle("mu-");
+
+  fParticleGun->SetParticleDefinition(particle);
+  fParticleGun->SetParticleTime(0.0*ns);
+  fParticleGun->SetParticlePosition(G4ThreeVector(-76.2*cm,0.0*cm,-21.59*cm));
+  //fParticleGun->SetParticleMomentumDirection(G4ThreeVector(1.,0.,0.));
+  //fParticleGun->SetParticleEnergy(2.0*GeV);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4BeamTestPrimaryGeneratorAction::~G4BeamTestPrimaryGeneratorAction()
+{
+  delete fParticleGun;
+  delete fGunMessenger;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4BeamTestPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
+{
+  fParticleGun->GeneratePrimaryVertex(anEvent);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4BeamTestPrimaryGeneratorAction::SetOptPhotonPolar()
+{
+ G4double angle = G4UniformRand() * 360.0*deg;
+ SetOptPhotonPolar(angle);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4BeamTestPrimaryGeneratorAction::SetOptPhotonPolar(G4double angle)
+{
+ if (fParticleGun->GetParticleDefinition()->GetParticleName()!="opticalphoton")
+   {
+     G4cout << "--> warning from PrimaryGeneratorAction::SetOptPhotonPolar() :"
+               "the particleGun is not an opticalphoton" << G4endl;
+     return;
+   }
+
+ G4ThreeVector normal (1., 0., 0.);
+ G4ThreeVector kphoton = fParticleGun->GetParticleMomentumDirection();
+ G4ThreeVector product = normal.cross(kphoton);
+ G4double modul2       = product*product;
+ 
+ G4ThreeVector e_perpend (0., 0., 1.);
+ if (modul2 > 0.) e_perpend = (1./std::sqrt(modul2))*product;
+ G4ThreeVector e_paralle    = e_perpend.cross(kphoton);
+ 
+ G4ThreeVector polar = std::cos(angle)*e_paralle + std::sin(angle)*e_perpend;
+ fParticleGun->SetParticlePolarization(polar);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/src/G4BeamTestPrimaryGeneratorMessenger.cxx b/src/G4BeamTestPrimaryGeneratorMessenger.cxx
new file mode 100644
index 0000000..a7ff777
--- /dev/null
+++ b/src/G4BeamTestPrimaryGeneratorMessenger.cxx
@@ -0,0 +1,53 @@
+#include "G4BeamTestPrimaryGeneratorMessenger.h"
+#include "G4BeamTestPrimaryGeneratorAction.h"
+
+#include "G4UIdirectory.hh"
+#include "G4UIcmdWithADoubleAndUnit.hh"
+#include "G4SystemOfUnits.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4BeamTestPrimaryGeneratorMessenger::
+  G4BeamTestPrimaryGeneratorMessenger(G4BeamTestPrimaryGeneratorAction* G4BeamTestGun)
+  : G4UImessenger(),
+    fG4BeamTestAction(G4BeamTestGun)
+{
+  fGunDir = new G4UIdirectory("/G4BeamTest/gun/");
+  fGunDir->SetGuidance("PrimaryGenerator control");
+
+  fPolarCmd =
+           new G4UIcmdWithADoubleAndUnit("/G4BeamTest/gun/optPhotonPolar",this);
+  fPolarCmd->SetGuidance("Set linear polarization");
+  fPolarCmd->SetGuidance("  angle w.r.t. (k,n) plane");
+  fPolarCmd->SetParameterName("angle",true);
+  fPolarCmd->SetUnitCategory("Angle");
+  fPolarCmd->SetDefaultValue(-360.0);
+  fPolarCmd->SetDefaultUnit("deg");
+  fPolarCmd->AvailableForStates(G4State_Idle);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4BeamTestPrimaryGeneratorMessenger::~G4BeamTestPrimaryGeneratorMessenger()
+{
+  delete fPolarCmd;
+  delete fGunDir;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4BeamTestPrimaryGeneratorMessenger::SetNewValue(
+                                        G4UIcommand* command, G4String newValue)
+{
+  if( command == fPolarCmd ) {
+      G4double angle = fPolarCmd->GetNewDoubleValue(newValue);
+      if ( angle == -360.0*deg ) {
+         fG4BeamTestAction->SetOptPhotonPolar();
+      } else {
+         fG4BeamTestAction->SetOptPhotonPolar(angle);
+      }
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
diff --git a/src/G4BeamTestRunAction.cxx b/src/G4BeamTestRunAction.cxx
new file mode 100644
index 0000000..bc0f542
--- /dev/null
+++ b/src/G4BeamTestRunAction.cxx
@@ -0,0 +1,131 @@
+// Make this appear first!
+#include "G4Timer.hh"
+
+#include "G4BeamTestRunAction.h"
+
+#include "G4Run.hh"
+
+#include "G4RunManager.hh"
+#include "G4ScoringManager.hh"
+#include "G4UnitsTable.hh"
+#include "G4LogicalVolume.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4PSNofCollision.hh"
+#include "G4SDManager.hh"
+#include "G4VSensitiveDetector.hh"
+
+#include "G4Run.hh"
+
+//#include "TH1F.h"  //This file seems to be ....
+
+//extern TH1F*  h_bragg9;
+//#include "B3PrimaryGeneratorAction.hh"
+//#include "B3Run.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4BeamTestRunAction::G4BeamTestRunAction()
+ : G4UserRunAction(),
+   fTimer(0)
+{
+  fTimer = new G4Timer;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4BeamTestRunAction::~G4BeamTestRunAction()
+{
+  delete fTimer;
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4Run* G4BeamTestRunAction::GenerateRun()
+{
+  return new G4Run;
+}
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4BeamTestRunAction::BeginOfRunAction(const G4Run* aRun)
+{
+  G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
+  fTimer->Start();
+
+  G4RunManager::GetRunManager()->SetRandomNumberStore(false);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+
+//
+
+void G4BeamTestRunAction::EndOfRunAction(const G4Run* aRun)
+{
+  fTimer->Stop();
+  G4cout << "number of event = " << aRun->GetNumberOfEvent()
+         << " " << *fTimer << G4endl;
+
+/*
+  int binmax = h_bragg9->GetMaximumBin(); double x = h_bragg9->GetXaxis()->GetBinCenter(binmax);
+
+  G4cout << " bin max JFGDSLKJF ==== " << binmax << "X VALUE DISTANCE SAUJDFHALSDF ==== "<< x << G4endl
+*/
+
+}
+
+//
+/*
+G4SDManager* SDman = G4SDManager::GetSDMpointer();
+G4cout << " Print out hit information" << G4endl;
+  G4SDManager* SDManager= G4SDManager::GetSDMpointer();
+  G4String detName = "inside_dom_lo_wintery_mix";
+  G4MultiFunctionalDetector* mfd=
+    (G4MultiFunctionalDetector*)SDManager-> FindSensitiveDetector(detName);
+  // std::cout << "number of mfd " << muonSD->GetNumberOfPrimitives() << std::endl;
+  // G4VPrimitiveScorer* scorer = muonSD->GetPrimitive(0);
+  // std::cout << muonSD->GetPrimitive(0) << std::endl;
+  // scorer-> PrintAll();
+  // muonSD-> DrawAll();
+  // G4cout << G4endl;
+  // G4PSNofCollision* muonSD2=
+  //   (G4PSNofCollision*)SDManager-> FindSensitiveDetector("inside_dom_up_winterymix");
+  // muonSD2-> PrintAll();
+    if ( mfd ){
+    //--- Loop over the registered primitive scorers.
+        for (G4int icol = 0; icol < mfd->GetNumberOfPrimitives(); icol++){
+    // Get Primitive Scorer object.
+            G4VPrimitiveScorer* scorer=mfd->GetPrimitive(icol);
+    // collection name and collectionID for HitsCollection,
+    // where type of HitsCollection is G4THitsMap in case of primitive scorer.
+    // The collection name is given by <MFD name>/<Primitive Scorer name>.
+            G4String collectionName = scorer->GetName();
+            G4String fullCollectionName = detName+"/"+collectionName;
+            G4int    collectionID = SDman->GetCollectionID(fullCollectionName);
+
+            if ( collectionID >= 0 ){
+                G4cout << "++ "<<fullCollectionName<< " id " << collectionID 
+                       << G4endl;
+    // Store obtained HitsCollection information into data members.
+    // And, creates new G4THitsMap for accumulating quantities during RUN.
+                std::cout << "fullCollectionName " << fullCollectionName << std::endl;
+                std::cout <<"collectionID " << collectionID << std::endl;
+                G4THitsMap<G4double> * hits =
+                    new G4THitsMap<G4double>(detName ,collectionName);
+                hits->PrintAllHits();
+            }else{
+                G4cout << "** collection " << fullCollectionName 
+                       << " not found. "<<G4endl;
+            }
+        }
+    }
+*/
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+
+
diff --git a/src/G4BeamTestRunManager.cxx b/src/G4BeamTestRunManager.cxx.backup
index 5ed2e05..5ed2e05 100644
--- a/src/G4BeamTestRunManager.cxx
+++ b/src/G4BeamTestRunManager.cxx.backup
diff --git a/src/G4BeamTestTank.cxx b/src/G4BeamTestTank.cxx
index b9d469e..8a5193d 100644
--- a/src/G4BeamTestTank.cxx
+++ b/src/G4BeamTestTank.cxx
@@ -27,13 +27,14 @@
 G4BeamTestTank::G4BeamTestTank()
 {
   // Get tank dimensions
-  tankThickness_ = 0.5*CLHEP::cm; // TODO(shivesh) : check thickness
+  tankThickness_ = 0.0*CLHEP::cm; // TODO(shivesh) : check thickness
   tankHeight_  = 76.83 * 2.54 * CLHEP::cm;
   innerRadius_ = 32 * 2.54 * CLHEP::cm;
   outerRadius_ = innerRadius_ + tankThickness_;
 
   // Get fill height
   fillHeight_ = 55.1 * 2.54 * CLHEP::cm;
+  airHeight_ = tankHeight_ - fillHeight_;
 
   // Set DOM dimensions
   glassOuterRadius_ = 6.5 * 2.54 * CLHEP::cm; // 6.5" outer glass sphere radius
@@ -81,6 +82,8 @@ G4VPhysicalVolume* G4BeamTestTank::InstallTank(G4VPhysicalVolume* mother, const
 
   // 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);
@@ -88,6 +91,7 @@ G4VPhysicalVolume* G4BeamTestTank::InstallTank(G4VPhysicalVolume* mother, const
       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;
   G4VPhysicalVolume* physAir_UNUSED =
       new G4PVPlacement(0, physAirPosition, logAir,
                         ("air_" + tankName).c_str(), tankLog_, false, 0);
@@ -104,6 +108,8 @@ G4VPhysicalVolume* G4BeamTestTank::InstallTank(G4VPhysicalVolume* mother, const
 
   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;
 
diff --git a/src/G4Interface.cxx b/src/G4Interface.cxx
index 9f44d3e..75d6f78 100644
--- a/src/G4Interface.cxx
+++ b/src/G4Interface.cxx
@@ -5,6 +5,8 @@
 #include "G4BeamTestPhysicsList.h"
 #include "G4BeamTestUserTrackingAction.h"
 #include "G4BeamTestUserSteppingAction.h"
+#include "G4BeamTestPrimaryGeneratorAction.h"
+#include "G4BeamTestRunAction.h"
 
 /* #include <icetray/I3Logging.h> */
 
@@ -73,46 +75,47 @@ void G4Interface::InitializeEvent()
 
   if(!eventInitialized_)
   {
-    runManager_.InitializeRun();
+    /* runManager_.InitializeRun(); */
+    runManager_.Initialize();
     eventInitialized_ = true;
   }
 }
 
 
-void G4Interface::InjectParticle(
-        const std::string& particleName, const G4ThreeVector& particlePosition,
-        const G4ThreeVector& particleDirection, const G4double particleEnergy
-        )
-{
-  if(!eventInitialized_)
-  {
-    /* log_fatal("No event initialized. Cannot inject particle!"); */
-    G4cout << "No event initialized. Cannot inject particle!" << G4endl;
-    return;
-  }
-
-  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
-  G4ParticleDefinition* particleDef = NULL;
-  particleDef = particleTable->FindParticle(particleName);
-
-  // Transform coordinates to world system
-  G4ThreeVector position = particlePosition - detector_->GetWorldOrigin();
-
-  G4ParticleGun gun(1);
-  gun.SetParticleDefinition(particleDef);
-  gun.SetParticleEnergy(particleEnergy);
-  gun.SetParticlePosition(position);
-  gun.SetParticleMomentumDirection(particleDirection);
-
-  G4cout << "Injecting %s: x=%.2f m, y=%.2f m, z=%.2f m, E=%.3f MeV",
-            particleName.c_str(),
-            position.x() / CLHEP::m,
-            position.y() / CLHEP::m,
-            position.z() / CLHEP::m,
-            gun.GetParticleEnergy() / CLHEP::MeV;
-
-  runManager_.InjectParticle(&gun);
-}
+/* void G4Interface::InjectParticle( */
+/*         const std::string& particleName, const G4ThreeVector& particlePosition, */
+/*         const G4ThreeVector& particleDirection, const G4double particleEnergy */
+/*         ) */
+/* { */
+/*   if(!eventInitialized_) */
+/*   { */
+/*     #<{(| log_fatal("No event initialized. Cannot inject particle!"); |)}># */
+/*     G4cout << "No event initialized. Cannot inject particle!" << G4endl; */
+/*     return; */
+/*   } */
+/*  */
+/*   G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable(); */
+/*   G4ParticleDefinition* particleDef = NULL; */
+/*   particleDef = particleTable->FindParticle(particleName); */
+/*  */
+/*   // Transform coordinates to world system */
+/*   G4ThreeVector position = particlePosition - detector_->GetWorldOrigin(); */
+/*  */
+/*   G4ParticleGun gun(1); */
+/*   gun.SetParticleDefinition(particleDef); */
+/*   gun.SetParticleEnergy(particleEnergy); */
+/*   gun.SetParticlePosition(position); */
+/*   gun.SetParticleMomentumDirection(particleDirection); */
+/*  */
+/*   G4cout << "Injecting %s: x=%.2f m, y=%.2f m, z=%.2f m, E=%.3f MeV", */
+/*             particleName.c_str(), */
+/*             position.x() / CLHEP::m, */
+/*             position.y() / CLHEP::m, */
+/*             position.z() / CLHEP::m, */
+/*             gun.GetParticleEnergy() / CLHEP::MeV; */
+/*  */
+/*   runManager_.InjectParticle(&gun); */
+/* } */
 
 
 void G4Interface::TerminateEvent()
@@ -124,7 +127,8 @@ void G4Interface::TerminateEvent()
 
   if(eventInitialized_)
   {
-    runManager_.TerminateRun();
+    /* runManager_.TerminateRun(); */
+      runManager_.RunTermination();
     eventInitialized_ = false;
   }
 }
@@ -147,6 +151,12 @@ void G4Interface::Initialize()
   G4cout << "Init physics list ..." << G4endl;
   runManager_.SetUserInitialization(new G4BeamTestPhysicsList());
 
+  G4cout << "Init PrimaryGeneratorAction ..." << G4endl;
+  runManager_.SetUserAction(new G4BeamTestPrimaryGeneratorAction());
+
+  G4cout << "Init RunAction ..." << G4endl;
+  runManager_.SetUserAction(new G4BeamTestRunAction());
+
   /* log_debug("Init UserTrackingAction ..."); */
   G4cout << "Init UserTrackingAction ..." << G4endl;
   runManager_.SetUserAction(new G4BeamTestUserTrackingAction());