PrimaryGeneratorAction: new isotropic generator

.gitignore

@@ -0,0 +1,3 @@
+*build*
+.idea
+.vscode

CMakeLists.txt

@@ -1,7 +1,7 @@
 project(restG4)
 
 CMAKE_MINIMUM_REQUIRED(VERSION 3.5)
-message("===============  restG4  ==============")
+message("==================================================  restG4  =================================================")
 set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/../../cmake ${CMAKE_MODULE_PATH})
 
 # Find REST
@@ -33,28 +33,20 @@ find_package(ROOT REQUIRED)
 message("-- Found ROOT version : ${ROOT_VERSION}")
 
 # Set include and libs
-set(external_libs ${Geant4_LIBRARIES} ${REST_LIBRARIES} ${ROOT_LIBRARIES} -L$ENV{ROOTSYS}/lib -lGui -lGeom -lGdml)
+set(external_libs ${Geant4_LIBRARIES} ${ROOT_LIBRARIES} ${REST_LIBRARIES})
 set(external_include_dirs ${PROJECT_SOURCE_DIR}/include ${Geant4_INCLUDE_DIR} ${ROOT_INCLUDE_DIRS} ${REST_INCLUDE_DIRS})
 
 include_directories(${incdir})
 
 file(GLOB sources ${PROJECT_SOURCE_DIR}/src/*.cc)
 file(GLOB headers ${PROJECT_SOURCE_DIR}/include/*.hh)
 
-#----------------------------------------------------------------------------
-# Add the executable, and link it to the Geant4 libraries
-#----------------------------------------------------------------------------
-#message( STATUS "external_libs G4: ${external_libs}" )
-#message( STATUS "external_include_dirs G4: ${external_include_dirs}" )
-#message( STATUS "sources G4: ${sources}" )
-#message( STATUS "headers G4: ${headers}" )
 add_executable(${PROJECT_NAME} restG4.cc ${sources} ${headers})
 target_link_libraries(${PROJECT_NAME} ${external_libs})
 target_include_directories(${PROJECT_NAME} SYSTEM PUBLIC ${external_include_dirs})
 
 target_compile_options(${PROJECT_NAME} PRIVATE)
 
-
 #----------------------------------------------------------------------------
 # Copy all scripts to the build directory, i.e. the directory in which we
 # build sphereSimLab. This is so that we can run the executable directly because it

include/PrimaryGeneratorAction.hh

@@ -35,21 +35,23 @@
 #ifndef PrimaryGeneratorAction_h
 #define PrimaryGeneratorAction_h 1
 
+#include <TH1D.h>
+
+#include <fstream>
+#include <iostream>
+
 #include "DetectorConstruction.hh"
 #include "G4IonTable.hh"
 #include "G4ParticleGun.hh"
 #include "G4VUserPrimaryGeneratorAction.hh"
 #include "globals.hh"
 
-#include <TH1D.h>
-
-#include <fstream>
-#include <iostream>
 using namespace std;
 
 const int nSpct = 3000;
 
 class G4Event;
+class G4SPSAngDistribution;
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
@@ -58,7 +60,6 @@ class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction {
     PrimaryGeneratorAction(DetectorConstruction* pDetector);
     ~PrimaryGeneratorAction();
 
-   public:
     virtual void GeneratePrimaries(G4Event*);
     G4ParticleGun* GetParticleGun() { return fParticleGun; };
 
@@ -106,6 +107,8 @@ class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction {
     G4ParticleGun* fParticleGun;
     DetectorConstruction* fDetector;
 
+    G4SPSAngDistribution* fG4AngDist;
+
     TH1D* fSpectrum;
     TH1D* fAngularDistribution;
 

src/PrimaryGeneratorAction.cc

@@ -5,14 +5,15 @@
 #include <TRestGeant4Event.h>
 #include <TRestGeant4Metadata.h>
 
-#include "G4Event.hh"
-#include "G4Geantino.hh"
-#include "G4IonTable.hh"
-#include "G4ParticleDefinition.hh"
-#include "G4ParticleTable.hh"
-#include "G4SystemOfUnits.hh"
-#include "G4UnitsTable.hh"
-#include "Randomize.hh"
+#include <G4Event.hh>
+#include <G4Geantino.hh>
+#include <G4IonTable.hh>
+#include <G4ParticleDefinition.hh>
+#include <G4ParticleTable.hh>
+#include <G4SPSAngDistribution.hh>
+#include <G4SystemOfUnits.hh>
+#include <G4UnitsTable.hh>
+#include <Randomize.hh>
 
 extern TRestGeant4Metadata* restG4Metadata;
 extern TRestGeant4Event* restG4Event;
@@ -24,18 +25,23 @@ Int_t face = 0;
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 PrimaryGeneratorAction::PrimaryGeneratorAction(DetectorConstruction* pDetector)
-    : G4VUserPrimaryGeneratorAction(), fParticleGun(0), fDetector(pDetector) {
+    : G4VUserPrimaryGeneratorAction(), fParticleGun(nullptr), fDetector(pDetector), fG4AngDist(nullptr) {
     G4int n_particle = 1;
     fParticleGun = new G4ParticleGun(n_particle);
 
     nCollections = restG4Metadata->GetPrimaryGenerator().GetNumberOfCollections();
 
     nBiasingVolumes = restG4Metadata->GetNumberOfBiasingVolumes();
+
+    fG4AngDist = new G4SPSAngDistribution();
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-PrimaryGeneratorAction::~PrimaryGeneratorAction() { delete fParticleGun; }
+PrimaryGeneratorAction::~PrimaryGeneratorAction() {
+    delete fParticleGun;
+    delete fG4AngDist;
+}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
@@ -665,20 +671,11 @@ void PrimaryGeneratorAction::SetParticlePosition(int n) {
 
 //_____________________________________________________________________________
 G4ThreeVector PrimaryGeneratorAction::GetIsotropicVector() {
-    G4double a, b, c;
-    G4double n;
-    do {
-        a = (G4UniformRand() - 0.5) / 0.5;
-        b = (G4UniformRand() - 0.5) / 0.5;
-        c = (G4UniformRand() - 0.5) / 0.5;
-        n = a * a + b * b + c * c;
-    } while (n > 1 || n == 0.0);
-
-    n = std::sqrt(n);
-    a /= n;
-    b /= n;
-    c /= n;
-    return G4ThreeVector(a, b, c);
+    G4String isotropicDistribution = "iso";
+    if (fG4AngDist->GetDistType() != isotropicDistribution) {
+        fG4AngDist->SetAngDistType(isotropicDistribution);
+    }
+    return fG4AngDist->GenerateOne();
 }
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 //