Skip to content
Snippets Groups Projects
Commit e14a5361 authored by Adrien Matta's avatar Adrien Matta :skull_crossbones:
Browse files

* Importing new physics list for neutron

parent 99f089b5
No related branches found
No related tags found
No related merge requests found
Pipeline #231982 passed
Stage: build-NPLib
Stage: build-NPSimulation
Stage: test
Hide whitespace changes
Inline Side-by-side
NPSimulation/Process/PhysicsList.cc
+ 173
− 117
View file @ e14a5361
......@@ -21,6 +21,8 @@
* *
*****************************************************************************/
#include <iostream>
#include <fstream>
// NPS
#include "PhysicsList.hh"
// NPL
......@@ -31,41 +33,52 @@
// Local physic directly implemented from the Hadronthrapy example
// Physic dedicated to the ion-ion inelastic processes
#include "Decay.hh"
#include "menate_R.hh"
#ifndef GEANT4_11
#include "NPIonIonInelasticPhysic.hh"
#endif
// G4
#include "G4BetheBlochIonGasModel.hh"
#include "G4BraggIonGasModel.hh"
#include "G4Material.hh"
#include "G4MaterialTable.hh"
#include "G4SystemOfUnits.hh"
#include "G4RunManager.hh"
#include "G4PhysListFactory.hh"
#include "G4VPhysicsConstructor.hh"
#include "G4LossTableManager.hh"
#include "G4UnitsTable.hh"
#include "G4ProcessManager.hh"
#include "G4ProcessVector.hh"
#include "G4FastSimulationManagerProcess.hh"
#include "G4StepLimiter.hh"
#include "G4BraggIonGasModel.hh"
#include "G4BetheBlochIonGasModel.hh"
#include "G4IonFluctuations.hh"
#include "G4IonParametrisedLossModel.hh"
#include "G4LevelManager.hh"
#include "G4LossTableManager.hh"
#include "G4UniversalFluctuation.hh"
#include "G4NuclearLevelData.hh"
#include "G4LevelManager.hh"
#include "G4PAIModel.hh"
#include "G4PAIPhotModel.hh"
#include "G4PhysListFactory.hh"
#include "G4ProcessManager.hh"
#include "G4RadioactiveDecay.hh"
#include "G4RunManager.hh"
#include "G4StepLimiter.hh"
#include "G4SystemOfUnits.hh"
#include "G4UImanager.hh"
#include "G4UnitsTable.hh"
#include "G4UniversalFluctuation.hh"
#include "G4VPhysicsConstructor.hh"
#include "menate_R.hh"
#include "G4ParticleTypes.hh"
#include "G4ParticleTable.hh"
#include "G4DecayPhysics.hh"
#include "G4RadioactiveDecayPhysics.hh"
#include "G4EmStandardPhysics.hh"
#include "G4EmExtraPhysics.hh"
#include "G4IonINCLXXPhysics.hh"
#include "G4StoppingPhysics.hh"
#include "G4HadronElasticPhysics.hh"
#include "G4HadronElasticPhysicsHP.hh"
#include "G4NeutronTrackingCut.hh"
#include "G4HadronPhysicsINCLXX.hh"
#ifdef USE_NEUTRONHP
#include "NeutronHPphysics.hh"
#endif
/////////////////////////////////////////////////////////////////////////////
PhysicsList::PhysicsList() : G4VUserPhysicsList() {
// UI commands to activate one or more High Energy Processes
// G4UImanager* UI = G4UImanager::GetUIpointer();
// UI->ApplyCommand("/physics_list/em/MuonNuclear true");
m_EmList = "Option4";
defaultCutValue = 1 * mm; // 0.2*mm;
opticalPhysicsList = NULL;
......@@ -89,101 +102,108 @@ PhysicsList::PhysicsList() : G4VUserPhysicsList() {
// Using the more accurate option4
emPhysicsList = NULL;
if (m_EmList == "Option4") {
emPhysicsList = new G4EmStandardPhysics_option4();
cout << "//// Using G4EmStandardPhysics_option4 Physics List ////" << endl;
if (m_EmList == "INCLXX_EM"){
cout << "Choosing to use INCLXX included EMList (probably standard)" << endl;
}
else if (m_EmList == "Option1") {
emPhysicsList = new G4EmStandardPhysics_option1();
else if(m_EmList== "Option1"){
emPhysicsList=new G4EmStandardPhysics_option1();
cout << "//// Using G4EmStandardPhysics_option1 Physics List ////" << endl;
}
else if (m_EmList == "Option2") {
emPhysicsList = new G4EmStandardPhysics_option2();
else if(m_EmList== "Option2"){
emPhysicsList=new G4EmStandardPhysics_option2();
cout << "//// Using G4EmStandardPhysics_option2 Physics List ////" << endl;
}
else if (m_EmList == "Option3") {
emPhysicsList = new G4EmStandardPhysics_option3();
else if(m_EmList== "Option3"){
emPhysicsList=new G4EmStandardPhysics_option3();
cout << "//// Using G4EmStandardPhysics_option3 Physics List ////" << endl;
}
else if (m_EmList == "Standard") {
emPhysicsList = new G4EmStandardPhysics();
else if(m_EmList == "Option4"){
emPhysicsList=new G4EmStandardPhysics_option4();
cout << "//// Using G4EmStandardPhysics_option4 Physics List ////" << endl;
}
else if(m_EmList== "Standard"){
emPhysicsList=new G4EmStandardPhysics();
cout << "//// Using G4EmStandardPhysics default EM constructor Physics List ////" << endl;
}
else if (m_EmList == "Livermore") {
emPhysicsList = new G4EmLivermorePhysics();
else if(m_EmList== "Livermore"){
emPhysicsList=new G4EmLivermorePhysics();
cout << "//// Using G4EmLivermorePhysics Physics List ////" << endl;
}
else if (m_EmList == "Penelope") {
emPhysicsList = new G4EmPenelopePhysics();
else if( m_EmList == "Penelope"){
emPhysicsList=new G4EmPenelopePhysics();
cout << "//// Using G4EmPenelopePhysics Physics List ////" << endl;
}
else {
std::cout << "\r\033[1;31mERROR: User given physics list " << m_EmList
<< " is not supported, option are Option4 Livermore Penelope\033[0m" << std::endl;
else{
std::cout << "\r\033[1;31mERROR: User given physics list " << m_EmList << " is not supported, option are Option4 Livermore Penelope\033[0m" <<std::endl;
exit(1);
}
emConfig = G4LossTableManager::Instance()->EmConfigurator();
RegisterPhysics(emPhysicsList);
// Hadronic physics
if (m_IonBinaryCascadePhysics) {
m_PhysList["IonBinaryCascadePhysics"] = new G4IonBinaryCascadePhysics();
if(m_IonBinaryCascadePhysics){
RegisterPhysics(new G4IonBinaryCascadePhysics());
cout << "//// Using G4IonBinaryCascadePhysics Physics List ////" << endl;
}
if (m_EmExtraPhysics)
m_PhysList["EmExtraPhysics"] = new G4EmExtraPhysics();
if (m_HadronElasticPhysics) {
m_PhysList["G4HadronElasticPhysics"] = new G4HadronElasticPhysics();
m_PhysList["G4IonElasticPhysics"] = new G4IonElasticPhysics();
if(m_EmExtraPhysics){
RegisterPhysics(new G4EmExtraPhysics());
}
if(m_HadronElasticPhysics){
RegisterPhysics(new G4HadronElasticPhysics());
RegisterPhysics(new G4IonElasticPhysics());
cout << "//// Using G4HadronElasticPhysics Physics List ////" << endl;
cout << "//// Using G4IonElasticPhysics Physics List ////" << endl;
}
#ifndef GEANT4_11
if (m_NPIonInelasticPhysics) {
m_PhysList["NPIonInelasticPhysics"] = new NPIonIonInelasticPhysic();
if(m_NPIonInelasticPhysics){
RegisterPhysics(new NPIonIonInelasticPhysic());
cout << "//// Using NPIonIonInelasticPhysic Physics List ////" << endl;
}
#endif
#ifdef USE_NEUTRONHP
if (m_NeutronHP) {
m_PhysList["NeutronHPphysics"] = new NeutronHPphysics("neutronHP");
if(m_NeutronHP){
RegisterPhysics(new NeutronHPphysics("neutronHP"));
cout << "//// Using NeutronHPPhysics ////" << endl;
}
#endif
if (m_StoppingPhysics) {
m_PhysList["StoppingPhysics"] = new G4StoppingPhysics();
if(m_StoppingPhysics){
RegisterPhysics(new G4StoppingPhysics());
cout << "//// Using G4StoppingPhysics Physics List ////" << endl;
}
if (m_HadronPhysicsINCLXX) {
m_PhysList["HadronPhysicsINCLXX"] = new G4HadronPhysicsINCLXX();
m_PhysList["IonPhysicsINCLXX"] = new G4IonINCLXXPhysics();
if(m_HadronPhysicsINCLXX){
RegisterPhysics(new G4HadronPhysicsINCLXX());
RegisterPhysics(new G4IonINCLXXPhysics());
cout << "//// Using INCLXX Physics List ////" << endl;
}
if (m_HadronPhysicsQGSP_BIC_HP) {
if(m_HadronPhysicsQGSP_BERT_HP){
RegisterPhysics(new G4HadronPhysicsQGSP_BERT_HP());
cout << "//// Using QGSP_BERT_HP Physics List ////" << endl;
}
if(m_HadronPhysicsQGSP_BIC_HP){
#if NPS_GEANT4_VERSION_MAJOR > 9
m_PhysList["HadronPhysicsQGSP_BIC_HP"] = new G4HadronPhysicsQGSP_BIC_HP();
RegisterPhysics(new G4HadronPhysicsQGSP_BIC_HP());
cout << "//// Using QGSP_BIC_HP Physics List ////" << endl;
#else
std::cout << "\r\032[1;31m Warning: physics list HadronPhysicsQGSP_BIC_HP require Geant4 10, process not activated "
"\033[0m"
<< std::endl;
std::cout << "\r\032[1;31m Warning: physics list HadronPhysicsQGSP_BIC_HP requires Geant4 10, process not activated \033[0m" <<std::endl;
#endif
}
if (m_HadronPhysicsQGSP_BERT_HP) {
m_PhysList["HadronPhysicsQGSP_BERT_HP"] = new G4HadronPhysicsQGSP_BERT_HP();
cout << "//// Using QGSP_BERT_HP Physics List ////" << endl;
if(m_HadronPhysicsQGSP_INCLXX_HP){
#if NPS_GEANT4_VERSION_MAJOR > 9
AddPackage("QGSP_INCLXX_HP");
#else
std::cout << "\r\032[1;31m Warning: physics list HadronPhysicsQGSP_INCLXX_HP requires Geant4 10, process not activated \033[0m" <<std::endl;
#endif
}
// Optical Photon for scintillator simulation
if (m_OpticalPhysics) {
opticalPhysicsList = new G4OpticalPhysics(0);
RegisterPhysics(new G4OpticalPhysics(0));
/*opticalPhysicsList->SetMaxNumPhotonsPerStep(100);
opticalPhysicsList->SetScintillationYieldFactor(0.1);
opticalPhysicsList->SetTrackSecondariesFirst(kScintillation,true);
......@@ -193,21 +213,21 @@ PhysicsList::PhysicsList() : G4VUserPhysicsList() {
}
// Drift electron for gazeous detector simulation
if (m_DriftElectronPhysics) {
driftElectronPhysicsList = new G4DriftElectronPhysics(0);
if(m_DriftElectronPhysics){
RegisterPhysics(new G4DriftElectronPhysics(0));
driftElectronPhysicsList->SetMaxNumDriftElectronPerStep(1e6);
cout << "//// Using Drift Electron Physics List ////" << endl;
}
if (m_IonGasModels) {
if(m_IonGasModels){
AddIonGasModels();
cout << "//// Using Ion Gas Model Physics List ////" << endl;
}
if (m_pai) {
if(m_pai){
AddPAIModel("pai");
cout << "//// Using PAI Model Physics List ////" << endl;
}
if (m_pai_photon) {
if(m_pai_photon){
AddPAIModel("pai_photon");
cout << "//// Using PAI PHOTON Model Physics List ////" << endl;
}
......@@ -216,16 +236,15 @@ PhysicsList::PhysicsList() : G4VUserPhysicsList() {
// Add Radioactive decay
// Gamma decay of known states
if (m_Decay) {
decay_List = new G4DecayPhysics();
radioactiveDecay_List = new G4RadioactiveDecayPhysics();
m_PhysList["decay_list"] = decay_List;
m_PhysList["radioactiveDecay"] = radioactiveDecay_List;
RegisterPhysics(new G4DecayPhysics());
RegisterPhysics(new G4RadioactiveDecayPhysics());
}
else {
decay_List = 0;
radioactiveDecay_List = 0;
}
emConfig = G4LossTableManager::Instance()->EmConfigurator();
}
////////////////////////////////////////////////////////////////////////////////
......@@ -240,6 +259,11 @@ void PhysicsList::ReadConfiguration(std::string filename) {
m_HadronPhysicsQGSP_BIC_HP = 0;
m_HadronPhysicsQGSP_BERT_HP = 0;
m_HadronPhysicsINCLXX = 0;
m_INCLXXPhysics = 0;
m_HadronPhysicsQGSP_INCLXX_HP = 0;
m_HadronPhysicsQGSP_INCLXX = 0;
m_HadronPhysicsFTFP_INCLXX_HP = 0;
m_HadronPhysicsFTFP_INCLXX = 0;
m_Decay = 0;
m_IonGasModels = 0;
m_pai = 0;
......@@ -286,6 +310,26 @@ void PhysicsList::ReadConfiguration(std::string filename) {
m_HadronPhysicsQGSP_BERT_HP = value;
else if (name == "HadronPhysicsINCLXX")
m_HadronPhysicsINCLXX = value;
else if (name == "HadronPhysicsQGSP_INCLXX_HP"){
m_HadronPhysicsQGSP_INCLXX_HP= value;
if(value)
m_INCLXXPhysics= true;
}
else if (name == "HadronPhysicsQGSP_INCLXX"){
m_HadronPhysicsQGSP_INCLXX= value;
if(value)
m_INCLXXPhysics= true;
}
else if (name == "HadronPhysicsFTFP_INCLXX_HP"){
m_HadronPhysicsFTFP_INCLXX_HP= value;
if(value)
m_INCLXXPhysics= true;
}
else if (name == "HadronPhysicsFTFP_INCLXX"){
m_HadronPhysicsFTFP_INCLXX= value;
if(value)
m_INCLXXPhysics= true;
}
else if (name == "Decay")
m_Decay = value;
else if (name == "IonGasModels")
......@@ -313,6 +357,22 @@ void PhysicsList::ReadConfiguration(std::string filename) {
/////////////////////////////////////////////////////////////////////////////
PhysicsList::~PhysicsList() {}
/////////////////////////////////////////////////////////////////////////////
void PhysicsList::AddPackage(const G4String& name)
{
G4PhysListFactory factory;
G4VModularPhysicsList* phys =factory.GetReferencePhysList(name);
G4int i=0;
const G4VPhysicsConstructor* elem= phys->GetPhysics(i);
G4VPhysicsConstructor* tmp = const_cast<G4VPhysicsConstructor*> (elem);
while (elem !=0)
{
RegisterPhysics(tmp);
elem= phys->GetPhysics(++i) ;
tmp = const_cast<G4VPhysicsConstructor*> (elem);
}
}
/////////////////////////////////////////////////////////////////////////////
void PhysicsList::AddIonGasModels() {
G4ParticleTable::G4PTblDicIterator* particleIterator = G4ParticleTable::GetParticleTable()->GetIterator();
......@@ -366,20 +426,10 @@ void PhysicsList::NewPAIModel(const G4ParticleDefinition* part, const G4String&
/////////////////////////////////////////////////////////////////////////////
void PhysicsList::ConstructParticle() {
if (m_OpticalPhysics) {
// G4OpticalPhoton::OpticalPhotonDefinition();
((G4VPhysicsConstructor*)opticalPhysicsList)->ConstructParticle();
}
if (m_DriftElectronPhysics) {
((G4VPhysicsConstructor*)driftElectronPhysicsList)->ConstructParticle();
if(m_INCLXXPhysics){
this->G4VModularPhysicsList::ConstructParticle();
}
if (decay_List) {
decay_List->ConstructParticle();
radioactiveDecay_List->ConstructParticle();
}
else {
// If decay is not activated we have to declare the particle ourself
G4He3::He3Definition();
......@@ -418,41 +468,47 @@ void PhysicsList::ConstructParticle() {
}
/////////////////////////////////////////////////////////////////////////////
void PhysicsList::ConstructProcess() {
// Transportation
AddTransportation();
// Electromagnetic physics
emPhysicsList->ConstructProcess();
if (opticalPhysicsList) {
((G4VPhysicsConstructor*)opticalPhysicsList)->ConstructProcess();
}
if (driftElectronPhysicsList) {
((G4VPhysicsConstructor*)driftElectronPhysicsList)->ConstructProcess();
}
// Hadronic physics
std::map<std::string, G4VPhysicsConstructor*>::iterator it;
for (it = m_PhysList.begin(); it != m_PhysList.end(); it++) {
it->second->ConstructProcess();
}
BiasCrossSectionByFactor(m_IonBinaryCascadePhysics);
em_parameters = G4EmParameters::Instance();
em_parameters->SetFluo(true);
em_parameters->SetAuger(true);
if(m_INCLXXPhysics){
//G4VModularPhysicsList::ReplacePhysics(new G4EmStandardPhysics_option4());
this->G4VModularPhysicsList::ConstructProcess();
}
else{
// Transportation
AddTransportation();
// Electromagnetic physics
emPhysicsList->ConstructProcess();
if (opticalPhysicsList) {
((G4VPhysicsConstructor*)opticalPhysicsList)->ConstructProcess();
}
if (driftElectronPhysicsList) {
((G4VPhysicsConstructor*)driftElectronPhysicsList)->ConstructProcess();
}
// Hadronic physics
for(auto it = m_PhysList.begin(), end = m_PhysList.end(); it!= end; it++){
it->second->ConstructProcess();
}
BiasCrossSectionByFactor(m_IonBinaryCascadePhysics);
em_parameters = G4EmParameters::Instance();
em_parameters->SetFluo(true);
em_parameters->SetAuger(true);
#if NPS_GEANT4_VERSION_MAJOR > 10
#if NPS_GEANT4_VERSION_MINOR > 2
em_parameters->SetDeexActiveRegion("DefaultRegionForTheWorld", true, true, true);
em_parameters->SetDeexcitationIgnoreCut(true);
em_parameters->SetDeexActiveRegion("DefaultRegionForTheWorld", true, true, true);
em_parameters->SetDeexcitationIgnoreCut(true);
#endif
#endif
AddParametrisation();
AddLevelData();
return;
}
AddParametrisation();
AddLevelData();
return;
}
/////////////////////////////////////////////////////////////////////////////
void PhysicsList::AddStepMax() {}
/////////////////////////////////////////////////////////////////////////////
......
NPSimulation/Process/PhysicsList.hh
+ 13
− 6
View file @ e14a5361
/*****************************************************************************
* Copyright (C) 2009-2016 this file is part of the NPTool Project *
* Copyright (C) 2009-2023 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
......@@ -9,7 +9,7 @@
* Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk *
* *
* Creation Date : January 2009 *
* Last update : October 2015 *
* Last update : October 2023 *
*---------------------------------------------------------------------------*
* Decription: *
* Modular Physics list calling Geant4 reference list *
......@@ -25,6 +25,7 @@
#define PhysicsList_h 1
#include "G4VUserPhysicsList.hh"
#include "G4VModularPhysicsList.hh"
#include "G4EmConfigurator.hh"
#include "globals.hh"
......@@ -47,6 +48,7 @@
#include "G4EmStandardPhysics_option2.hh"
#include "G4EmStandardPhysics_option3.hh"
#include "G4EmStandardPhysics_option4.hh"
#include "G4EmStandardPhysics.hh"
#include "G4EmExtraPhysics.hh"
#include "G4StoppingPhysics.hh"
......@@ -57,7 +59,6 @@
#include "G4DriftElectronPhysics.hh"
//Hadronique
#include "G4IonElasticPhysics.hh"
#include "G4HadronElasticPhysics.hh"
#include "G4HadronInelasticQBBC.hh"
......@@ -67,13 +68,14 @@
#if NPS_GEANT4_VERSION_MAJOR > 9
#include "G4HadronPhysicsQGSP_BIC_HP.hh"
#include "G4HadronPhysicsQGSP_BERT_HP.hh"
#endif
#include "G4HadronPhysicsQGSP_BIC.hh"
#include "INCLXXPhysicsListHelper.hh"
#endif
#include "G4HadronPhysicsQGSP_BERT.hh"
class G4VPhysicsConstructor;
class PhysicsList: public G4VUserPhysicsList{
class PhysicsList: public G4VModularPhysicsList{
public:
PhysicsList();
virtual ~PhysicsList();
......@@ -84,11 +86,11 @@ class PhysicsList: public G4VUserPhysicsList{
void ConstructProcess();
void AddStepMax();
void AddParametrisation();
void AddPackage(const G4String& name);
void BiasCrossSectionByFactor(double factor);
void AddIonGasModels();
void AddPAIModel(const G4String& modname);
void NewPAIModel(const G4ParticleDefinition* part, const G4String& modname,const G4String& procname);
void AddPackage(const G4String& name);
void AddLevelData();
private:
std::map<std::string,G4VPhysicsConstructor*> m_PhysList;
......@@ -115,6 +117,11 @@ class PhysicsList: public G4VUserPhysicsList{
double m_HadronPhysicsQGSP_BIC_HP;
double m_HadronPhysicsQGSP_BERT_HP;
double m_HadronPhysicsINCLXX;
bool m_INCLXXPhysics;
double m_HadronPhysicsQGSP_INCLXX_HP;
double m_HadronPhysicsQGSP_INCLXX;
double m_HadronPhysicsFTFP_INCLXX_HP;
double m_HadronPhysicsFTFP_INCLXX;
double m_Decay;
double m_IonGasModels;
double m_pai;
......
NPSimulation/ressources/input/PhysicsListOption.txt.in
+ 4
− 0
View file @ e14a5361
......@@ -6,6 +6,10 @@ EmExtraPhysics 0
HadronElasticPhysics 0
StoppingPhysics 0
OpticalPhysics 0
HadronElasticPhysics 0
HadronPhysicsQGSP_BIC_HP 0
HadronPhysicsQGSP_BERT_HP 0
HadronPhysicsQGSP_INCLXX_HP 0
HadronPhysicsINCLXX 0
HadronPhysicsQGSP_BIC_HP 0
Decay 0
Projects/Minos/Minos.detector
+ 3
− 1
View file @ e14a5361
......@@ -19,4 +19,6 @@ Minos
BaseLine= 250
Sampling= 10
ZOffset= 0 mm
XML= /Users/matta/mrdc/xml/s034/MINOS.xml
TargetZOffset= 0 mm
ZRotation= 0 mm
XML= ../S034/db/MINOS.xml
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment