fix: physics

include/GeneAnalysisManager.hh

@@ -44,7 +44,6 @@ private:
     static GeneAnalysisManager *instance;
     struct detEvent mydata;
 
-    // G4int evtNo;
     G4String analysisFileName;
 
     TFile *theTFile;

include/GeneDetectorConstructionMessenger.hh

@@ -23,7 +23,6 @@ private:
     G4UIdirectory* DetectorDir;
     G4UIcmdWithADoubleAndUnit* ZoffsetCmd;
     G4UIcmdWithADoubleAndUnit* ZactiveshiftCmd;
-    G4UIcmdWithADouble* PE_BfractionCmd;
     G4UIcmdWithAString* DetReactionCmd;
 };
 

include/GeneHe3detSD.hh

@@ -22,7 +22,6 @@ public:
 
 private:
     GeneHe3detHitsCollection* scintillatorCollection;
-    G4int HitID;
 };
 
 #endif

include/GenePhysicsList.hh

@@ -1,36 +0,0 @@
-#ifndef GenePhysicsList_hh
-#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

@@ -16,7 +16,6 @@ public:
 
     void GeneratePrimaries(G4Event* anEvent);
     void SetBeamEnergy(G4double);
-    void SetAngType(G4String&);
 
 private:
     G4ParticleGun* particleGun;

include/GenePrimaryGeneratorActionMessenger.hh

@@ -21,8 +21,6 @@ private:
 
     G4UIdirectory* PrimaryDir;
     G4UIcmdWithADoubleAndUnit* BeamEnergyCmd;
-    G4UIcmdWithAString* AngularCmd;
-    G4UIcmdWithAString* ReactionCmd;
 };
 
 #endif

init_vis.mac

@@ -1,16 +1,10 @@
-# Macro file for the initialization of example B1
+/control/verbose 1
-# in interactive session
-#
-# Set some default verbose
-/control/verbose 2
-/control/saveHistory
 /run/verbose 2
-#
+/event/verbose 0
-# Change the default number of threads (in multi-threaded mode)
+/tracking/verbose 0
-#/run/numberOfThreads 4
+
-#
 # Initialize kernel
 /run/initialize
-#
+
 # Visualization setting
 /control/execute vis.mac

main.cpp

@@ -1,3 +1,4 @@
+// clang-format off
 #include "TROOT.h"
 #include "time.h"
 
@@ -9,11 +10,14 @@
 #include "G4UIterminal.hh"
 #include "G4UnitsTable.hh"
 #include "G4VisExecutive.hh"
+#include "G4OpticalPhysics.hh"
+#include "QGSP_BERT_HP.hh"
+
 #include "GeneDetectorConstruction.hh"
 #include "GeneEventAction.hh"
-#include "GenePhysicsList.hh"
 #include "GenePrimaryGeneratorAction.hh"
 #include "GeneRunAction.hh"
