Merge pull request #12 from rest-for-physics/revise_generators

Revise generators

Author: jgalan

.gitlab-ci.yml

@@ -1,4 +1,4 @@
-image: lobis/root-geant4-garfieldpp:cxx14_ROOTv6-25-01_Geant4v10.4.3
+image: nkx1231/root6-geant4-garfield:0.8
 
 stages:
   - build
@@ -16,7 +16,7 @@ before_script:
     - apt update && apt install -y wget
 
 # This is the one we use to run validations tests. We save the installation as an artifact.
-build_root_6.22_Geant4_10.4.2:
+build:
   type: build
   script:
     - echo "**${CI_PROJECT_DIR}**"
@@ -26,7 +26,9 @@ build_root_6.22_Geant4_10.4.2:
     - git submodule init source/libraries/geant4
     - git submodule update source/libraries/geant4
     - cd source/libraries/geant4/
-    - git checkout master
+      # We attempt to pull a branch from geant4lib named as the one at restG4. 
+      # If not it will pull master by default
+    - ../../../scripts/checkoutRemoteBranch.sh ${CI_COMMIT_BRANCH}
     - cd ../../../
     - git submodule init source/packages/restG4
     - git submodule update source/packages/restG4

include/PrimaryGeneratorAction.hh

@@ -2,16 +2,19 @@
 #ifndef PrimaryGeneratorAction_h
 #define PrimaryGeneratorAction_h 1
 
+#include <TF3.h>
 #include <TH1D.h>
 
+#include <fstream>
+#include <iostream>
+
 #include <G4IonTable.hh>
 #include <G4ParticleGun.hh>
 #include <G4VUserPrimaryGeneratorAction.hh>
-#include <fstream>
 #include <globals.hh>
-#include <iostream>
 
 #include "DetectorConstruction.hh"
+#include "TRestGeant4Particle.h"
 
 using namespace std;
 
