feat: add Ne22 target

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:
@@ -17,8 +17,6 @@ public:
     G4VPhysicalVolume* Construct();
     inline void SetZoffset(G4double val) { Zoffset = val; };
     inline void SetZactiveshift(G4double val) { Zactiveshift = val; };
-    inline void SetBfraction(G4double val) { PE_Bfraction = val; };
-    inline void SetDetReaction(G4String val) { SReaction = val; };
 
 private:
     G4LogicalVolume* logicalWorld;
@@ -26,8 +24,6 @@ private:
 
     G4double Zoffset;
     G4double Zactiveshift;
-    G4double PE_Bfraction;
-    G4String SReaction;
     GeneDetectorConstructionMessenger* pMessenger;
 };
 

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
-#define GENEPHYSICSLIST_HH 1
+#ifndef GenePhysicsList_hh
+#define GenePhysicsList_hh 1
 
 #include "G4VUserPhysicsList.hh"
 #include "globals.hh"
 
-#endif
+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;
@@ -17,7 +17,6 @@ public:
     void GeneratePrimaries(G4Event* anEvent);
     void SetBeamEnergy(G4double);
     void SetAngType(G4String&);
-    void SetReaction(G4String&);
 
 private:
     G4ParticleGun* particleGun;
@@ -31,7 +30,6 @@ private:
     G4double p4_cm;
 
     TH1F* hAng;
-    G4int iReaction;
     GenePrimaryGeneratorActionMessenger* pMessenger;
 };
 

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);
-
-    G4VModularPhysicsList *physicsList = new QGSP_BERT_HP;
-    physicsList->SetVerboseLevel(1);
-    runManager->SetUserInitialization(physicsList);
+    runManager->SetUserInitialization(new GenePhysicsList);
 
     runManager->SetUserAction(new GenePrimaryGeneratorAction);
     runManager->SetUserAction(new GeneRunAction);

run.mac

@@ -12,13 +12,9 @@
 # 2nd energy
 # 3rd angular
 
-# reaction setting
-/Gene/PrimaryGA/SetReaction C13an           # V51pn or C13an
-/Gene/Detector/SetDetReaction C13an         # V51pn or C13an, default: C13an
-
 # energy setting
-/Gene/PrimaryGA/SetBeamEnergy 1.0 MeV       # alpha energy for 13C(a,n)16O reaction
-											# proton energy for 51V(p,n)51Cr reaction
+/Gene/PrimaryGA/SetBeamEnergy 1.0 MeV       # alpha energy for 22Ne(a,n)25Mg reaction
+
 # angular setting
 /Gene/PrimaryGA/SetAngular ISO              # ISO or ENDF
 
@@ -30,7 +26,7 @@
 /Gene/Detector/SetBfraction 0.054           # perCent, default: 0.054%
 
 # user setting for runNo of output file
-/Gene/Run/RunNo 00_C13an_ISO_1000keV_zoffset0mm_zshift0mm_B054_100000evts_Coll60mm
+/Gene/Run/RunNo 00_Ne22an_ISO_1000keV_zoffset0mm_zshift0mm_B054_100000evts_Coll60mm
 
 /run/initialize
 

src/GeneAnalysisManager.cc

@@ -2,7 +2,7 @@
 
 GeneAnalysisManager* GeneAnalysisManager::instance = 0;
 
