feat: physics list

include/GeneAnalysisManager.hh

@@ -5,7 +5,7 @@
 #include "TROOT.h"
 #include "TTree.h"
 
-#include <globals.hh>
+#include "globals.hh"
 
 struct detEvent {
     G4double Thetacm;

include/GeneDetectorConstruction.hh

@@ -6,7 +6,7 @@ class G4VPhysicalVolume;
 class GeneDetectorConstructionMessenger;
 
 #include "G4VUserDetectorConstruction.hh"
-#include <globals.hh>
+#include "globals.hh"
 
 class GeneDetectorConstruction : public G4VUserDetectorConstruction {
 public:

include/GeneDetectorConstructionMessenger.hh

@@ -2,7 +2,7 @@
 #define GeneDetectorConstructionMessenger_hh 1
 
 #include "G4UImessenger.hh"
-#include <globals.hh>
+#include "globals.hh"
 
 class GeneDetectorConstruction;
 class G4UIdirectory;

include/GeneEventAction.hh

@@ -2,7 +2,7 @@
 #define GeneEventAction_hh 1
 
 #include "G4UserEventAction.hh"
-#include <globals.hh>
+#include "globals.hh"
 
 class GeneEventAction : public G4UserEventAction {
 public:

include/GeneHe3detSD.hh

@@ -3,7 +3,7 @@
 
 #include "G4VSensitiveDetector.hh"
 #include "GeneHe3detHit.hh"
-#include <globals.hh>
+#include "globals.hh"
 
 class G4Step;
 class G4HCofThisEvent;

include/GenePhysicsList.hh

@@ -1,7 +1,36 @@
-#ifndef GENEPHYSICSLIST_HH
+#ifndef GenePhysicsList_hh
-#define GENEPHYSICSLIST_HH 1
+#define GenePhysicsList_hh 1
 
 #include "G4VUserPhysicsList.hh"
 #include "globals.hh"
 
+class GenePhysicsList : public G4VUserPhysicsList {
+public:
+    GenePhysicsList();
+    virtual ~GenePhysicsList();
+    virtual void SetCuts();
+
+protected:
+    // construct particles and physics processes
+    virtual void ConstructParticle();
+    virtual void ConstructProcess();
+
+    virtual void ConstructGeneral();
+    virtual void ConstructEM();
+    virtual void ConstructOp();
+    virtual void ConstructHad();
+
+private:
+    G4int OpVerblevel;
+
+    G4double cutForGamma;
+    G4double cutForElectron;
+    G4double cutForPositron;
+
+    void ConstructMyBosons();
+    void ConstructMyLeptons();
+    void ConstructMyHadrons();
+    void ConstructMyShortLiveds();
+};
+
 #endif

include/GenePrimaryGeneratorAction.hh

@@ -2,7 +2,7 @@
 #define GenePrimaryGeneratorAction_hh 1
 
 #include "G4VUserPrimaryGeneratorAction.hh"
-#include <globals.hh>
+#include "globals.hh"
 
 class G4ParticleGun;
 class G4Event;

include/GenePrimaryGeneratorActionMessenger.hh

@@ -2,7 +2,7 @@
 #define GenePrimaryGeneratorActionMessenger_hh 1
 
 #include "G4UImessenger.hh"
-#include <globals.hh>
+#include "globals.hh"
 
 class GenePrimaryGeneratorAction;
 class G4UIdirectory;

include/GeneRunAction.hh

@@ -2,7 +2,7 @@
 #define GeneRunAction_hh 1
 
 #include "G4UserRunAction.hh"
-#include <globals.hh>
+#include "globals.hh"
 
 class G4String;
 class GeneRunActionMessenger;

include/GeneRunActionMessenger.hh

@@ -2,7 +2,7 @@
 #define GeneRunActionMessenger_hh 1
 
 #include "G4UImessenger.hh"
-#include <globals.hh>
+#include "globals.hh"
 
 class GeneRunAction;
 class G4UIdirectory;

main.cpp

@@ -11,27 +11,28 @@
 #include "G4VisExecutive.hh"
 #include "GeneDetectorConstruction.hh"
 #include "GeneEventAction.hh"
+#include "GenePhysicsList.hh"
 #include "GenePrimaryGeneratorAction.hh"
 #include "GeneRunAction.hh"
-#include "QGSP_BERT_HP.hh"
-#include "Randomize.hh"
 
 int main(int argc, char **argv) {
+    // random engine
+    // CLHEP::HepJamesRandom randomEngine;
+    G4long seeds;
+    seeds = time(NULL);
+    CLHEP::Ranlux64Engine randomEngine;
+    CLHEP::HepRandom::setTheEngine(&randomEngine);
+    CLHEP::HepRandom::setTheSeed(seeds);
+    CLHEP::HepRandom::showEngineStatus();
+
     G4UIExecutive *ui = nullptr;
     if (argc == 1) ui = new G4UIExecutive(argc, argv);
 
-    gROOT->GetInterpreter();
-    CLHEP::HepRandom::setTheEngine(new CLHEP::RanecuEngine);
-    CLHEP::HepRandom::setTheSeed(time(NULL));
-
     // run manager
     G4RunManager *runManager = new G4RunManager;
 
     runManager->SetUserInitialization(new GeneDetectorConstruction);
-
+    runManager->SetUserInitialization(new GenePhysicsList);
-    G4VModularPhysicsList *physicsList = new QGSP_BERT_HP;
-    physicsList->SetVerboseLevel(1);
-    runManager->SetUserInitialization(physicsList);
 
     runManager->SetUserAction(new GenePrimaryGeneratorAction);
     runManager->SetUserAction(new GeneRunAction);

src/GeneDetectorConstruction.cc

@@ -18,7 +18,7 @@
 #include "G4VisAttributes.hh"
 #include "GeneDetectorConstructionMessenger.hh"
 #include "GeneHe3detSD.hh"
-#include <globals.hh>
+#include "globals.hh"
 
 GeneDetectorConstruction::GeneDetectorConstruction() {
     Zoffset = 0;
@@ -106,7 +106,12 @@ G4VPhysicalVolume* GeneDetectorConstruction::Construct() {
     ej_200->AddElement(ele_H, 8.5 * perCent);
     ej_200->AddElement(ele_C, 91.5 * perCent);
 
-    // Material properties table 物质光学属性表
+    /*
+        Material properties table 物质光学属性表
+        Ref:
+            1. Geant4 Book For Application Developers, TrackingAndPhysics, physicsProcess, optical-photon-processes
+            2. Eljen-Catalog: Ej-200
+    */
     G4MaterialPropertiesTable* ej_200_MPT = new G4MaterialPropertiesTable();
     // Optical properties 光学属性
     const G4int nEntries = 21;
@@ -130,12 +135,12 @@ G4VPhysicalVolume* GeneDetectorConstruction::Construct() {
     // Intrinsic properties 本征属性
     ej_200_MPT->AddConstProperty("RESOLUTIONSCALE", 10.0);                // 本征分辨率
     ej_200_MPT->AddConstProperty("SCINTILLATIONYIELD", 10000. / MeV);     // 本征光产额
-    ej_200_MPT->AddConstProperty("YIELDRATIO", 1.0);                      // 快慢响应之比
+    ej_200_MPT->AddConstProperty("YIELDRATIO", 1.0);                      // 快成分占比
     ej_200_MPT->AddConstProperty("FASTSCINTILLATIONRISETIME", 0.9 * ns);  // 快成分上升时间
     ej_200_MPT->AddConstProperty("FASTTIMECONSTANT", 2.1 * ns);           // 快成分衰减时间
     ej_200->SetMaterialPropertiesTable(ej_200_MPT);
 
-    // ********** Boundary Optical Properties **********
+    // Boundary Optical Properties 边界光学属性
     G4OpticalSurface* ej_200_OpticalSurface = new G4OpticalSurface("ej_200_OpticalSurface");
     ej_200_OpticalSurface->SetType(dielectric_metal);
     ej_200_OpticalSurface->SetFinish(polished);

src/GenePhysicsList.cc

@@ -0,0 +1,746 @@
+#include "GenePhysicsList.hh"
+
+#include "G4ChipsKaonMinusInelasticXS.hh"
+#include "G4ChipsKaonPlusInelasticXS.hh"
+#include "G4ChipsKaonZeroInelasticXS.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4ParticleTable.hh"
+#include "G4ParticleTypes.hh"
+#include "G4ParticleWithCuts.hh"
+#include "G4ProcessManager.hh"
+#include "G4ProcessVector.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4UserLimits.hh"
+#include "G4ios.hh"
+#include "globals.hh"
+
+#include <iomanip>
+
+GenePhysicsList::GenePhysicsList() : G4VUserPhysicsList() {
+    defaultCutValue = 1.0 * micrometer;
+    cutForGamma = defaultCutValue;
+    cutForElectron = 1 * nanometer;
+    cutForPositron = defaultCutValue;
+
+    SetVerboseLevel(1);
+    OpVerblevel = 0;
+}
+
+GenePhysicsList::~GenePhysicsList() {}
+
+// ********** Construct Particles **********
+void GenePhysicsList::ConstructParticle() {
+    ConstructMyBosons();
+    ConstructMyLeptons();
+    ConstructMyHadrons();
+    ConstructMyShortLiveds();
+}
+
+// Boson 玻色子
+void GenePhysicsList::ConstructMyBosons() {
+    G4Geantino::GeantinoDefinition();
+    G4ChargedGeantino::ChargedGeantinoDefinition();
+    G4Gamma::GammaDefinition();
+    G4OpticalPhoton::OpticalPhotonDefinition();
+}
+
+// Lepoton 轻子
+void GenePhysicsList::ConstructMyLeptons() {
+    G4Electron::ElectronDefinition();
+    G4Positron::PositronDefinition();
+    G4MuonPlus::MuonPlusDefinition();
+    G4MuonMinus::MuonMinusDefinition();
+
+    G4NeutrinoE::NeutrinoEDefinition();
+    G4AntiNeutrinoE::AntiNeutrinoEDefinition();
+    G4NeutrinoMu::NeutrinoMuDefinition();
+    G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
+}
+
+#include "G4BaryonConstructor.hh"
+#include "G4IonConstructor.hh"
+#include "G4MesonConstructor.hh"
+
+// Hadron 强子
+void GenePhysicsList::ConstructMyHadrons() {
+    // Baryon 重子
+    G4BaryonConstructor bConstructor;
+    bConstructor.ConstructParticle();
+
+    // Meson 介子
+    G4MesonConstructor mConstructor;
+    mConstructor.ConstructParticle();
+
+    // Ion 离子
+    G4IonConstructor iConstructor;
+    iConstructor.ConstructParticle();
+}
+
+#include "G4ShortLivedConstructor.hh"
+
+// ShortLived 短寿命
+void GenePhysicsList::ConstructMyShortLiveds() {
+    G4ShortLivedConstructor slConstructor;
+    slConstructor.ConstructParticle();
+}
+
+// ********** Construct Processes **********
+void GenePhysicsList::ConstructProcess() {
+    AddTransportation();
+    ConstructGeneral();
+    ConstructEM();
+    ConstructHad();
+    ConstructOp();
+}
+
+// clang-format off
+// ******** Electromagnetic Processes ********
+// ****** gamma ******
+// 光电效应
+#include "G4PhotoElectricEffect.hh"
+#include "G4LivermorePhotoElectricModel.hh"
+// 康普顿散射
+#include "G4ComptonScattering.hh"
+#include "G4LivermoreComptonModel.hh"
+// 瑞利散射
+#include "G4RayleighScattering.hh"
+#include "G4LivermoreRayleighModel.hh"
+// 电子对
+#include "G4GammaConversion.hh"
+#include "G4LivermoreGammaConversionModel.hh"
+
+// ****** e- ******
+// 韧致辐射
+#include "G4eBremsstrahlung.hh"
+#include "G4LivermoreBremsstrahlungModel.hh"
+// 电离
+#include "G4eIonisation.hh"
+#include "G4LivermoreIonisationModel.hh"
+// 多重散射
+#include "G4eMultipleScattering.hh"
+
+// ****** e+ ******
+// 韧致辐射
+#include "G4eBremsstrahlung.hh"
+// 电离
+#include "G4eIonisation.hh"
+// 正负电子湮灭
+#include "G4eplusAnnihilation.hh"
+
+// ****** muno ******
+// 韧致辐射
+#include "G4MuBremsstrahlung.hh"
+// 电离
+#include "G4MuIonisation.hh"
+// 电子对
+#include "G4MuPairProduction.hh"
+// μ-捕获
+#include "G4MuonMinusCapture.hh"
+
+// ****** other ******
+// 离子参数化损失模型
+#include "G4IonParametrisedLossModel.hh"
+// 韧致辐射
+#include "G4hBremsstrahlung.hh"
+// 电离
+#include "G4hIonisation.hh"
+// 多重散射
+#include "G4hMultipleScattering.hh"
+// 离子电离
+#include "G4ionIonisation.hh"
+
+// em process options to allow msc step-limitation to be switched off
+#include "G4EmParameters.hh"
+#include "G4LossTableManager.hh"
+#include "G4UAtomicDeexcitation.hh"
+#include "G4VAtomDeexcitation.hh"
+// clang-format on
+
+void GenePhysicsList::ConstructEM() {
+    G4EmParameters* param = G4EmParameters::Instance();
+    // a finer grid of the physic tables in order to improve precision , 100 GeV with 200 bins
+    param->SetMaxEnergy(100 * GeV);
+    param->SetNumberOfBinsPerDecade(20);
+    param->SetMscStepLimitType(fMinimal);
+    // fluorescence 荧光
+    param->SetFluo(true);
+    // Particle (photons, electrons and ions) Induced X-ray Emission (PIXE)
+    param->SetPixe(true);
+    // Auger 俄歇电子
+    param->SetAuger(true);
+
+    G4LossTableManager* man = G4LossTableManager::Instance();
+    // Atomic Deexcitation 原子退激
+    G4VAtomDeexcitation* ad = man->AtomDeexcitation();
+    if (!ad) man->SetAtomDeexcitation(new G4UAtomicDeexcitation());
+
+    auto particleIterator = GetParticleIterator();
+    particleIterator->reset();
+    while ((*particleIterator)()) {
+        G4ParticleDefinition* particle = particleIterator->value();
+        G4ProcessManager* pmanager = particle->GetProcessManager();
+        G4String particleName = particle->GetParticleName();
+        G4String particleType = particle->GetParticleType();
+        G4double charge = particle->GetPDGCharge();
+
+        if (particleName == "gamma") {
+            G4PhotoElectricEffect* thePhotoElectricEffect = new G4PhotoElectricEffect();
+            thePhotoElectricEffect->SetEmModel(new G4LivermorePhotoElectricModel());
+            pmanager->AddDiscreteProcess(thePhotoElectricEffect);
+
+            G4ComptonScattering* theComptonScattering = new G4ComptonScattering();
+            theComptonScattering->SetEmModel(new G4LivermoreComptonModel());
+            pmanager->AddDiscreteProcess(theComptonScattering);
+
+            G4RayleighScattering* theRayleighScattering = new G4RayleighScattering();
+            theRayleighScattering->SetEmModel(new G4LivermoreRayleighModel());
+            pmanager->AddDiscreteProcess(theRayleighScattering);
+
+            G4GammaConversion* theGammaConversion = new G4GammaConversion();
+            theGammaConversion->SetEmModel(new G4LivermoreGammaConversionModel());
+            pmanager->AddDiscreteProcess(theGammaConversion);
+        } else if (particleName == "e-") {
+            G4eMultipleScattering* msc = new G4eMultipleScattering();
+            msc->SetStepLimitType(fUseDistanceToBoundary);
+            pmanager->AddProcess(msc,
+                                 -1,   // G4int ordAtRestDoIt
+                                 1,    // G4int ordAlongSteptDoIt
+                                 -1);  // G4int ordPostStepDoIt
+
+            // Ionisation
+            G4eIonisation* eIonisation = new G4eIonisation();
+            eIonisation->SetEmModel(new G4LivermoreIonisationModel());
+            eIonisation->SetStepFunction(0.2, 100 * um);  // improved precision in tracking
+            pmanager->AddProcess(eIonisation, -1, 2, 2);
+
+            // Bremsstrahlung
+            G4eBremsstrahlung* eBremsstrahlung = new G4eBremsstrahlung();
+            eBremsstrahlung->SetEmModel(new G4LivermoreBremsstrahlungModel());
+            pmanager->AddProcess(eBremsstrahlung, -1, -3, 3);
+        } else if (particleName == "e+") {
+            // positron
+            G4eMultipleScattering* msc = new G4eMultipleScattering();
+            msc->SetStepLimitType(fUseDistanceToBoundary);
+            pmanager->AddProcess(msc, -1, 1, 1);
+
+            // Ionisation
+            G4eIonisation* eIonisation = new G4eIonisation();
+            eIonisation->SetStepFunction(0.2, 100 * um);
+            pmanager->AddProcess(eIonisation, -1, 2, 2);
+
+            // Bremsstrahlung (use default, no low-energy available)
+            pmanager->AddProcess(new G4eBremsstrahlung(), -1, -1, 3);
+
+            // Annihilation
+            pmanager->AddProcess(new G4eplusAnnihilation(), 0, -1, 4);
+        } else if (particleName == "mu+" || particleName == "mu-") {
+            // muon
+            pmanager->AddProcess(new G4eMultipleScattering, -1, 1, 1);
+            pmanager->AddProcess(new G4MuIonisation(), -1, 2, 2);
+            pmanager->AddProcess(new G4MuBremsstrahlung(), -1, -1, 3);
+            pmanager->AddProcess(new G4MuPairProduction(), -1, -1, 4);
+            if (particleName == "mu-") pmanager->AddProcess(new G4MuonMinusCapture(), 0, -1, -1);
+        } else if (particleName == "proton" || particleName == "pi+" || particleName == "pi-") {
+            // multiple scattering
+            pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
+
+            // ionisation
+            G4hIonisation* hIonisation = new G4hIonisation();
+            hIonisation->SetStepFunction(0.2, 50 * um);
+            pmanager->AddProcess(hIonisation, -1, 2, 2);
+
+            // bremmstrahlung
+            pmanager->AddProcess(new G4hBremsstrahlung, -1, -3, 3);
+        } else if (particleName == "alpha" || particleName == "deuteron" || particleName == "triton" ||
+                   particleName == "He3") {
+            // multiple scattering
+            pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
+
+            // ionisation
+            G4ionIonisation* ionIoni = new G4ionIonisation();
+            ionIoni->SetStepFunction(0.1, 20 * um);
+            pmanager->AddProcess(ionIoni, -1, 2, 2);
+        } else if (particleName == "GenericIon") {
+            // OBJECT may be dynamically created as either a GenericIon or nucleus
+            // G4Nucleus exists and therefore has particle type nucleus
+            // genericIon:
+
+            // multiple scattering
+            pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
+
+            // ionisation
+            G4ionIonisation* ionIoni = new G4ionIonisation();
+            ionIoni->SetEmModel(new G4IonParametrisedLossModel());
+            ionIoni->SetStepFunction(0.1, 20 * um);
+            pmanager->AddProcess(ionIoni, -1, 2, 2);
+        }
+
+        else if ((!particle->IsShortLived()) && (charge != 0.0) && (particle->GetParticleName() != "chargedgeantino")) {
+            // all others charged particles except geantino
+            G4hMultipleScattering* aMultipleScattering = new G4hMultipleScattering();
+            G4hIonisation* ahadronIon = new G4hIonisation();
+
+            // multiple scattering
+            pmanager->AddProcess(aMultipleScattering, -1, 1, 1);
+
+            // ionisation
+            pmanager->AddProcess(ahadronIon, -1, 2, 2);
+        }
+    }
+}
+
+// ******** Optical Processes ********
+#include "G4OpAbsorption.hh"
+#include "G4OpBoundaryProcess.hh"
+#include "G4OpRayleigh.hh"
+#include "G4Scintillation.hh"
+
+void GenePhysicsList::ConstructOp() {
+    // default scintillation process
+    G4Scintillation* theScintProcessDef = new G4Scintillation("Scintillation");
+    theScintProcessDef->DumpPhysicsTable();
+    theScintProcessDef->SetTrackSecondariesFirst(true);
+    theScintProcessDef->SetScintillationYieldFactor(1.0);
+    theScintProcessDef->SetScintillationExcitationRatio(0.0);
+    theScintProcessDef->SetVerboseLevel(OpVerblevel);
+
+    // scintillation process for alpha:
+    G4Scintillation* theScintProcessAlpha = new G4Scintillation("Scintillation");
+    theScintProcessAlpha->SetTrackSecondariesFirst(true);
+    theScintProcessAlpha->SetScintillationYieldFactor(1.1);
+    theScintProcessAlpha->SetScintillationExcitationRatio(1.0);
+    theScintProcessAlpha->SetVerboseLevel(OpVerblevel);
+
+    // scintillation process for heavy nuclei
+    G4Scintillation* theScintProcessNuc = new G4Scintillation("Scintillation");
+    theScintProcessNuc->SetTrackSecondariesFirst(true);
+    theScintProcessNuc->SetScintillationYieldFactor(0.2);
+    theScintProcessNuc->SetScintillationExcitationRatio(1.0);
+    theScintProcessNuc->SetVerboseLevel(OpVerblevel);
+
+    // optical processes
+    G4OpAbsorption* theAbsorptionProcess = new G4OpAbsorption();
+    G4OpRayleigh* theRayleighScatteringProcess = new G4OpRayleigh();
+    G4OpBoundaryProcess* theBoundaryProcess = new G4OpBoundaryProcess();
+    theAbsorptionProcess->SetVerboseLevel(OpVerblevel);
+    theRayleighScatteringProcess->SetVerboseLevel(OpVerblevel);
+    theBoundaryProcess->SetVerboseLevel(OpVerblevel);
+
+    auto particleIterator = GetParticleIterator();
+    particleIterator->reset();
+    while ((*particleIterator)()) {
+        G4ParticleDefinition* particle = particleIterator->value();
+        G4ProcessManager* pmanager = particle->GetProcessManager();
+        G4String particleName = particle->GetParticleName();
+        if (theScintProcessDef->IsApplicable(*particle)) {
+            if (particle->GetParticleName() == "GenericIon") {
+                pmanager->AddProcess(theScintProcessNuc);
+                pmanager->SetProcessOrderingToLast(theScintProcessNuc, idxAtRest);
+                pmanager->SetProcessOrderingToLast(theScintProcessNuc, idxPostStep);
+            } else if (particle->GetParticleName() == "alpha") {
+                pmanager->AddProcess(theScintProcessAlpha);
+                pmanager->SetProcessOrderingToLast(theScintProcessAlpha, idxAtRest);
+                pmanager->SetProcessOrderingToLast(theScintProcessAlpha, idxPostStep);
+            } else {
+                pmanager->AddProcess(theScintProcessDef);
+                pmanager->SetProcessOrderingToLast(theScintProcessDef, idxAtRest);
+                pmanager->SetProcessOrderingToLast(theScintProcessDef, idxPostStep);
+            }
+        }
+
+        if (particleName == "opticalphoton") {
+            pmanager->AddDiscreteProcess(theAbsorptionProcess);
+            pmanager->AddDiscreteProcess(theRayleighScatteringProcess);
+            pmanager->AddDiscreteProcess(theBoundaryProcess);
+        }
+    }
+}
+
+// ******** Hadronic Processes ********
+// ****** Elastic ******
+#include "G4ChipsElasticModel.hh"
+#include "G4ElasticHadrNucleusHE.hh"
+#include "G4HadronElasticProcess.hh"
+
+// ****** Inelastic ******
+#include "G4AlphaInelasticProcess.hh"
+#include "G4AntiNeutronInelasticProcess.hh"
+#include "G4AntiProtonInelasticProcess.hh"
+#include "G4DeuteronInelasticProcess.hh"
+#include "G4KaonMinusInelasticProcess.hh"
+#include "G4KaonPlusInelasticProcess.hh"
+#include "G4KaonZeroLInelasticProcess.hh"
+#include "G4KaonZeroSInelasticProcess.hh"
+#include "G4NeutronInelasticProcess.hh"
+#include "G4PionMinusInelasticProcess.hh"
+#include "G4PionPlusInelasticProcess.hh"
+#include "G4ProtonInelasticProcess.hh"
+#include "G4TritonInelasticProcess.hh"
+
+// ****** High energy FTFP model and Bertini cascade ******
+#include "G4CascadeInterface.hh"
+#include "G4ExcitedStringDecay.hh"
+#include "G4FTFModel.hh"
+#include "G4GeneratorPrecompoundInterface.hh"
+#include "G4LundStringFragmentation.hh"
+#include "G4PreCompoundModel.hh"
+#include "G4TheoFSGenerator.hh"
+
+// ****** Cross sections ******
+#include "G4AntiNuclElastic.hh"
+#include "G4BGGNucleonInelasticXS.hh"
+#include "G4BGGPionElasticXS.hh"
+#include "G4ComponentAntiNuclNuclearXS.hh"
+#include "G4ComponentGGNuclNuclXsc.hh"
+#include "G4CrossSectionDataSetRegistry.hh"
+#include "G4CrossSectionElastic.hh"
+#include "G4CrossSectionInelastic.hh"
+#include "G4CrossSectionPairGG.hh"
+#include "G4HadronCaptureProcess.hh"
+#include "G4HadronElastic.hh"
+#include "G4PiNuclearCrossSection.hh"
+#include "G4VCrossSectionDataSet.hh"
+
+// ****** Neutron high-precision models: <20 MeV ******
+#include "G4ParticleHPCapture.hh"
+#include "G4ParticleHPCaptureData.hh"
+#include "G4ParticleHPElastic.hh"
+#include "G4ParticleHPElasticData.hh"
+#include "G4ParticleHPInelastic.hh"
+#include "G4ParticleHPInelasticData.hh"
+
+// ****** Stopping ******
+#include "G4AntiProtonAbsorptionFritiof.hh"
+#include "G4KaonMinusAbsorptionBertini.hh"
+#include "G4PiMinusAbsorptionBertini.hh"
+
+void GenePhysicsList::ConstructHad() {
+    // Elastic models
+    const G4double elastic_elimitPi = 1.0 * GeV;
+
+    G4HadronElastic* elastic_lhep0 = new G4HadronElastic();
+    G4HadronElastic* elastic_lhep1 = new G4HadronElastic();
+    elastic_lhep1->SetMaxEnergy(elastic_elimitPi);
+    G4ChipsElasticModel* elastic_chip = new G4ChipsElasticModel();
+    G4ElasticHadrNucleusHE* elastic_he = new G4ElasticHadrNucleusHE();
+    elastic_he->SetMinEnergy(elastic_elimitPi);
+
+    // Inelastic scattering
+    const G4double theFTFMin0 = 0.0 * GeV;
+    const G4double theFTFMin1 = 4.0 * GeV;
+    const G4double theFTFMax = 100.0 * TeV;
+    const G4double theBERTMin0 = 0.0 * GeV;
+    const G4double theBERTMin1 = 19.0 * MeV;
+    const G4double theBERTMax = 5.0 * GeV;
+    const G4double theHPMin = 0.0 * GeV;
+    const G4double theHPMax = 20.0 * MeV;
+
+    G4FTFModel* theStringModel = new G4FTFModel;
+    G4ExcitedStringDecay* theStringDecay = new G4ExcitedStringDecay(new G4LundStringFragmentation);
+    theStringModel->SetFragmentationModel(theStringDecay);
+    G4PreCompoundModel* thePreEquilib = new G4PreCompoundModel(new G4ExcitationHandler);
+    G4GeneratorPrecompoundInterface* theCascade = new G4GeneratorPrecompoundInterface(thePreEquilib);
+
+    G4TheoFSGenerator* theFTFModel0 = new G4TheoFSGenerator("FTFP");
+    theFTFModel0->SetHighEnergyGenerator(theStringModel);
+    theFTFModel0->SetTransport(theCascade);
+    theFTFModel0->SetMinEnergy(theFTFMin0);
+    theFTFModel0->SetMaxEnergy(theFTFMax);
+
+    G4TheoFSGenerator* theFTFModel1 = new G4TheoFSGenerator("FTFP");
+    theFTFModel1->SetHighEnergyGenerator(theStringModel);
+    theFTFModel1->SetTransport(theCascade);
+    theFTFModel1->SetMinEnergy(theFTFMin1);
+    theFTFModel1->SetMaxEnergy(theFTFMax);
+
+    G4CascadeInterface* theBERTModel0 = new G4CascadeInterface;
+    theBERTModel0->SetMinEnergy(theBERTMin0);
+    theBERTModel0->SetMaxEnergy(theBERTMax);
+
+    G4CascadeInterface* theBERTModel1 = new G4CascadeInterface;
+    theBERTModel1->SetMinEnergy(theBERTMin1);
+    theBERTModel1->SetMaxEnergy(theBERTMax);
+
+    G4VCrossSectionDataSet* thePiData = new G4CrossSectionPairGG(new G4PiNuclearCrossSection, 91 * GeV);
+    G4VCrossSectionDataSet* theAntiNucleonData = new G4CrossSectionInelastic(new G4ComponentAntiNuclNuclearXS);
+    G4ComponentGGNuclNuclXsc* ggNuclNuclXsec = new G4ComponentGGNuclNuclXsc();
+    G4VCrossSectionDataSet* theGGNuclNuclData = new G4CrossSectionInelastic(ggNuclNuclXsec);
+
+    auto particleIterator = GetParticleIterator();
+    particleIterator->reset();
+    while ((*particleIterator)()) {
+        G4ParticleDefinition* particle = particleIterator->value();
+        G4ProcessManager* pmanager = particle->GetProcessManager();
+        G4String particleName = particle->GetParticleName();
+
+        if (particleName == "pi+") {
+            // Elastic scattering
+            G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
+            theElasticProcess->AddDataSet(new G4BGGPionElasticXS(particle));
+            theElasticProcess->RegisterMe(elastic_lhep1);
+            theElasticProcess->RegisterMe(elastic_he);
+            pmanager->AddDiscreteProcess(theElasticProcess);
+            // Inelastic scattering
+            G4PionPlusInelasticProcess* theInelasticProcess = new G4PionPlusInelasticProcess("inelastic");
+            theInelasticProcess->AddDataSet(thePiData);
+            theInelasticProcess->RegisterMe(theFTFModel1);
+            theInelasticProcess->RegisterMe(theBERTModel0);
+            pmanager->AddDiscreteProcess(theInelasticProcess);
+        }
+
+        else if (particleName == "pi-") {
+            // Elastic scattering
+            G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
+            theElasticProcess->AddDataSet(new G4BGGPionElasticXS(particle));
+            theElasticProcess->RegisterMe(elastic_lhep1);
+            theElasticProcess->RegisterMe(elastic_he);
+            pmanager->AddDiscreteProcess(theElasticProcess);
+            // Inelastic scattering
+            G4PionMinusInelasticProcess* theInelasticProcess = new G4PionMinusInelasticProcess("inelastic");
+            theInelasticProcess->AddDataSet(thePiData);
+            theInelasticProcess->RegisterMe(theFTFModel1);
+            theInelasticProcess->RegisterMe(theBERTModel0);
+            pmanager->AddDiscreteProcess(theInelasticProcess);
+            // Absorption
+            pmanager->AddRestProcess(new G4PiMinusAbsorptionBertini, ordDefault);
+        }
+
+        else if (particleName == "kaon+") {
+            // Elastic scattering
+            G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
+            theElasticProcess->RegisterMe(elastic_lhep0);
+            pmanager->AddDiscreteProcess(theElasticProcess);
+            // Inelastic scattering
+            G4KaonPlusInelasticProcess* theInelasticProcess = new G4KaonPlusInelasticProcess("inelastic");
+            theInelasticProcess->AddDataSet(G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(
+                G4ChipsKaonPlusInelasticXS::Default_Name()));
+            theInelasticProcess->RegisterMe(theFTFModel1);
+            theInelasticProcess->RegisterMe(theBERTModel0);
+            pmanager->AddDiscreteProcess(theInelasticProcess);
+        }
+
+        else if (particleName == "kaon0S") {
+            // Elastic scattering
+            G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
+            theElasticProcess->RegisterMe(elastic_lhep0);
+            pmanager->AddDiscreteProcess(theElasticProcess);
+            // Inelastic scattering
+            G4KaonZeroSInelasticProcess* theInelasticProcess = new G4KaonZeroSInelasticProcess("inelastic");
+            theInelasticProcess->AddDataSet(G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(
+                G4ChipsKaonZeroInelasticXS::Default_Name()));
+            theInelasticProcess->RegisterMe(theFTFModel1);
+            theInelasticProcess->RegisterMe(theBERTModel0);
+            pmanager->AddDiscreteProcess(theInelasticProcess);
+        }
+
+        else if (particleName == "kaon0L") {
+            // Elastic scattering
+            G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
+            theElasticProcess->RegisterMe(elastic_lhep0);
+            pmanager->AddDiscreteProcess(theElasticProcess);
+            // Inelastic scattering
+            G4KaonZeroLInelasticProcess* theInelasticProcess = new G4KaonZeroLInelasticProcess("inelastic");
+            theInelasticProcess->AddDataSet(G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(
+                G4ChipsKaonZeroInelasticXS::Default_Name()));
+            theInelasticProcess->RegisterMe(theFTFModel1);
+            theInelasticProcess->RegisterMe(theBERTModel0);
+            pmanager->AddDiscreteProcess(theInelasticProcess);
+        }
+
+        else if (particleName == "kaon-") {
+            // Elastic scattering
+            G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
+            theElasticProcess->RegisterMe(elastic_lhep0);
+            pmanager->AddDiscreteProcess(theElasticProcess);
+            // Inelastic scattering
+            G4KaonMinusInelasticProcess* theInelasticProcess = new G4KaonMinusInelasticProcess("inelastic");
+            theInelasticProcess->AddDataSet(G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(
+                G4ChipsKaonMinusInelasticXS::Default_Name()));
+            theInelasticProcess->RegisterMe(theFTFModel1);
+            theInelasticProcess->RegisterMe(theBERTModel0);
+            pmanager->AddDiscreteProcess(theInelasticProcess);
+            pmanager->AddRestProcess(new G4KaonMinusAbsorptionBertini, ordDefault);
+        }
+
+        else if (particleName == "proton") {
+            // Elastic scattering
+            G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
+            theElasticProcess->AddDataSet(G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(
+                G4ChipsProtonElasticXS::Default_Name()));
+            theElasticProcess->RegisterMe(elastic_chip);
+            pmanager->AddDiscreteProcess(theElasticProcess);
+            // Inelastic scattering
+            G4ProtonInelasticProcess* theInelasticProcess = new G4ProtonInelasticProcess("inelastic");
+            theInelasticProcess->AddDataSet(new G4BGGNucleonInelasticXS(G4Proton::Proton()));
+            theInelasticProcess->RegisterMe(theFTFModel1);
+            theInelasticProcess->RegisterMe(theBERTModel0);
+            pmanager->AddDiscreteProcess(theInelasticProcess);
+        } else if (particleName == "anti_proton") {
+            // Elastic scattering
+            const G4double elastic_elimitAntiNuc = 100.0 * MeV;
+            G4AntiNuclElastic* elastic_anuc = new G4AntiNuclElastic();
+            elastic_anuc->SetMinEnergy(elastic_elimitAntiNuc);
+            G4CrossSectionElastic* elastic_anucxs = new G4CrossSectionElastic(elastic_anuc->GetComponentCrossSection());
+            G4HadronElastic* elastic_lhep2 = new G4HadronElastic();
+            elastic_lhep2->SetMaxEnergy(elastic_elimitAntiNuc);
+            G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
+            theElasticProcess->AddDataSet(elastic_anucxs);
+            theElasticProcess->RegisterMe(elastic_lhep2);
+            theElasticProcess->RegisterMe(elastic_anuc);
+            pmanager->AddDiscreteProcess(theElasticProcess);
+            // Inelastic scattering
+            G4AntiProtonInelasticProcess* theInelasticProcess = new G4AntiProtonInelasticProcess("inelastic");
+            theInelasticProcess->AddDataSet(theAntiNucleonData);
+            theInelasticProcess->RegisterMe(theFTFModel0);
+            pmanager->AddDiscreteProcess(theInelasticProcess);
+            // Absorption
+            pmanager->AddRestProcess(new G4AntiProtonAbsorptionFritiof, ordDefault);
+        }
+
+        else if (particleName == "neutron") {
+            // elastic scattering
+            G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
+            theElasticProcess->AddDataSet(G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(
+                G4ChipsNeutronElasticXS::Default_Name()));
+            G4HadronElastic* elastic_neutronChipsModel = new G4ChipsElasticModel();
+            elastic_neutronChipsModel->SetMinEnergy(19.0 * MeV);
+            theElasticProcess->RegisterMe(elastic_neutronChipsModel);
+            G4ParticleHPElastic* theElasticNeutronHP = new G4ParticleHPElastic;
+            theElasticNeutronHP->SetMinEnergy(theHPMin);
+            theElasticNeutronHP->SetMaxEnergy(theHPMax);
+            theElasticProcess->RegisterMe(theElasticNeutronHP);
+            theElasticProcess->AddDataSet(new G4ParticleHPElasticData);
+            pmanager->AddDiscreteProcess(theElasticProcess);
+            // inelastic scattering
+            G4NeutronInelasticProcess* theInelasticProcess = new G4NeutronInelasticProcess("inelastic");
+            theInelasticProcess->AddDataSet(new G4BGGNucleonInelasticXS(G4Neutron::Neutron()));
+            theInelasticProcess->RegisterMe(theFTFModel1);
+            theInelasticProcess->RegisterMe(theBERTModel1);
+            G4ParticleHPInelastic* theNeutronInelasticHPModel = new G4ParticleHPInelastic;
+            theNeutronInelasticHPModel->SetMinEnergy(theHPMin);
+            theNeutronInelasticHPModel->SetMaxEnergy(theHPMax);
+            theInelasticProcess->RegisterMe(theNeutronInelasticHPModel);
+            theInelasticProcess->AddDataSet(new G4ParticleHPInelasticData);
+            pmanager->AddDiscreteProcess(theInelasticProcess);
+            // capture
+            G4HadronCaptureProcess* theCaptureProcess = new G4HadronCaptureProcess;
+            G4ParticleHPCapture* theLENeutronCaptureModel = new G4ParticleHPCapture;
+            theLENeutronCaptureModel->SetMinEnergy(theHPMin);
+            theLENeutronCaptureModel->SetMaxEnergy(theHPMax);
+            theCaptureProcess->RegisterMe(theLENeutronCaptureModel);
+            theCaptureProcess->AddDataSet(new G4ParticleHPCaptureData);
+            pmanager->AddDiscreteProcess(theCaptureProcess);
+
+        } else if (particleName == "anti_neutron") {
+            // Elastic scattering
+            G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
+            theElasticProcess->RegisterMe(elastic_lhep0);
+            pmanager->AddDiscreteProcess(theElasticProcess);
+            // Inelastic scattering (include annihilation on-fly)
+            G4AntiNeutronInelasticProcess* theInelasticProcess = new G4AntiNeutronInelasticProcess("inelastic");
+            theInelasticProcess->AddDataSet(theAntiNucleonData);
+            theInelasticProcess->RegisterMe(theFTFModel0);
+            pmanager->AddDiscreteProcess(theInelasticProcess);
+        }
+
+        else if (particleName == "deuteron") {
+            // Elastic scattering
+            G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
+            theElasticProcess->RegisterMe(elastic_lhep0);
+            pmanager->AddDiscreteProcess(theElasticProcess);
+            // Inelastic scattering
+            G4DeuteronInelasticProcess* theInelasticProcess = new G4DeuteronInelasticProcess("inelastic");
+            theInelasticProcess->AddDataSet(theGGNuclNuclData);
+            theInelasticProcess->RegisterMe(theFTFModel1);
+            theInelasticProcess->RegisterMe(theBERTModel0);
+            pmanager->AddDiscreteProcess(theInelasticProcess);
+        }
+
+        else if (particleName == "triton") {
+            // Elastic scattering
+            G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
+            theElasticProcess->RegisterMe(elastic_lhep0);
+            pmanager->AddDiscreteProcess(theElasticProcess);
+            // Inelastic scattering
+            G4TritonInelasticProcess* theInelasticProcess = new G4TritonInelasticProcess("inelastic");
+            theInelasticProcess->AddDataSet(theGGNuclNuclData);
+            theInelasticProcess->RegisterMe(theFTFModel1);
+            theInelasticProcess->RegisterMe(theBERTModel0);
+            pmanager->AddDiscreteProcess(theInelasticProcess);
+        } else if (particleName == "alpha") {
+            // Elastic scattering
+            G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
+            theElasticProcess->RegisterMe(elastic_lhep0);
+            pmanager->AddDiscreteProcess(theElasticProcess);
+            // Inelastic scattering
+            G4AlphaInelasticProcess* theInelasticProcess = new G4AlphaInelasticProcess("inelastic");
+            theInelasticProcess->AddDataSet(theGGNuclNuclData);
+            theInelasticProcess->RegisterMe(theFTFModel1);
+            theInelasticProcess->RegisterMe(theBERTModel0);
+            pmanager->AddDiscreteProcess(theInelasticProcess);
+        }
+    }
+}
+
+// ******** Decays ********
+#include "G4Decay.hh"
+#include "G4IonTable.hh"
+#include "G4Ions.hh"
+#include "G4RadioactiveDecay.hh"
+
+void GenePhysicsList::ConstructGeneral() {
+    // Add Decay Process
+    G4Decay* theDecayProcess = new G4Decay();
+    auto particleIterator = GetParticleIterator();
+    particleIterator->reset();
+    while ((*particleIterator)()) {
+        G4ParticleDefinition* particle = particleIterator->value();
+        G4ProcessManager* pmanager = particle->GetProcessManager();
+
+        if (theDecayProcess->IsApplicable(*particle) && !particle->IsShortLived()) {
+            pmanager->AddProcess(theDecayProcess);
+            // set ordering for PostStepDoIt and AtRestDoIt
+            pmanager->SetProcessOrdering(theDecayProcess, idxPostStep);
+            pmanager->SetProcessOrdering(theDecayProcess, idxAtRest);
+        }
+    }
+
+    // Declare radioactive decay to the GenericIon in the IonTable.
+    const G4IonTable* theIonTable = G4ParticleTable::GetParticleTable()->GetIonTable();
+    G4RadioactiveDecay* theRadioactiveDecay = new G4RadioactiveDecay();
+
+    for (G4int i = 0; i < theIonTable->Entries(); i++) {
+        G4String particleName = theIonTable->GetParticle(i)->GetParticleName();
+        G4String particleType = theIonTable->GetParticle(i)->GetParticleType();
+
+        if (particleName == "GenericIon") {
+            G4ProcessManager* pmanager = theIonTable->GetParticle(i)->GetProcessManager();
+            pmanager->SetVerboseLevel(verboseLevel);
+            pmanager->AddProcess(theRadioactiveDecay);
+            pmanager->SetProcessOrdering(theRadioactiveDecay, idxPostStep);
+            pmanager->SetProcessOrdering(theRadioactiveDecay, idxAtRest);
+        }
+    }
+}
+
+// ******** Cuts ********
+void GenePhysicsList::SetCuts() {
+    if (verboseLevel > 1) G4cout << "GenePhysicsList::SetCuts:";
+
+    if (verboseLevel > 0) {
+        G4cout << "GenePhysicsList::SetCuts:";
+        G4cout << "CutLength : " << G4BestUnit(defaultCutValue, "Length") << G4endl;
+    }
+
+    // special for low energy physics
+    G4double lowlimit = 250 * eV;
+    G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange(lowlimit, 100. * GeV);
+
+    // set cut values for gamma at first and for e- second and next for e+,
+    // because some processes for e+/e- need cut values for gamma
+    SetCutValue(cutForGamma, "gamma");
+    SetCutValue(cutForElectron, "e-");
+    SetCutValue(cutForPositron, "e+");
+
+    if (verboseLevel > 0) DumpCutValuesTable();
+}

src/GenePrimaryGeneratorAction.cc

@@ -11,7 +11,7 @@
 #include "GenePrimaryGeneratorAction.hh"
 #include "GenePrimaryGeneratorActionMessenger.hh"
 #include "Randomize.hh"
-#include <globals.hh>
+#include "globals.hh"
 
 #include "fstream"