Add Cerenkov Processes

8 files changed, 387 insertions, 41 deletions

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());