-GeneAnalysisManager::GeneAnalysisManager() : analysisFileName("Simulation_C13an.root"), theTFile(0), tr(0) {
+GeneAnalysisManager::GeneAnalysisManager() : analysisFileName("Simulation_Ne22an.root"), theTFile(0), tr(0) {
     // evtNo = 0;
     memset(&mydata, 0x00, sizeof(mydata));
 }
@@ -43,7 +43,7 @@ void GeneAnalysisManager::book() {
     theTFile = new TFile(analysisFileName, "RECREATE");
 
     // create a tree
-    tr = new TTree("SimData1", "C13an simulation data");
+    tr = new TTree("SimData1", "Ne22an simulation data");
     tr->Branch("Thetacm", &mydata.Thetacm, "Thetacm/D");
     tr->Branch("Thetalab", &mydata.Thetalab, "Thetalab/D");
     tr->Branch("Ndets", &mydata.Ndets, "Ndets/I");

src/GeneDetectorConstruction.cc

@@ -18,13 +18,11 @@
 #include "G4VisAttributes.hh"
 #include "GeneDetectorConstructionMessenger.hh"
 #include "GeneHe3detSD.hh"
-#include <globals.hh>
+#include "globals.hh"
 
 GeneDetectorConstruction::GeneDetectorConstruction() {
     Zoffset = 0;
     Zactiveshift = 0;
-    PE_Bfraction = 0.054;
-    SReaction = "C13an";
     pMessenger = new GeneDetectorConstructionMessenger(this);
 }
 
@@ -43,18 +41,8 @@ G4VPhysicalVolume* GeneDetectorConstruction::Construct() {
     G4Element* ele_C = manager->FindOrBuildElement("C");
     // Boron 硼
     G4Material* mat_B = manager->FindOrBuildMaterial("G4_B");
-    // Copper 铜
-    G4Material* mat_Cu = manager->FindOrBuildMaterial("G4_Cu");
-    // Gold 金
-    G4Material* mat_Au = manager->FindOrBuildMaterial("G4_Au");
-    // Aluminum 铝
-    G4Material* mat_Al = manager->FindOrBuildMaterial("G4_Al");
-    // Tantalum 钽
-    G4Material* mat_Ta = manager->FindOrBuildMaterial("G4_Ta");
     // Stainless stell 不锈钢
     G4Material* mat_stell = manager->FindOrBuildMaterial("G4_STAINLESS-STEEL");
-    // Water 水
-    G4Material* mat_Water = manager->FindOrBuildMaterial("G4_WATER");
     // Air 空气
     G4Material* mat_Air = manager->FindOrBuildMaterial("G4_AIR");
     // Polyethylene 聚乙烯
@@ -62,15 +50,10 @@ G4VPhysicalVolume* GeneDetectorConstruction::Construct() {
 
     // shielding 屏蔽层 PE(B) 含硼聚乙烯
     density = 0.938 * g / cm3;
-    G4Material* shield_PE = new G4Material("PE_B", density, ncomponents = 2);
+    G4Material* shield_PE = new G4Material("BPE", density, ncomponents = 2);
     shield_PE->AddMaterial(mat_PE, fractionmass = 93 * perCent);
     shield_PE->AddMaterial(mat_B, fractionmass = 7 * perCent);
 
-    // moderator 慢化剂 PE(B) 含硼聚乙烯
-    G4Material* moderator_PE = new G4Material("PE_B_impurity", density, ncomponents = 2);
-    moderator_PE->AddMaterial(mat_PE, fractionmass = (100.0 - PE_Bfraction) * perCent);
-    moderator_PE->AddMaterial(mat_B, fractionmass = PE_Bfraction * perCent);
-
     // Vacuum 真空
     temperature = 77 * kelvin;
     pressure = 1.0e-5 * pascal;
@@ -90,14 +73,15 @@ G4VPhysicalVolume* GeneDetectorConstruction::Construct() {
     G4Material* work_gas = new G4Material("Workgas", density, ncomponents = 1, kStateGas, temperature, pressure);
     work_gas->AddElement(ele_He, 100. * perCent);
 
-    // C13 Target C13靶
-    G4Isotope* iso_C12 = new G4Isotope(name = "C12", z = 6.0, n = 12.0, a = 12.0 * g / mole);
-    G4Isotope* iso_C13 = new G4Isotope(name = "C13", z = 6.0, n = 13.0, a = 13.003 * g / mole);
-    G4Element* ele_C13 = new G4Element(name = "Carbon", symbol = "C13", ncomponents = 2);
-    ele_C13->AddIsotope(iso_C12, 1.0 * perCent);
-    ele_C13->AddIsotope(iso_C13, 99.0 * perCent);
-    G4Material* target_C = new G4Material("targetC", density = 2.3 * g / cm3, ncomponents = 1);
-    target_C->AddElement(ele_C13, 100 * perCent);
+    // Ne22 Target Ne22靶
+    temperature = 300 * kelvin;
+    pressure = 100 * pascal;
+    density = (0.981146 * kg / m3) * (pressure / atmosphere) * (273.15 * kelvin / temperature);
+    G4Isotope* iso_Ne22 = new G4Isotope(name = "Ne22", z = 10, n = 12, a = 21.99138 * g / mole);
+    G4Element* ele_Ne22 = new G4Element(name = "Neon-22", symbol = "Ne22", ncomponents = 1);
+    ele_Ne22->AddIsotope(iso_Ne22, 100 * perCent);
+    G4Material* target_Ne = new G4Material("targetNe", density, ncomponents = 1, kStateGas, temperature, pressure);
+    target_Ne->AddElement(ele_Ne22, 100 * perCent);
 
     // ********** Sciintillator **********
     // Plastic scintillator 塑料闪烁体 EJ-200
@@ -106,7 +90,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 +119,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);
@@ -155,14 +144,13 @@ G4VPhysicalVolume* GeneDetectorConstruction::Construct() {
     // shielding PE
     G4double BPE_xyz = 60.0 * cm;
     G4Box* BPE_box = new G4Box("BPE_box", 0.5 * BPE_xyz, 0.5 * BPE_xyz, 0.5 * BPE_xyz);
-    // moderator PE
-    G4double PE_xyz = 50.0 * cm;
-    G4Box* PE_box = new G4Box("PE_box", 0.5 * PE_xyz, 0.5 * PE_xyz, 0.5 * PE_xyz);
+    // Plastic scintillator
+    G4double PS_xyz = 50.0 * cm;
+    G4Box* PS_box = new G4Box("PS_box", 0.5 * PS_xyz, 0.5 * PS_xyz, 0.5 * PS_xyz);
 
     // Beam hole
-    G4double rMin0 = 0.;
-    G4double rMax0 = 5.0 * cm;
-    G4Tubs* PE_hole = new G4Tubs("PE_hole", rMin0, rMax0, 0.5 * BPE_xyz, 0, 360. * deg);
+    G4double rMax = 5.0 * cm;
+    G4Tubs* PS_hole = new G4Tubs("PS_hole", 0, rMax, 0.5 * BPE_xyz, 0, 360. * deg);
 
     // He3 Tube holes
     G4double active_len = 30 * cm;  // length
@@ -202,48 +190,48 @@ G4VPhysicalVolume* GeneDetectorConstruction::Construct() {
     He3tub_log->SetVisAttributes(G4VisAttributes::GetInvisible);
     // He3active_log -> SetVisAttributes(G4VisAttributes::GetInvisible);
 
-    G4SubtractionSolid* solid_PE = new G4SubtractionSolid("solid_PE", PE_box, PE_hole);
-    G4SubtractionSolid* sub_BPE = new G4SubtractionSolid("sub_BPE", BPE_box, PE_box);
-    G4SubtractionSolid* solid_BPE = new G4SubtractionSolid("solid_BPE", sub_BPE, PE_hole);
+    G4SubtractionSolid* solid_PS = new G4SubtractionSolid("solid_PS", PS_box, PS_hole);
+    G4SubtractionSolid* sub_BPE = new G4SubtractionSolid("sub_BPE", BPE_box, PS_box);
+    G4SubtractionSolid* solid_BPE = new G4SubtractionSolid("solid_BPE", sub_BPE, PS_hole);
 
     char tmp[50];
     G4double posx, posy;
     for (G4int i = 0; i < 12; i++) {
         posx = 8.5 * cm * std::cos(i * 30.0 * deg);
         posy = 8.5 * cm * std::sin(i * 30.0 * deg);
-        sprintf(tmp, "solid_PE%d", i);
-        solid_PE = new G4SubtractionSolid(tmp, solid_PE, He3tube_Hole,
+        sprintf(tmp, "solid_PS%d", i);
+        solid_PS = new G4SubtractionSolid(tmp, solid_PS, He3tube_Hole,
                                           G4Transform3D(rot0, G4ThreeVector(posx, posy, tube_hole_z)));
         sprintf(tmp, "solid_BPE%d", i);
         solid_BPE = new G4SubtractionSolid(tmp, solid_BPE, He3tube_Hole,
                                            G4Transform3D(rot0, G4ThreeVector(posx, posy, tube_hole_z)));
         posx = 13.0 * cm * std::cos(i * 30.0 * deg + 15.0 * deg);
         posy = 13.0 * cm * std::sin(i * 30.0 * deg + 15.0 * deg);
-        sprintf(tmp, "solid_PE%d", i + 12);
-        solid_PE = new G4SubtractionSolid(tmp, solid_PE, He3tube_Hole,
+        sprintf(tmp, "solid_PS%d", i + 12);
+        solid_PS = new G4SubtractionSolid(tmp, solid_PS, He3tube_Hole,
                                           G4Transform3D(rot0, G4ThreeVector(posx, posy, tube_hole_z)));
         sprintf(tmp, "solid_BPE%d", i + 12);
         solid_BPE = new G4SubtractionSolid(tmp, solid_BPE, He3tube_Hole,
                                            G4Transform3D(rot0, G4ThreeVector(posx, posy, tube_hole_z)));
     }
     // z offset for position measurement
-    G4LogicalVolume* PE_log = new G4LogicalVolume(solid_PE, moderator_PE, "logic_PE");
-    new G4PVPlacement(0, G4ThreeVector(0, 0, Zoffset), "PE_phys", PE_log, physicalWorld, false, 1, true);
+    G4LogicalVolume* PS_log = new G4LogicalVolume(solid_PS, ej_200, "logic_PS");
+    new G4PVPlacement(0, G4ThreeVector(0, 0, Zoffset), "PS_phys", PS_log, physicalWorld, false, 1, true);
 
     G4LogicalVolume* BPE_log = new G4LogicalVolume(solid_BPE, shield_PE, "logic_BPE");
     new G4PVPlacement(0, G4ThreeVector(0, 0, Zoffset), "BPE_phys", BPE_log, physicalWorld, false, 2, true);
 
-    G4VisAttributes* PE_vis = new G4VisAttributes(G4Colour(0.0, 0.8, 0.8, 0.5));
-    PE_vis->SetForceSolid(true);
-    // PE_vis->SetForceWireframe(true);
-    PE_log->SetVisAttributes(PE_vis);
+    G4VisAttributes* PS_vis = new G4VisAttributes(G4Colour(0.0, 0.8, 0.8, 0.5));
+    PS_vis->SetForceSolid(true);
+    // PS_vis->SetForceWireframe(true);
+    PS_log->SetVisAttributes(PS_vis);
 
     // G4VisAttributes* BPE_vis= new G4VisAttributes(G4Colour(1.0,1.0,1.0,0.1));
     // BPE_vis->SetForceSolid(true);
     // BPE_vis->SetForceWireframe(true);
     // BPE_log->SetVisAttributes(BPE_vis);
 
-    // PE_log -> SetVisAttributes(G4VisAttributes::GetInvisible);
+    // PS_log -> SetVisAttributes(G4VisAttributes::GetInvisible);
 
     G4VisAttributes* BPE_vis = new G4VisAttributes(G4Colour(0.1, 0.3, 1, 0.8));
     BPE_log->SetVisAttributes(BPE_vis);
@@ -274,329 +262,26 @@ G4VPhysicalVolume* GeneDetectorConstruction::Construct() {
                           physicalWorld, false, 300 + i + 12, true);
     }
 
-    G4VisAttributes* C13_vis = new G4VisAttributes(G4Colour(1.0, 1.0, 1.0, 1.0));
-    C13_vis->SetForceSolid(true);
-    G4VisAttributes* CuSub_vis = new G4VisAttributes(G4Colour(0.1, 0.5, 0.7, 1.0));
-    CuSub_vis->SetForceSolid(true);
     G4VisAttributes* tube_vis = new G4VisAttributes(G4Colour(0.8, 0.8, 0.8, 1.0));
     tube_vis->SetForceSolid(true);
-    G4VisAttributes* water_vis = new G4VisAttributes(G4Colour(0.0, 0.0, 0.1, 1.0));
-    water_vis->SetForceSolid(true);
+    G4VisAttributes* Ne22_vis = new G4VisAttributes(G4Colour(1.0, 1.0, 1.0, 1.0));
+    Ne22_vis->SetForceSolid(true);
 
-    if (SReaction == "C13an") {
-        G4double C13thickness = 2.0 * mm;
-        G4Tubs* C13tub = new G4Tubs("C13tub", 0, 20. * mm, 0.5 * C13thickness, 0, 360. * deg);
+    // Target Ne22
+    G4double rBeam = rMax - 5 * mm;
+    G4double beam_len = 60.0 * cm;
+    G4double target_len = beam_len - 1 * cm;
+    G4Tubs* Ne22tub = new G4Tubs("Ne22tub", 0, rBeam, 0.5 * target_len, 0, 360. * deg);
+    G4LogicalVolume* Ne22log = new G4LogicalVolume(Ne22tub, target_Ne, "Ne22log");
+    new G4PVPlacement(0, G4ThreeVector(), "Ne22phys", Ne22log, physicalWorld, false, 3, true);
+    Ne22log->SetVisAttributes(Ne22_vis);
 
-        G4LogicalVolume* C13log = new G4LogicalVolume(C13tub, target_C, "C13log");
-        new G4PVPlacement(0, G4ThreeVector(), "C13phys", C13log, physicalWorld, false, 3, true);
-        C13log->SetVisAttributes(C13_vis);
-        // C13log -> SetVisAttributes(G4VisAttributes::GetInvisible);
-
-        // Cu_Substrate
-        G4double Cu0thickness = 1.5 * mm;
-        G4double Cu0z = 0.5 * (C13thickness + Cu0thickness);
-        G4Tubs* CuSubstrate_tub0 =
-            new G4Tubs("CuSubstrate_tub0", 18. * mm, 25. * mm, 0.5 * Cu0thickness, 0, 360. * deg);
-        G4LogicalVolume* CuSubstrate_log0 = new G4LogicalVolume(CuSubstrate_tub0, mat_Cu, "CuSubstrate_log0");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, Cu0z), "CuSubstrate_phys0", CuSubstrate_log0, physicalWorld, false, 4,
-                          true);
-
-        G4double Cu1thickness = C13thickness;
-        G4double Cu1z = 0.0 * mm;
-        G4Tubs* CuSubstrate_tub1 =
-            new G4Tubs("CuSubstrate_tub1", 20. * mm, 24.25 * mm, 0.5 * Cu1thickness, 0, 360. * deg);
-        G4LogicalVolume* CuSubstrate_log1 = new G4LogicalVolume(CuSubstrate_tub1, mat_Cu, "CuSubstrate_log1");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, Cu1z), "CuSubstrate_phys1", CuSubstrate_log1, physicalWorld, false, 5,
-                          true);
-
-        G4double Cu2thickness = 3.5 * mm;
-        G4double Cu2z = -0.5 * (C13thickness + Cu2thickness);
-        G4Tubs* CuSubstrate_tub2 = new G4Tubs("CuSubstrate_tub2", 0, 24.25 * mm, 0.5 * Cu2thickness, 0, 360. * deg);
-        G4LogicalVolume* CuSubstrate_log2 = new G4LogicalVolume(CuSubstrate_tub2, mat_Cu, "CuSubstrate_log2");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, Cu2z), "CuSubstrate_phys2", CuSubstrate_log2, physicalWorld, false, 6,
-                          true);
-
-        G4double Cu3thickness = 2.0 * mm;
-        G4double Cu3z = Cu2z - 0.5 * Cu2thickness - 0.5 * Cu3thickness;
-        G4Tubs* CuSubstrate_tub3 = new G4Tubs("CuSubstrate_tub3", 0, 29. * mm, 0.5 * Cu3thickness, 0, 360. * deg);
-        G4LogicalVolume* CuSubstrate_log3 = new G4LogicalVolume(CuSubstrate_tub3, mat_Cu, "CuSubstrate_log3");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, Cu3z), "CuSubstrate_phys3", CuSubstrate_log3, physicalWorld, false, 7,
-                          true);
-
-        CuSubstrate_log0->SetVisAttributes(CuSub_vis);
-        CuSubstrate_log1->SetVisAttributes(CuSub_vis);
-        CuSubstrate_log2->SetVisAttributes(CuSub_vis);
-        CuSubstrate_log3->SetVisAttributes(CuSub_vis);
-
-        // Target flange
-        G4double flange0_thickness = 9.0 * mm;
-        G4double flange0z = Cu2z - 0.5 * Cu2thickness + 0.5 * flange0_thickness;
-        G4Tubs* flange_tub0 = new G4Tubs("flange_tub0", 25. * mm, 42.33 * mm, 0.5 * flange0_thickness, 0, 360. * deg);
-        G4LogicalVolume* flange_log0 = new G4LogicalVolume(flange_tub0, mat_stell, "flange_log0");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, flange0z), "flange_phys0", flange_log0, physicalWorld, false, 8, true);
-
-        G4double flange1_thickness = 7.0 * mm;
-        G4double flange1z = flange0z - 0.5 * (flange0_thickness + flange1_thickness);
-        G4Tubs* flange_tub1 = new G4Tubs("flange_tub1", 29. * mm, 44.45 * mm, 0.5 * flange1_thickness, 0, 360. * deg);
-        G4LogicalVolume* flange_log1 = new G4LogicalVolume(flange_tub1, mat_stell, "flange_log1");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, flange1z), "flange_phys1", flange_log1, physicalWorld, false, 9, true);
-
-        // BeamLine_stainlesss
-        G4double BeampipeLength = 60.0 * cm;
-        G4double Beampipez = flange0z - 0.5 * flange0_thickness + 0.5 * BeampipeLength;
-        G4Tubs* Beampipe_tub =
-            new G4Tubs("Beampipe_tub", 42.33 * mm, 44.45 * mm, 0.5 * BeampipeLength, 0. * deg, 360. * deg);
-        G4LogicalVolume* Beampipe_log = new G4LogicalVolume(Beampipe_tub, mat_stell, "Beampipe_log");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, Beampipez), "Beampipe_phys", Beampipe_log, physicalWorld, false, 10,
-                          true);
-
-        flange_log0->SetVisAttributes(tube_vis);
-        flange_log1->SetVisAttributes(tube_vis);
-        Beampipe_log->SetVisAttributes(tube_vis);
-
-        // Cold trap
-        G4double CuTrapLength = 50.0 * cm;
-        G4double CuTrapz = 0.5 * CuTrapLength + 12.0 * mm + 60.0 * mm;
-        G4Tubs* CuTrap_tub = new G4Tubs("CuTrap_tub", 28.0 * mm, 30.0 * mm, 0.5 * CuTrapLength, 0, 360. * deg);
-        G4LogicalVolume* CuTrap_log = new G4LogicalVolume(CuTrap_tub, mat_Cu, "CuTrap_log");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, CuTrapz), "CuTrap_phys", CuTrap_log, physicalWorld, false, 11, true);
-
-        CuTrap_log->SetVisAttributes(CuSub_vis);
-
-        // Suppressor
-        G4double CuSupLength = 1.0 * mm;
-        G4double CuSupz = 0.5 * CuSupLength + 8.0 * mm + 60.0 * mm;
-        G4Tubs* CuSup_tub = new G4Tubs("CuSup_tub", 15.0 * mm, 30.0 * mm, 0.5 * CuSupLength, 0, 360. * deg);
-        G4LogicalVolume* CuSup_log = new G4LogicalVolume(CuSup_tub, mat_Cu, "CuSup_log");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, CuSupz), "CuSup_phys", CuSup_log, physicalWorld, false, 12, true);
-
-        // Collimator
-        G4double AuCollLength = 3.0 * um;
-        G4double AuCollz = -0.5 * AuCollLength + 12.0 * mm + 60.0 * mm;
-        G4Tubs* AuColl_tub = new G4Tubs("AuColl_tub", 5.0 * mm, 28.0 * mm, 0.5 * AuCollLength, 0, 360. * deg);
-        G4LogicalVolume* AuColl_log = new G4LogicalVolume(AuColl_tub, mat_Au, "AuColl_log");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, AuCollz), "AuColl_phys", AuColl_log, physicalWorld, false, 13, true);
-
-        G4double CuColl0Length = 2.0 * mm;
-        G4double CuColl0z = 0.5 * CuColl0Length + 12.0 * mm + 60.0 * mm;
-        G4Tubs* CuColl_tub0 = new G4Tubs("CuColl_tub0", 5.0 * mm, 28.0 * mm, 0.5 * CuColl0Length, 0, 360. * deg);
-        G4LogicalVolume* CuColl_log0 = new G4LogicalVolume(CuColl_tub0, mat_Cu, "CuColl_log0");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, CuColl0z), "CuColl_phys0", CuColl_log0, physicalWorld, false, 14,
-                          true);
-
-        G4double CuColl1Length = 3.0 * mm;
-        G4double CuColl1z = CuColl0z + 0.5 * CuColl0Length + 0.5 * CuColl1Length;
-        G4Tubs* CuColl_tub1 = new G4Tubs("CuColl_tub1", 25.0 * mm, 28.0 * mm, 0.5 * CuColl1Length, 0, 360. * deg);
-        G4LogicalVolume* CuColl_log1 = new G4LogicalVolume(CuColl_tub1, mat_Cu, "CuColl_log1");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, CuColl1z), "CuColl_phys1", CuColl_log1, physicalWorld, false, 15,
-                          true);
-
-        G4double VacuumLength = 70.0 * cm;
-        G4double Vacuumz = 0.5 * VacuumLength + Cu2z - 0.5 * Cu2thickness;
-        G4Tubs* Vacuum_tub = new G4Tubs("Vacuum_tub", 0, 42.33 * mm, 0.5 * VacuumLength, 0, 360. * deg);
-        G4double VacuumTflangez = -0.5 * VacuumLength + 0.5 * flange0_thickness;
-        G4SubtractionSolid* Vacuum_sub =
-            new G4SubtractionSolid("VacuumTflange_sub", Vacuum_tub, flange_tub0,
-                                   G4Transform3D(rot0, G4ThreeVector(0, 0, VacuumTflangez)));  // Target flange0
-        G4double VacuumCu2z = -0.5 * VacuumLength + 0.5 * Cu2thickness;
-        Vacuum_sub = new G4SubtractionSolid("Vacuum_sub2", Vacuum_sub, CuSubstrate_tub2,
-                                            G4Transform3D(rot0, G4ThreeVector(0, 0, VacuumCu2z)));  // Cu2
-        G4double VacuumCu1z = VacuumCu2z + 0.5 * Cu2thickness + 0.5 * Cu1thickness;
-        Vacuum_sub = new G4SubtractionSolid("Vacuum_sub1", Vacuum_sub, CuSubstrate_tub1,
-                                            G4Transform3D(rot0, G4ThreeVector(0, 0, VacuumCu1z)));  // Cu1
-        G4double VacuumCu0z = VacuumCu1z + 0.5 * Cu1thickness + 0.5 * Cu0thickness;
-        Vacuum_sub = new G4SubtractionSolid("Vacuum_sub0", Vacuum_sub, CuSubstrate_tub0,
-                                            G4Transform3D(rot0, G4ThreeVector(0, 0, VacuumCu0z)));  // Cu0
-        G4double VacuumC13z = VacuumCu1z;
-        Vacuum_sub = new G4SubtractionSolid("VacuumC13_sub", Vacuum_sub, C13tub,
-                                            G4Transform3D(rot0, G4ThreeVector(0, 0, VacuumC13z)));  // C13 target
-        G4double VacuumCuTrapz = VacuumC13z + 12.0 * mm + 0.5 * CuTrapLength + 60.0 * mm;
-        Vacuum_sub = new G4SubtractionSolid("VacuumCuTrap_sub", Vacuum_sub, CuTrap_tub,
-                                            G4Transform3D(rot0, G4ThreeVector(0, 0, VacuumCuTrapz)));  // Cu Trap
-        G4double VacuumCuSupz = VacuumC13z + 8.0 * mm + 0.5 * CuSupLength + 60.0 * mm;
-        Vacuum_sub = new G4SubtractionSolid("VacuumCuSup_sub", Vacuum_sub, CuSup_tub,
-                                            G4Transform3D(rot0, G4ThreeVector(0, 0, VacuumCuSupz)));  // Cu Suppressor
-        G4double VacuumAuCollz = VacuumC13z + 12.0 * mm - 0.5 * AuCollLength + 60.0 * mm;
-        Vacuum_sub = new G4SubtractionSolid("VacuumAuColl_sub", Vacuum_sub, AuColl_tub,
-                                            G4Transform3D(rot0, G4ThreeVector(0, 0, VacuumAuCollz)));  // Au Collimator
-        G4double VacuumCuColl0z = VacuumC13z + 12.0 * mm + 0.5 * CuColl0Length + 60.0 * mm;
-        Vacuum_sub =
-            new G4SubtractionSolid("VacuumCuColl_sub0", Vacuum_sub, CuColl_tub0,
-                                   G4Transform3D(rot0, G4ThreeVector(0, 0, VacuumCuColl0z)));  // Cu Collimator0
-        G4double VacuumCuColl1z = VacuumCuColl0z + 0.5 * CuColl0Length + 0.5 * CuColl1Length;
-        Vacuum_sub =
-            new G4SubtractionSolid("VacuumCuColl_sub1", Vacuum_sub, CuColl_tub1,
-                                   G4Transform3D(rot0, G4ThreeVector(0, 0, VacuumCuColl1z)));  // Cu Collimator1
-        G4LogicalVolume* Vacuum_log = new G4LogicalVolume(Vacuum_sub, Vacuum, "Vacuum_log");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, Vacuumz), "Vacuum_phys", Vacuum_log, physicalWorld, false, 16, true);
-
-        // Water_cold
-        G4double WaterBack_thickness = 1.0 * mm;
-        G4double WaterBackz = Cu3z - 0.5 * Cu3thickness - 0.5 * WaterBack_thickness;
-        G4Tubs* WaterBack_tub =
-            new G4Tubs("WaterBack_tub", 0, 29.0 * mm, 0.5 * WaterBack_thickness, 0. * deg, 360. * deg);
-        G4LogicalVolume* WaterBack_log = new G4LogicalVolume(WaterBack_tub, mat_Water, "WaterBack_log");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, WaterBackz), "WaterBack_phys", WaterBack_log, physicalWorld, false, 17,
-                          true);
-        WaterBack_log->SetVisAttributes(water_vis);
-
-        // Al Substrate
-        G4double Al0thickness = 4.0 * mm;
-        G4double Al0z = WaterBackz - 0.5 * WaterBack_thickness - 0.5 * Al0thickness;
-        G4Tubs* AlSubstrate_tub0 =
-            new G4Tubs("AlSubstrate_tub0", 0, 29.0 * mm, 0.5 * Al0thickness, 0. * deg, 360. * deg);
-        G4LogicalVolume* AlSubstrate_log0 = new G4LogicalVolume(AlSubstrate_tub0, mat_Al, "AlSubstrate_log0");
-        G4PVPlacement* Al0phys = new G4PVPlacement(0, G4ThreeVector(0, 0, Al0z), "AlSubstrate_phys0", AlSubstrate_log0,
-                                                   physicalWorld, false, 18, true);
-        AlSubstrate_log0->SetVisAttributes(tube_vis);
-
-        G4double Al1thickness = 4.0 * mm;
-        G4double Al1z = Al0z - 0.5 * Al0thickness - 0.5 * Al1thickness;
-        G4Tubs* AlSubstrate_tub1 =
-            new G4Tubs("AlSubstrate_tub1", 0, 44.45 * mm, 0.5 * Al1thickness, 0. * deg, 360. * deg);
-        G4LogicalVolume* AlSubstrate_log1 = new G4LogicalVolume(AlSubstrate_tub1, mat_Al, "AlSubstrate_log1");
-        G4PVPlacement* Al1phys = new G4PVPlacement(0, G4ThreeVector(0, 0, Al1z), "AlSubstrate_phys1", AlSubstrate_log1,
-                                                   physicalWorld, false, 19, true);
-        AlSubstrate_log1->SetVisAttributes(tube_vis);
-
-        // Waterpipe
-        G4double Waterpipe0_thickness = 4.0 * mm;
-        G4Tubs* Waterpipe_tub0 =
-            new G4Tubs("Waterpipe_tub0", 0, 2.5 * mm, 0.5 * Waterpipe0_thickness, 0. * deg, 360. * deg);
-        G4LogicalVolume* Waterpipe_log0 = new G4LogicalVolume(Waterpipe_tub0, mat_Water, "Waterpipe_log0");
-        new G4PVPlacement(0, G4ThreeVector(0, 21.5 * mm, 0), "Waterpipe_phys00", Waterpipe_log0, Al0phys, false, 20,
-                          true);
-        new G4PVPlacement(0, G4ThreeVector(0, -21.5 * mm, 0), "Waterpipe_phys01", Waterpipe_log0, Al0phys, false, 21,
-                          true);
-        new G4PVPlacement(0, G4ThreeVector(0, 21.5 * mm, 0), "Waterpipe_phys02", Waterpipe_log0, Al1phys, false, 22,
-                          true);
-        new G4PVPlacement(0, G4ThreeVector(0, -21.5 * mm, 0), "Waterpipe_phys03", Waterpipe_log0, Al1phys, false, 23,
-                          true);
-
-        G4double Waterpipe1_thickness = BPE_xyz / 2.0 + 10.0 * cm;
-        G4double Waterpipe1z = Al1z - 0.5 * Al1thickness - 0.5 * Waterpipe1_thickness;
-        G4Tubs* Waterpipe_tub1 =
-            new G4Tubs("Waterpipe_tub1", 0, 2.5 * mm, 0.5 * Waterpipe1_thickness, 0. * deg, 360. * deg);
-        G4LogicalVolume* Waterpipe_log1 = new G4LogicalVolume(Waterpipe_tub1, mat_Water, "Waterpipe_log1");
-        new G4PVPlacement(0, G4ThreeVector(0, 21.5 * mm, Waterpipe1z), "Waterpipe_phys1", Waterpipe_log1, physicalWorld,
-                          false, 24, true);
-        new G4PVPlacement(0, G4ThreeVector(0, -21.5 * mm, Waterpipe1z), "Waterpipe_phys2", Waterpipe_log1,
-                          physicalWorld, false, 25, true);
-        Waterpipe_log0->SetVisAttributes(water_vis);
-        Waterpipe_log1->SetVisAttributes(water_vis);
-
-    } else if (SReaction == "V51pn") {
-        // V51 on Ta target
-        G4double Tathickness = 1.0 * mm;
-        G4Tubs* Tatub = new G4Tubs("Tatub", 0, 15. * mm, 0.5 * Tathickness, 0, 360. * deg);
-
-        G4LogicalVolume* Talog = new G4LogicalVolume(Tatub, mat_Ta, "Talog");
-        new G4PVPlacement(0, G4ThreeVector(), "Taphys", Talog, physicalWorld, false, 3, true);
-        Talog->SetVisAttributes(C13_vis);
-
-        // Cu_Substrate
-        G4double Cu0thickness = 8.0 * mm;
-        G4double Cu0z = 0.5 * (Tathickness + Cu0thickness);
-        G4Tubs* CuSubstrate_tub0 =
-            new G4Tubs("CuSubstrate_tub0", 11. * mm, 29. * mm, 0.5 * Cu0thickness, 0, 360. * deg);
-        G4LogicalVolume* CuSubstrate_log0 = new G4LogicalVolume(CuSubstrate_tub0, mat_Cu, "CuSubstrate_log0");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, Cu0z), "CuSubstrate_phys1", CuSubstrate_log0, physicalWorld, false, 4,
-                          true);
-
-        CuSubstrate_log0->SetVisAttributes(CuSub_vis);
-
-        // Target flange
-        G4double flange0_thickness = 9.0 * mm;
-        G4double flange0z = Cu0z + 0.5 * Cu0thickness + 0.5 * flange0_thickness;
-        G4Tubs* flange_tub0 = new G4Tubs("flange_tub0", 25. * mm, 42.33 * mm, 0.5 * flange0_thickness, 0, 360. * deg);
-        G4LogicalVolume* flange_log0 = new G4LogicalVolume(flange_tub0, mat_stell, "flange_log0");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, flange0z), "flange_phys0", flange_log0, physicalWorld, false, 8, true);
-
-        G4double flange1_thickness = 7.0 * mm;
-        G4double flange1z = flange0z - 0.5 * (flange0_thickness + flange1_thickness);
-        G4Tubs* flange_tub1 = new G4Tubs("flange_tub1", 29. * mm, 44.45 * mm, 0.5 * flange1_thickness, 0, 360. * deg);
-        G4LogicalVolume* flange_log1 = new G4LogicalVolume(flange_tub1, mat_stell, "flange_log1");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, flange1z), "flange_phys1", flange_log1, physicalWorld, false, 9, true);
-
-        // BeamLine_stainlesss
-        G4double BeampipeLength = 60.0 * cm;
-        G4double Beampipez = flange0z - 0.5 * flange0_thickness + 0.5 * BeampipeLength;
-        G4Tubs* Beampipe_tub =
-            new G4Tubs("Beampipe_tub", 42.33 * mm, 44.45 * mm, 0.5 * BeampipeLength, 0. * deg, 360. * deg);
-        G4LogicalVolume* Beampipe_log = new G4LogicalVolume(Beampipe_tub, mat_stell, "Beampipe_log");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, Beampipez), "Beampipe_phys", Beampipe_log, physicalWorld, false, 10,
-                          true);
-
-        flange_log0->SetVisAttributes(tube_vis);
-        flange_log1->SetVisAttributes(tube_vis);
-        Beampipe_log->SetVisAttributes(tube_vis);
-
-        G4double VacuumLength = 70.0 * cm;
-        G4double Vacuumz = 0.5 * VacuumLength + Cu0z - 0.5 * Cu0thickness;
-        G4Tubs* Vacuum_tub = new G4Tubs("Vacuum_tub", 0, 42.33 * mm, 0.5 * VacuumLength, 0, 360. * deg);
-        G4double VacuumCu0z = -0.5 * VacuumLength + 0.5 * Cu0thickness;
-        G4Tubs* CuSubstrate_tub1 = new G4Tubs("CuSubstrate_tub1", 11. * mm, 50 * mm, 0.5 * Cu0thickness, 0, 360. * deg);
-        G4SubtractionSolid* Vacuum_sub = new G4SubtractionSolid(
-            "Vacuum_sub0", Vacuum_tub, CuSubstrate_tub1, G4Transform3D(rot0, G4ThreeVector(0, 0, VacuumCu0z)));  // Cu0
-        G4double Vacuumflange0z = VacuumCu0z + 0.5 * Cu0thickness + 0.5 * flange0_thickness;
-        Vacuum_sub =
-            new G4SubtractionSolid("Vacuumflange0_sub", Vacuum_sub, flange_tub0,
-                                   G4Transform3D(rot0, G4ThreeVector(0, 0, Vacuumflange0z)));  // Target flange0
-        G4LogicalVolume* Vacuum_log = new G4LogicalVolume(Vacuum_sub, Vacuum, "Vacuum_log");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, Vacuumz), "Vacuum_phys", Vacuum_log, physicalWorld, false, 11, true);
-
-        Vacuum_log->SetVisAttributes(G4VisAttributes::GetInvisible);
-        // Water_cold
-        G4double WaterBack_thickness = 1.3 * mm;
-        G4double WaterBackz = -0.5 * Tathickness - 0.5 * WaterBack_thickness;
-        G4Tubs* WaterBack_tub =
-            new G4Tubs("WaterBack_tub", 0, 12.0 * mm, 0.5 * WaterBack_thickness, 0. * deg, 360. * deg);
-        G4LogicalVolume* WaterBack_log = new G4LogicalVolume(WaterBack_tub, mat_Water, "WaterBack_log");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, WaterBackz), "WaterBack_phys", WaterBack_log, physicalWorld, false, 12,
-                          true);
-
-        WaterBack_log->SetVisAttributes(water_vis);
-
-        // Fe backflange
-        G4double Fe0thickness = 1.3 * mm;
-        G4double Fe0z = WaterBackz;
-        G4Tubs* FeBack_tub0 =
-            new G4Tubs("FeBack_tub0", 12.0 * mm, 44.45 * mm, 0.5 * Fe0thickness, 0. * deg, 360. * deg);
-        G4LogicalVolume* FeBack_log0 = new G4LogicalVolume(FeBack_tub0, mat_stell, "FeBack_log0");
-        new G4PVPlacement(0, G4ThreeVector(0, 0, Fe0z), "FeBack_phys0", FeBack_log0, physicalWorld, false, 13, true);
-
-        G4double Fe1thickness = 6.7 * mm;
-        G4double Fe1z = Fe0z - 0.5 * Fe0thickness - 0.5 * Fe1thickness;
-        G4Tubs* FeBack_tub1 = new G4Tubs("FeBack_tub1", 0, 44.45 * mm, 0.5 * Fe1thickness, 0. * deg, 360. * deg);
-        G4LogicalVolume* FeBack_log1 = new G4LogicalVolume(FeBack_tub1, mat_stell, "FeBack_log1");
-        G4PVPlacement* Fe1phys = new G4PVPlacement(0, G4ThreeVector(0, 0, Fe1z), "FeBack_phys1", FeBack_log1,
-                                                   physicalWorld, false, 14, true);
-
-        FeBack_log0->SetVisAttributes(tube_vis);
-        FeBack_log1->SetVisAttributes(tube_vis);
-
-        // Waterpipe
-        G4double Waterpipe1_thickness = Fe1thickness;
-        G4Tubs* Waterpipe_tub1 =
-            new G4Tubs("Waterpipe_tub1", 0, 1.6 * mm, 0.5 * Waterpipe1_thickness, 0. * deg, 360. * deg);
-        G4LogicalVolume* Waterpipe_log1 = new G4LogicalVolume(Waterpipe_tub1, mat_Water, "Waterpipe_log1");
-        new G4PVPlacement(0, G4ThreeVector(0, 9.0 * mm, 0), "Waterpipe_phys01", Waterpipe_log1, Fe1phys, false, 17,
-                          true);
-        new G4PVPlacement(0, G4ThreeVector(0, -9.0 * mm, 0), "Waterpipe_phys02", Waterpipe_log1, Fe1phys, false, 18,
-                          true);
-
-        G4double Waterpipe2_thickness = BPE_xyz / 2.0 + 10.0 * cm;
-        G4double Waterpipe2z = Fe1z - 0.5 * Fe1thickness - 0.5 * Waterpipe2_thickness;
-        G4Tubs* Waterpipe_tub2 =
-            new G4Tubs("Waterpipe_tub2", 0, 2.0 * mm, 0.5 * Waterpipe2_thickness, 0. * deg, 360. * deg);
-        G4LogicalVolume* Waterpipe_log2 = new G4LogicalVolume(Waterpipe_tub2, mat_Water, "Waterpipe_log2");
-        new G4PVPlacement(0, G4ThreeVector(0, 9.0 * mm, Waterpipe2z), "Waterpipe_phys20", Waterpipe_log2, physicalWorld,
-                          false, 19, true);
-        new G4PVPlacement(0, G4ThreeVector(0, -9.0 * mm, Waterpipe2z), "Waterpipe_phys21", Waterpipe_log2,
-                          physicalWorld, false, 20, true);
-        Waterpipe_log1->SetVisAttributes(water_vis);
-        Waterpipe_log2->SetVisAttributes(water_vis);
-    }
+    // Beam Line Tube Stainlesss 不锈钢管
+    G4Tubs* beamtub = new G4Tubs("Beam_tub", 0, rMax, 0.5 * beam_len, 0. * deg, 360. * deg);
+    G4SubtractionSolid* beamtub_sub = new G4SubtractionSolid("Beam_sub", beamtub, Ne22tub);
+    G4LogicalVolume* beam_log = new G4LogicalVolume(beamtub_sub, mat_stell, "beam_log");
+    new G4PVPlacement(0, G4ThreeVector(0, 0, 0), "beam_phys", beam_log, physicalWorld, false, 10, true);
+    beam_log->SetVisAttributes(tube_vis);
 
     G4SDManager* SDman = G4SDManager::GetSDMpointer();
     GeneHe3detSD* detSD = new GeneHe3detSD("/GeneHe3det/");