@@ -28,54 +31,21 @@ class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction {
     virtual void GeneratePrimaries(G4Event*);
     G4ParticleGun* GetParticleGun() { return fParticleGun; };
 
-    void SetSpectrum(TH1D* spt, double eMin = 0, double eMax = 0) {
-        TString xLabel = (TString)spt->GetXaxis()->GetTitle();
-
-        if (xLabel.Contains("MeV")) {
-            energyFactor = 1.e3;
-        } else if (xLabel.Contains("GeV")) {
-            energyFactor = 1.e6;
-        } else {
-            energyFactor = 1.;
-        }
-
-        fSpectrum = spt;
-        fSpectrumIntegral = fSpectrum->Integral();
-
-        startEnergyBin = 1;
-        endEnergyBin = fSpectrum->GetNbinsX();
-
-        if (eMin > 0) {
-            for (int i = startEnergyBin; i <= endEnergyBin; i++) {
-                if (fSpectrum->GetBinCenter(i) > eMin) {
-                    startEnergyBin = i;
-                    break;
-                }
-            }
-        }
-
-        if (eMax > 0) {
-            for (int i = startEnergyBin; i <= endEnergyBin; i++) {
-                if (fSpectrum->GetBinCenter(i) > eMax) {
-                    endEnergyBin = i;
-                    break;
-                }
-            }
-        }
-
-        fSpectrumIntegral = fSpectrum->Integral(startEnergyBin, endEnergyBin);
-    }
+    void SetSpectrum(TH1D* spt, double eMin = 0, double eMax = 0);
+    void SetGeneratorSpatialDensity(TString str);
 
     void SetAngularDistribution(TH1D* ang) { fAngularDistribution = ang; }
 
    private:
+    vector<TRestGeant4Particle> fTempParticles;
+
     G4ParticleGun* fParticleGun;
     DetectorConstruction* fDetector;
-
     G4ParticleDefinition* fParticle = nullptr;
 
     TH1D* fSpectrum;
     TH1D* fAngularDistribution;
+    TF3* fGeneratorSpatialDensityFunction;
 
     Int_t startEnergyBin;
     Int_t endEnergyBin;
@@ -88,15 +58,26 @@ class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction {
     Double_t lastEnergy;
 
     void SetParticlePosition();
-    void SetParticlePosition(int n);
-    G4ParticleDefinition* SetParticleDefinition(int n);
-    void SetParticleEnergy(int n);
-    void SetParticleDirection(int n);
+    G4ParticleDefinition* SetParticleDefinition(Int_t particlesourceindex, TRestGeant4Particle p);
+    void SetParticleEnergy(Int_t particlesourceindex, TRestGeant4Particle p);
+    void SetParticleDirection(Int_t particlesourceindex, TRestGeant4Particle p);
 
     G4ThreeVector GetIsotropicVector();
     Double_t GetAngle(G4ThreeVector x, G4ThreeVector y);
     Double_t GetCosineLowRandomThetaAngle();
 
+    void GenPositionOnGDMLVolume(double& x, double& y, double& z);
+    void GenPositionOnGDMLSurface(double& x, double& y, double& z);
+    void GenPositionOnBoxVolume(double& x, double& y, double& z);
+    void GenPositionOnBoxSurface(double& x, double& y, double& z);
+    void GenPositionOnSphereVolume(double& x, double& y, double& z);
+    void GenPositionOnSphereSurface(double& x, double& y, double& z);
+    void GenPositionOnCylinderVolume(double& x, double& y, double& z);
+    void GenPositionOnCylinderSurface(double& x, double& y, double& z);
+    void GenPositionOnPoint(double& x, double& y, double& z);
+    void GenPositionOnWall(double& x, double& y, double& z);
+    void GenPositionOnPlate(double& x, double& y, double& z);
+
     G4String fParType;
     G4String fGenType;
     G4double fParEnergy;

restG4.cc

@@ -182,12 +182,12 @@ int main(int argc, char** argv) {
     // set user action classes
     PrimaryGeneratorAction* prim = new PrimaryGeneratorAction(det);
 
-    if (restG4Metadata->GetParticleSource(0).GetEnergyDistType() == "TH1D") {
-        TString fileFullPath = (TString)restG4Metadata->GetParticleSource(0).GetSpectrumFilename();
+    if (restG4Metadata->GetParticleSource(0)->GetEnergyDistType() == "TH1D") {
+        TString fileFullPath = (TString)restG4Metadata->GetParticleSource(0)->GetSpectrumFilename();
 
         TFile fin(fileFullPath);
 
-        TString sptName = restG4Metadata->GetParticleSource(0).GetSpectrumName();
+        TString sptName = restG4Metadata->GetParticleSource(0)->GetSpectrumName();
 
         TH1D* h = (TH1D*)fin.Get(sptName);
         ;
@@ -200,22 +200,22 @@ int main(int argc, char** argv) {
 
         initialEnergySpectrum = *h;
 
-        Double_t minEnergy = restG4Metadata->GetParticleSource(0).GetMinEnergy();
+        Double_t minEnergy = restG4Metadata->GetParticleSource(0)->GetMinEnergy();
         if (minEnergy < 0) minEnergy = 0;
 
-        Double_t maxEnergy = restG4Metadata->GetParticleSource(0).GetMaxEnergy();
+        Double_t maxEnergy = restG4Metadata->GetParticleSource(0)->GetMaxEnergy();
         if (maxEnergy < 0) maxEnergy = 0;
 
         // We set the initial spectrum energy provided from TH1D
         prim->SetSpectrum(&initialEnergySpectrum, minEnergy, maxEnergy);
     }
 
-    if (restG4Metadata->GetParticleSource(0).GetAngularDistType() == "TH1D") {
-        TString fileFullPath = (TString)restG4Metadata->GetParticleSource(0).GetAngularFilename();
+    if (restG4Metadata->GetParticleSource(0)->GetAngularDistType() == "TH1D") {
+        TString fileFullPath = (TString)restG4Metadata->GetParticleSource(0)->GetAngularFilename();
 
         TFile fin(fileFullPath);
 
-        TString sptName = restG4Metadata->GetParticleSource(0).GetAngularName();
+        TString sptName = restG4Metadata->GetParticleSource(0)->GetAngularName();
         TH1D* h = (TH1D*)fin.Get(sptName);
         ;
         if (h == NULL) {
@@ -311,13 +311,13 @@ int main(int argc, char** argv) {
         while (biasing) {
             // Definning isotropic gammas using the spectrum and angular distribution
             // previously obtained
-            restG4Metadata->RemoveSources();
+            restG4Metadata->RemoveParticleSources();
 
-            TRestGeant4ParticleSource src;
-            src.SetParticleName("gamma");
-            src.SetEnergyDistType("TH1D");
-            src.SetAngularDistType("TH1D");
-            restG4Metadata->AddSource(src);
+            TRestGeant4ParticleSource* src = new TRestGeant4ParticleSource();
+            src->SetParticleName("gamma");
+            src->SetEnergyDistType("TH1D");
+            src->SetAngularDistType("TH1D");
+            restG4Metadata->AddParticleSource(src);
 
             // We set the spectrum from previous biasing volume inside the primary
             // generator
@@ -331,8 +331,8 @@ int main(int argc, char** argv) {
             // the biasing volume
             restG4Metadata->SetGeneratorType(
                 restG4Metadata->GetBiasingVolume(biasing).GetBiasingVolumeType());
-            restG4Metadata->SetGeneratorSize(
-                restG4Metadata->GetBiasingVolume(biasing).GetBiasingVolumeSize());
+            double size = restG4Metadata->GetBiasingVolume(biasing).GetBiasingVolumeSize();
+            restG4Metadata->SetGeneratorSize(TVector3(size,size,size));
             // restG4Metadata->GetBiasingVolume( biasing-1 ).PrintBiasingVolume();
 
             // Definning biasing the number of event to be re-launched