diff --git a/include/GeneAnalysisManager.hh b/include/GeneAnalysisManager.hh index 8b14909..3a39b4d 100644 --- a/include/GeneAnalysisManager.hh +++ b/include/GeneAnalysisManager.hh @@ -44,7 +44,6 @@ private: static GeneAnalysisManager *instance; struct detEvent mydata; - // G4int evtNo; G4String analysisFileName; TFile *theTFile; diff --git a/include/GeneDetectorConstructionMessenger.hh b/include/GeneDetectorConstructionMessenger.hh index 52536c1..85733f6 100644 --- a/include/GeneDetectorConstructionMessenger.hh +++ b/include/GeneDetectorConstructionMessenger.hh @@ -23,7 +23,6 @@ private: G4UIdirectory* DetectorDir; G4UIcmdWithADoubleAndUnit* ZoffsetCmd; G4UIcmdWithADoubleAndUnit* ZactiveshiftCmd; - G4UIcmdWithADouble* PE_BfractionCmd; G4UIcmdWithAString* DetReactionCmd; }; diff --git a/include/GeneHe3detSD.hh b/include/GeneHe3detSD.hh index 387bde4..cd9f900 100644 --- a/include/GeneHe3detSD.hh +++ b/include/GeneHe3detSD.hh @@ -22,7 +22,6 @@ public: private: GeneHe3detHitsCollection* scintillatorCollection; - G4int HitID; }; #endif diff --git a/include/GenePhysicsList.hh b/include/GenePhysicsList.hh deleted file mode 100644 index 398882f..0000000 --- a/include/GenePhysicsList.hh +++ /dev/null @@ -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 diff --git a/include/GenePrimaryGeneratorAction.hh b/include/GenePrimaryGeneratorAction.hh index 7e3e2c2..51037d6 100644 --- a/include/GenePrimaryGeneratorAction.hh +++ b/include/GenePrimaryGeneratorAction.hh @@ -16,7 +16,6 @@ public: void GeneratePrimaries(G4Event* anEvent); void SetBeamEnergy(G4double); - void SetAngType(G4String&); private: G4ParticleGun* particleGun; diff --git a/include/GenePrimaryGeneratorActionMessenger.hh b/include/GenePrimaryGeneratorActionMessenger.hh index bc9bfc9..6bd5073 100644 --- a/include/GenePrimaryGeneratorActionMessenger.hh +++ b/include/GenePrimaryGeneratorActionMessenger.hh @@ -21,8 +21,6 @@ private: G4UIdirectory* PrimaryDir; G4UIcmdWithADoubleAndUnit* BeamEnergyCmd; - G4UIcmdWithAString* AngularCmd; - G4UIcmdWithAString* ReactionCmd; }; #endif diff --git a/init_vis.mac b/init_vis.mac index 8e46032..2ed5342 100644 --- a/init_vis.mac +++ b/init_vis.mac @@ -1,16 +1,10 @@ -# Macro file for the initialization of example B1 -# in interactive session -# -# Set some default verbose -/control/verbose 2 -/control/saveHistory +/control/verbose 1 /run/verbose 2 -# -# Change the default number of threads (in multi-threaded mode) -#/run/numberOfThreads 4 -# +/event/verbose 0 +/tracking/verbose 0 + # Initialize kernel /run/initialize -# + # Visualization setting /control/execute vis.mac diff --git a/main.cpp b/main.cpp index 7023b6d..25cca0b 100644 --- a/main.cpp +++ b/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; } diff --git a/run.mac b/run.mac index b2c2ad4..0db4fd2 100644 --- a/run.mac +++ b/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 -# 2nd energy -# 3rd angular +# 1st energy +# 2nd offset +# 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 diff --git a/src/GeneAnalysisManager.cc b/src/GeneAnalysisManager.cc index 9f7095f..54b7416 100644 --- a/src/GeneAnalysisManager.cc +++ b/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 <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(); diff --git a/src/GeneHe3detSD.cc b/src/GeneHe3detSD.cc index 3dd784d..a2dcd5b 100644 --- a/src/GeneHe3detSD.cc +++ b/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]); - - HitID = -1; + static G4int HCID = -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(); diff --git a/src/GenePhysicsList.cc b/src/GenePhysicsList.cc deleted file mode 100644 index 1cabc3e..0000000 --- a/src/GenePhysicsList.cc +++ /dev/null @@ -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 - -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(); -} diff --git a/src/GenePrimaryGeneratorAction.cc b/src/GenePrimaryGeneratorAction.cc index 89889ad..0eedb67 100644 --- a/src/GenePrimaryGeneratorAction.cc +++ b/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; diff --git a/src/GenePrimaryGeneratorActionMessenger.cc b/src/GenePrimaryGeneratorActionMessenger.cc index 80f5f51..be34b08 100644 --- a/src/GenePrimaryGeneratorActionMessenger.cc +++ b/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); - } }