src/GeneDetectorConstructionMessenger.cc

@@ -24,24 +24,12 @@ GeneDetectorConstructionMessenger::GeneDetectorConstructionMessenger(GeneDetecto
     ZactiveshiftCmd->SetUnitCategory("Length");
     ZactiveshiftCmd->SetRange("Zactiveshift>-20 && Zactiveshift<20");
     ZactiveshiftCmd->AvailableForStates(G4State_PreInit, G4State_Idle);
-
-    PE_BfractionCmd = new G4UIcmdWithADouble("/Gene/Detector/SetBfraction", this);
-    PE_BfractionCmd->SetGuidance("Set Boron fraction of PE moderator");
-    PE_BfractionCmd->SetParameterName("PE_Bfraction", false);
-    PE_BfractionCmd->SetRange("PE_Bfraction>=0");
-    PE_BfractionCmd->AvailableForStates(G4State_PreInit, G4State_Idle);
-
-    DetReactionCmd = new G4UIcmdWithAString("/Gene/Detector/SetDetReaction", this);
-    DetReactionCmd->SetGuidance("Set type of reaction for detector setup.");
-    DetReactionCmd->SetParameterName("DetReaction", false);
-    DetReactionCmd->AvailableForStates(G4State_PreInit);
 }
 
 GeneDetectorConstructionMessenger::~GeneDetectorConstructionMessenger() {
     delete ZoffsetCmd;
     delete ZactiveshiftCmd;
     delete PE_BfractionCmd;
-    delete DetReactionCmd;
     delete DetectorDir;
 }
 
@@ -52,10 +40,4 @@ void GeneDetectorConstructionMessenger::SetNewValue(G4UIcommand* command, G4Stri
     if (command == ZactiveshiftCmd) {
         pDetectorConstruction->SetZactiveshift(ZactiveshiftCmd->GetNewDoubleValue(newValue));
     }
-    if (command == PE_BfractionCmd) {
-        pDetectorConstruction->SetBfraction(PE_BfractionCmd->GetNewDoubleValue(newValue));
-    }
-    if (command == DetReactionCmd) {
-        pDetectorConstruction->SetDetReaction(newValue);
-    }
 }

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"
 
@@ -38,52 +38,20 @@ GenePrimaryGeneratorAction::~GenePrimaryGeneratorAction() {
     delete pMessenger;
 }
 
-void GenePrimaryGeneratorAction::SetReaction(G4String& name) {
-    if (name == "V51pn")
-        iReaction = 0;
-    else if (name == "C13an")
-        iReaction = 1;
-    else {
-        G4cout << "Unknown reaction type!!!!" << G4endl;
-        G4cout << "Unknown reaction type!!!!" << G4endl;
-        G4cout << "Unknown reaction type!!!!" << G4endl;
-        G4cout << "Unknown reaction type!!!!" << G4endl;
-        G4cout << "Unknown reaction type!!!!" << G4endl;
-        G4cin >> iReaction;
-    }
-}
-
 // ???
 void GenePrimaryGeneratorAction::SetBeamEnergy(G4double val) {
     fBeamEnergy = val / MeV;
 
     G4double umass = 931.494;
     G4double massHe4 = 4.00260325413 * umass;
-    G4double massC13 = 13.00335483521 * umass;
-    G4double massO16 = 15.9949146196 * umass;
+    G4double massNe22 = 21.991385114 * umass;
+    G4double massMg25 = 24.98583697 * umass;
     G4double massn = 939.56563;
-    G4double massp = 938.78294;
-    G4double massV51 = 50.9439569 * umass;
-    G4double massCr51 = 50.9447647 * umass;
 
-    if (iReaction == 0) {
-        mass[0] = massp;
-        mass[1] = massV51;
-        mass[2] = massCr51;
-        mass[3] = massn;
-    } else if (iReaction == 1) {
-        mass[0] = massHe4;
-        mass[1] = massC13;
-        mass[2] = massO16;
-        mass[3] = massn;
-    } else {
-        G4cout << "Unknown reaction type!!!!" << G4endl;
-        G4cout << "Unknown reaction type!!!!" << G4endl;
-        G4cout << "Unknown reaction type!!!!" << G4endl;
-        G4cout << "Unknown reaction type!!!!" << G4endl;
-        G4cout << "Unknown reaction type!!!!" << G4endl;
-        G4cin >> iReaction;
-    }
+    mass[0] = massHe4;
+    mass[1] = massNe22;
+    mass[2] = massMg25;
+    mass[3] = massn;
 
     G4double eb = fBeamEnergy + mass[0];
     G4double pb2 = fBeamEnergy * fBeamEnergy + 2.0 * fBeamEnergy * mass[0];
@@ -96,7 +64,6 @@ void GenePrimaryGeneratorAction::SetBeamEnergy(G4double val) {
     G4double t_cm = e_cm - mass[2] - mass[3];
     if (t_cm < 0.0) {
         G4cout << "No solution!" << G4endl;
-        G4cin >> iReaction;
         return;
     }
     G4double t_cm2 = t_cm * t_cm;

src/GenePrimaryGeneratorActionMessenger.cc

@@ -22,11 +22,6 @@ GenePrimaryGeneratorActionMessenger::GenePrimaryGeneratorActionMessenger(GenePri
     AngularCmd->SetGuidance("Set type of angular distribution.");
     AngularCmd->SetParameterName("Angular", false);
     AngularCmd->AvailableForStates(G4State_PreInit);
-
-    ReactionCmd = new G4UIcmdWithAString("/Gene/PrimaryGA/SetReaction", this);
-    ReactionCmd->SetGuidance("Set type of reaction.");
-    ReactionCmd->SetParameterName("Reaction", false);
-    ReactionCmd->AvailableForStates(G4State_PreInit);
 }
 
 GenePrimaryGeneratorActionMessenger::~GenePrimaryGeneratorActionMessenger() {
@@ -44,8 +39,4 @@ void GenePrimaryGeneratorActionMessenger::SetNewValue(G4UIcommand* command, G4St
     if (command == AngularCmd) {
         pPrimaryGeneratorAction->SetAngType(newValue);
     }
-
-    if (command == ReactionCmd) {
-        pPrimaryGeneratorAction->SetReaction(newValue);
-    }
 }

vis.mac

@@ -21,7 +21,7 @@
 # 设置颜色
 /vis/viewer/set/style surface
 # 设置 World 不可见
-/vis/geometry/set/visibility world 0 false
+/vis/geometry/set/visibility World 0 false
 
 # 坐标轴
 /vis/scene/add/scale   # Simple scale line