+// clang-format on
 
 int main(int argc, char **argv) {
     // random engine
@@ -32,7 +36,9 @@ int main(int argc, char **argv) {
     G4RunManager *runManager = new G4RunManager;
 
     runManager->SetUserInitialization(new GeneDetectorConstruction);
-    runManager->SetUserInitialization(new GenePhysicsList);
+    G4VModularPhysicsList *physics = new QGSP_BERT_HP;
+    physics->RegisterPhysics(new G4OpticalPhysics);
+    runManager->SetUserInitialization(physics);
 
     runManager->SetUserAction(new GenePrimaryGeneratorAction);
     runManager->SetUserAction(new GeneRunAction);
@@ -48,7 +54,7 @@ int main(int argc, char **argv) {
         G4String fileName = argv[1];
         UIManager->ApplyCommand(command + fileName);
     } else {
-        UIManager->ApplyCommand("/control/execute vis.mac");
+        UIManager->ApplyCommand("/control/execute init_vis.mac");
         ui->SessionStart();
         delete ui;
     }

run.mac

@@ -1,35 +1,25 @@
-# Macro file for the initialization phase of "Neutron.cc"
-# Sets some default verbose and initializes the graphic.
-#
-#
-# Some Settings
 /control/verbose 1
-/run/verbose 0
+/run/verbose 2
+/event/verbose 0
 /tracking/verbose 0
 
 # user setting
-# 1st reaction
+# 1st energy
-# 2nd energy
+# 2nd offset
-# 3rd angular
+# 3rd set particle
 
 # energy setting
 /Gene/PrimaryGA/SetBeamEnergy 1.0 MeV       # alpha energy for 22Ne(a,n)25Mg reaction
 
-# angular setting
-/Gene/PrimaryGA/SetAngular ISO              # ISO or ENDF
-
 # detector offset setting
 /Gene/Detector/SetZoffset 0.0 mm           # default: 0 mm
 /Gene/Detector/SetZactiveshift 0.0 mm      # default: 0 mm
 
-# Boron fraction in PE moderator
-/Gene/Detector/SetBfraction 0.054           # perCent, default: 0.054%
-
 # user setting for runNo of output file
-/Gene/Run/RunNo 00_Ne22an_ISO_1000keV_zoffset0mm_zshift0mm_B054_100000evts_Coll60mm
+# /Gene/Run/RunNo 00_Ne22an_1000keV_zoffset0mm_zshift0mm_10000evts
 
 /run/initialize
 
 # gsp particle
 /gun/particle neutron
-/run/beamOn 1000
+/run/beamOn 10000

src/GeneAnalysisManager.cc

@@ -3,7 +3,6 @@
 GeneAnalysisManager* GeneAnalysisManager::instance = 0;
 
 GeneAnalysisManager::GeneAnalysisManager() : analysisFileName("Simulation_Ne22an.root"), theTFile(0), tr(0) {
-    // evtNo = 0;
     memset(&mydata, 0x00, sizeof(mydata));
 }
 
@@ -34,8 +33,6 @@ void GeneAnalysisManager::FillAng(G4double thetacm, G4double thetalab) {
     mydata.Thetalab = thetalab;
 }
 
-// void GeneAnalysisManager::EventNo() { evtNo++; }
-
 void GeneAnalysisManager::book() {
     // delete all associated variables created via new, moreover it delete itself.
     if (theTFile != 0) delete theTFile;
@@ -62,6 +59,5 @@ void GeneAnalysisManager::save() {
         theTFile->Write();
         theTFile->Close();
     }
-    // evtNo = 0;
     // 4cout<<" --> reactNum: "<< reactNum <<G4endl;
 }

src/GeneDetectorConstructionMessenger.cc

@@ -29,7 +29,6 @@ GeneDetectorConstructionMessenger::GeneDetectorConstructionMessenger(GeneDetecto
 GeneDetectorConstructionMessenger::~GeneDetectorConstructionMessenger() {
     delete ZoffsetCmd;
     delete ZactiveshiftCmd;
-    delete PE_BfractionCmd;
     delete DetectorDir;
 }
 

src/GeneEventAction.cc

@@ -23,15 +23,14 @@ GeneEventAction::~GeneEventAction() {}
 void GeneEventAction::BeginOfEventAction(const G4Event* evt) {
     G4int evtNb = evt->GetEventID();
     if (evtNb % 1000 == 0) G4cout << "Event No." << evtNb << G4endl;
-    G4SDManager* SDman = G4SDManager::GetSDMpointer();
 
+    G4SDManager* SDman = G4SDManager::GetSDMpointer();
     if (GeneHe3detSDCollID == -1) {
         GeneHe3detSDCollID = SDman->GetCollectionID("GeneHe3detHitCollection");
     }
 }
 
 void GeneEventAction::EndOfEventAction(const G4Event* evt) {
-    // G4int event_id = evt->GetEventID();
     G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
     GeneHe3detHitsCollection* GeneHe3detHC = 0;
     GeneAnalysisManager* analysisManager = GeneAnalysisManager::GetInstance();

src/GeneHe3detSD.cc

@@ -17,21 +17,20 @@ GeneHe3detSD::GeneHe3detSD(G4String name) : G4VSensitiveDetector(name) {
 
 GeneHe3detSD::~GeneHe3detSD() {}
 
-void GeneHe3detSD::Initialize(G4HCofThisEvent*) {
+void GeneHe3detSD::Initialize(G4HCofThisEvent* HCE) {
     scintillatorCollection = new GeneHe3detHitsCollection(SensitiveDetectorName, collectionName[0]);
-
+    static G4int HCID = -1;
-    HitID = -1;
+    if (HCID < 0) HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
+    HCE->AddHitsCollection(HCID, scintillatorCollection);
 }
 
 G4bool GeneHe3detSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
     G4double edep = aStep->GetTotalEnergyDeposit();
+    if (edep == 0.) return false;
+
     G4ParticleDefinition* particleType = aStep->GetTrack()->GetDefinition();
     G4String particleName = particleType->GetParticleName();
-    // G4double stepl = 0.;
-    // if (particleType->GetPDGCharge() != 0.)
-    // stepl = aStep->GetStepLength();
     if (particleName == "opticalphoton") return false;
-    if ((edep == 0.)) return false;
 
     G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
     G4TouchableHistory* theTouchable = (G4TouchableHistory*)(preStepPoint->GetTouchable());
@@ -46,17 +45,11 @@ G4bool GeneHe3detSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
     newHit->SetParticle(particleName);
     newHit->SetParticleEnergy(aStep->GetPreStepPoint()->GetKineticEnergy());
 
-    HitID = scintillatorCollection->insert(newHit);
+    scintillatorCollection->insert(newHit);
-
     return true;
 }
 
 void GeneHe3detSD::EndOfEvent(G4HCofThisEvent* HCE) {
-    G4String HCname = collectionName[0];
-    static G4int HCID = -1;
-    if (HCID < 0) HCID = G4SDManager::GetSDMpointer()->GetCollectionID(HCname);
-    HCE->AddHitsCollection(HCID, scintillatorCollection);
-
     G4int nHits = scintillatorCollection->entries();
     if (verboseLevel >= 1) G4cout << "     Si collection: " << nHits << " hits" << G4endl;
     if (verboseLevel >= 2) scintillatorCollection->PrintAllHits();

src/GenePhysicsList.cc

@@ -1,746 +0,0 @@
-#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

@@ -27,6 +27,12 @@ GenePrimaryGeneratorAction::GenePrimaryGeneratorAction() {
 
     // 通过角分布采样获取随机数
     hAng = new TH1F("hAng", "hAng;theta(deg);prob.", 180, 0, 180);
+    G4double thcm, pthcm;
+    for (G4int i = 0; i < 180; i++) {
+        thcm = pi * i / 180.0;
+        pthcm = 2.0 * pi * std::sin(thcm);
+        hAng->SetBinContent(i + 1, pthcm);
+    }
 
     fBeamEnergy = 0;
     pMessenger = new GenePrimaryGeneratorActionMessenger(this);
@@ -75,65 +81,6 @@ void GenePrimaryGeneratorAction::SetBeamEnergy(G4double val) {
     p4_cm = std::sqrt(p4_cm2);
 }
 
-// ???
-void GenePrimaryGeneratorAction::SetAngType(G4String& name) {
-    if (name == "ENDF") {
-        std::ifstream infile("xt_nov21.a");
-        G4int iN = 0;
-        if (infile.fail()) {
-            G4cout << "ifstream: Failure to open input file!!!" << G4endl;
-            G4cout << "ifstream: Failure to open input file!!!" << G4endl;
-            G4cout << "ifstream: Failure to open input file!!!" << G4endl;
-            G4cout << "ifstream: Failure to open input file!!!" << G4endl;
-            G4cout << "ifstream: Failure to open input file!!!" << G4endl;
-            G4cin >> iN;
-            return;
-        }
-        G4double Elab[5000], A1[5000], A2[5000], A3[5000];
-        G4double felab, fecm, fa0, fa1, fa2, fa3;
-        while (infile >> felab >> fecm >> fa0 >> fa1 >> fa2 >> fa3) {
-            Elab[iN] = felab;
-            A1[iN] = fa1;
-            A2[iN] = fa2;
-            A3[iN] = fa3;
-            iN++;
-        }
-        infile.close();
-        TGraph gr1(iN, Elab, A1);
-        TGraph gr2(iN, Elab, A2);
-        TGraph gr3(iN, Elab, A3);
-        fa1 = gr1.Eval(fBeamEnergy, 0, "S");
-        fa2 = gr2.Eval(fBeamEnergy, 0, "S");
-        fa3 = gr3.Eval(fBeamEnergy, 0, "S");
-        G4double thcm, costh, pthcm;
-        for (G4int i = 0; i < 180; i++) {
-            thcm = pi * i / 180.0;
-            costh = std::cos(thcm);
-            pthcm = 1.0 + fa1 * costh + fa2 * 0.5 * (3.0 * costh * costh - 1.0) +
-                    fa3 * 0.5 * (5.0 * std::pow(costh, 3.0) - 3.0 * costh);
-            pthcm = pthcm * 2.0 * pi * std::sin(thcm);
-            if (pthcm < 0.0) pthcm = 0.0;
-            hAng->SetBinContent(i + 1, pthcm);
-        }
-    } else if (name == "ISO") {
-        G4double thcm, pthcm;
-        for (G4int i = 0; i < 180; i++) {
-            thcm = pi * i / 180.0;
-            pthcm = 2.0 * pi * std::sin(thcm);
-            hAng->SetBinContent(i + 1, pthcm);
-        }
-    } else {
-        G4cout << "Wrong Angular type!!!!" << G4endl;
-        G4cout << "Wrong Angular type!!!!" << G4endl;
-        G4cout << "Wrong Angular type!!!!" << G4endl;
-        G4cout << "Wrong Angular type!!!!" << G4endl;
-        G4cout << "Wrong Angular type!!!!" << G4endl;
-        G4int iN;
-        G4cin >> iN;
-        return;
-    }
-}
-
 void GenePrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) {
     G4double thetacmr, thetalabr;
     thetacmr = hAng->GetRandom() / 180.0 * pi;

src/GenePrimaryGeneratorActionMessenger.cc

@@ -17,26 +17,15 @@ GenePrimaryGeneratorActionMessenger::GenePrimaryGeneratorActionMessenger(GenePri
     BeamEnergyCmd->SetUnitCategory("Energy");
     BeamEnergyCmd->SetRange("BeamEnergy>0.0");
     BeamEnergyCmd->AvailableForStates(G4State_PreInit, G4State_Idle);
-
-    AngularCmd = new G4UIcmdWithAString("/Gene/PrimaryGA/SetAngular", this);
-    AngularCmd->SetGuidance("Set type of angular distribution.");
-    AngularCmd->SetParameterName("Angular", false);
-    AngularCmd->AvailableForStates(G4State_PreInit);
 }
 
 GenePrimaryGeneratorActionMessenger::~GenePrimaryGeneratorActionMessenger() {
     delete BeamEnergyCmd;
     delete PrimaryDir;
-    delete AngularCmd;
-    delete ReactionCmd;
 }
 
 void GenePrimaryGeneratorActionMessenger::SetNewValue(G4UIcommand* command, G4String newValue) {
     if (command == BeamEnergyCmd) {
         pPrimaryGeneratorAction->SetBeamEnergy(BeamEnergyCmd->GetNewDoubleValue(newValue));
     }
-
-    if (command == AngularCmd) {
-        pPrimaryGeneratorAction->SetAngType(newValue);
-    }
 }