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

* Adding custom level decay option 

        - User can add a level decay file in geant4 format
        - this is done in the event generator file
        - the level scheme file should be saved in the tree -> to be done
parent d05242f3
No related branches found
No related tags found
No related merge requests found
Pipeline #46581 passed
Stage: build-NPLib
Stage: build-NPSimulation
Stage: test
Expand all Hide whitespace changes
Inline Side-by-side
NPSimulation/Process/PhysicsList.cc
+ 407
− 371
View file @ ad528b71
...@@ -23,6 +23,10 @@ ...@@ -23,6 +23,10 @@
//NPS //NPS
#include "PhysicsList.hh" #include "PhysicsList.hh"
//NPL
#include "NPInputParser.h"
#include "NPOptionManager.h"
#include "NPNucleus.h"
// Local physic directly implemented from the Hadronthrapy example // Local physic directly implemented from the Hadronthrapy example
// Physic dedicated to the ion-ion inelastic processes // Physic dedicated to the ion-ion inelastic processes
...@@ -44,7 +48,8 @@ ...@@ -44,7 +48,8 @@
#include "G4IonParametrisedLossModel.hh" #include "G4IonParametrisedLossModel.hh"
#include "G4UniversalFluctuation.hh" #include "G4UniversalFluctuation.hh"
#include "G4UImanager.hh" #include "G4UImanager.hh"
#include "G4NuclearLevelData.hh"
#include "G4LevelManager.hh"
#include "G4PAIModel.hh" #include "G4PAIModel.hh"
#include "G4PAIPhotModel.hh" #include "G4PAIPhotModel.hh"
#include "menate_R.hh" #include "menate_R.hh"
...@@ -55,227 +60,227 @@ PhysicsList::PhysicsList() : G4VUserPhysicsList(){ ...@@ -55,227 +60,227 @@ PhysicsList::PhysicsList() : G4VUserPhysicsList(){
// UI commands to activate one or more High Energy Processes // UI commands to activate one or more High Energy Processes
// G4UImanager* UI = G4UImanager::GetUIpointer(); // G4UImanager* UI = G4UImanager::GetUIpointer();
// UI->ApplyCommand("/physics_list/em/MuonNuclear true"); // UI->ApplyCommand("/physics_list/em/MuonNuclear true");
m_EmList = "Option4"; m_EmList = "Option4";
defaultCutValue = 1*mm;//0.2*mm; defaultCutValue = 1*mm;//0.2*mm;
opticalPhysicsList = NULL; opticalPhysicsList = NULL;
driftElectronPhysicsList = NULL; driftElectronPhysicsList = NULL;
ReadConfiguration("PhysicsListOption.txt"); ReadConfiguration("PhysicsListOption.txt");
G4LossTableManager::Instance(); G4LossTableManager::Instance();
SetVerboseLevel(0); SetVerboseLevel(0);
// ****** Definition of defaults for the physics processes ***** // ****** Definition of defaults for the physics processes *****
// ****** in case no physics is called by the macro file ***** // ****** in case no physics is called by the macro file *****
// //
// The default physics corresponds to the actual QGSP_BIC_HP list // The default physics corresponds to the actual QGSP_BIC_HP list
// but with the following differences: // but with the following differences:
// --> G4EmStandardPhysics_option4 for the electromagnetic processes // --> G4EmStandardPhysics_option4 for the electromagnetic processes
// is used n place of the less accurate G4EmStandardPhysics // is used n place of the less accurate G4EmStandardPhysics
// --> The G4RadioactiveDecayPhysics is add // --> The G4RadioactiveDecayPhysics is add
// --> G4HadronPhysicsQGSP_BIC is used in place of G4HadronPhysicsQGSP_BIC_HP // --> G4HadronPhysicsQGSP_BIC is used in place of G4HadronPhysicsQGSP_BIC_HP
// --> G4HadronElasticPhysics is used in place of G4HadronElasticPhysics_HP // --> G4HadronElasticPhysics is used in place of G4HadronElasticPhysics_HP
// Elecromagnetic physics // Elecromagnetic physics
// Using the more accurate option4 // Using the more accurate option4
emPhysicsList=NULL; emPhysicsList=NULL;
if(m_EmList == "Option4"){ if(m_EmList == "Option4"){
emPhysicsList = new G4EmStandardPhysics_option4(); emPhysicsList = new G4EmStandardPhysics_option4();
cout << "//// Using G4EmStandardPhysics_option4 Physics List ////" << endl; cout << "//// Using G4EmStandardPhysics_option4 Physics List ////" << endl;
} }
else if(m_EmList== "Option1"){ else if(m_EmList== "Option1"){
emPhysicsList = new G4EmStandardPhysics_option1(); emPhysicsList = new G4EmStandardPhysics_option1();
cout << "//// Using G4EmStandardPhysics_option1 Physics List ////" << endl; cout << "//// Using G4EmStandardPhysics_option1 Physics List ////" << endl;
} }
else if(m_EmList== "Option2"){ else if(m_EmList== "Option2"){
emPhysicsList = new G4EmStandardPhysics_option2(); emPhysicsList = new G4EmStandardPhysics_option2();
cout << "//// Using G4EmStandardPhysics_option2 Physics List ////" << endl; cout << "//// Using G4EmStandardPhysics_option2 Physics List ////" << endl;
} }
else if(m_EmList== "Option3"){ else if(m_EmList== "Option3"){
emPhysicsList = new G4EmStandardPhysics_option3(); emPhysicsList = new G4EmStandardPhysics_option3();
cout << "//// Using G4EmStandardPhysics_option3 Physics List ////" << endl; cout << "//// Using G4EmStandardPhysics_option3 Physics List ////" << endl;
} }
else if(m_EmList== "Standard"){ else if(m_EmList== "Standard"){
emPhysicsList = new G4EmStandardPhysics(); emPhysicsList = new G4EmStandardPhysics();
cout << "//// Using G4EmStandardPhysics default EM constructor Physics List ////" << endl; cout << "//// Using G4EmStandardPhysics default EM constructor Physics List ////" << endl;
} }
else if(m_EmList== "Livermore"){ else if(m_EmList== "Livermore"){
emPhysicsList = new G4EmLivermorePhysics(); emPhysicsList = new G4EmLivermorePhysics();
cout << "//// Using G4EmLivermorePhysics Physics List ////" << endl; cout << "//// Using G4EmLivermorePhysics Physics List ////" << endl;
} }
else if( m_EmList == "Penelope"){ else if( m_EmList == "Penelope"){
emPhysicsList = new G4EmPenelopePhysics(); emPhysicsList = new G4EmPenelopePhysics();
cout << "//// Using G4EmPenelopePhysics Physics List ////" << endl; cout << "//// Using G4EmPenelopePhysics Physics List ////" << endl;
} }
else{ 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; 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); exit(1);
} }
emConfig = G4LossTableManager::Instance()->EmConfigurator(); emConfig = G4LossTableManager::Instance()->EmConfigurator();
// Hadronic physics // Hadronic physics
if(m_IonBinaryCascadePhysics){ if(m_IonBinaryCascadePhysics){
m_PhysList["IonBinaryCascadePhysics"] = new G4IonBinaryCascadePhysics(); m_PhysList["IonBinaryCascadePhysics"] = new G4IonBinaryCascadePhysics();
cout << "//// Using G4IonBinaryCascadePhysics Physics List ////" << endl; cout << "//// Using G4IonBinaryCascadePhysics Physics List ////" << endl;
} }
if(m_EmExtraPhysics) if(m_EmExtraPhysics)
m_PhysList["EmExtraPhysics"]=new G4EmExtraPhysics(); m_PhysList["EmExtraPhysics"]=new G4EmExtraPhysics();
if(m_HadronElasticPhysics){ if(m_HadronElasticPhysics){
m_PhysList["G4HadronElasticPhysics"]=new G4HadronElasticPhysics(); m_PhysList["G4HadronElasticPhysics"]=new G4HadronElasticPhysics();
m_PhysList["G4IonElasticPhysics"]=new G4IonElasticPhysics(); m_PhysList["G4IonElasticPhysics"]=new G4IonElasticPhysics();
cout << "//// Using G4HadronElasticPhysics Physics List ////" << endl; cout << "//// Using G4HadronElasticPhysics Physics List ////" << endl;
cout << "//// Using G4IonElasticPhysics Physics List ////" << endl; cout << "//// Using G4IonElasticPhysics Physics List ////" << endl;
} }
if(m_NPIonInelasticPhysics){ if(m_NPIonInelasticPhysics){
m_PhysList["NPIonInelasticPhysics"]=new NPIonIonInelasticPhysic(); m_PhysList["NPIonInelasticPhysics"]=new NPIonIonInelasticPhysic();
cout << "//// Using NPIonIonInelasticPhysic Physics List ////" << endl; cout << "//// Using NPIonIonInelasticPhysic Physics List ////" << endl;
} }
if(m_StoppingPhysics){ if(m_StoppingPhysics){
m_PhysList["StoppingPhysics"]=new G4StoppingPhysics(); m_PhysList["StoppingPhysics"]=new G4StoppingPhysics();
cout << "//// Using G4StoppingPhysics Physics List ////" << endl; cout << "//// Using G4StoppingPhysics Physics List ////" << endl;
} }
if(m_HadronPhysicsINCLXX){ if(m_HadronPhysicsINCLXX){
m_PhysList["HadronPhysicsINCLXX"] = new G4HadronPhysicsINCLXX(); m_PhysList["HadronPhysicsINCLXX"] = new G4HadronPhysicsINCLXX();
m_PhysList["IonPhysicsINCLXX"] = new G4IonINCLXXPhysics(); m_PhysList["IonPhysicsINCLXX"] = new G4IonINCLXXPhysics();
cout << "//// Using INCLXX Physics List ////" << endl; cout << "//// Using INCLXX Physics List ////" << endl;
} }
if(m_HadronPhysicsQGSP_BIC_HP){ if(m_HadronPhysicsQGSP_BIC_HP){
#if NPS_GEANT4_VERSION_MAJOR > 9 #if NPS_GEANT4_VERSION_MAJOR > 9
m_PhysList["HadronPhysicsQGSP_BIC_HP"] = new G4HadronPhysicsQGSP_BIC_HP(); m_PhysList["HadronPhysicsQGSP_BIC_HP"] = new G4HadronPhysicsQGSP_BIC_HP();
cout << "//// Using QGSP_BIC_HP Physics List ////" << endl; cout << "//// Using QGSP_BIC_HP Physics List ////" << endl;
#else #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 require Geant4 10, process not activated \033[0m" <<std::endl;
#endif #endif
} }
// Optical Photon for scintillator simulation // Optical Photon for scintillator simulation
if(m_OpticalPhysics){ if(m_OpticalPhysics){
opticalPhysicsList = new G4OpticalPhysics(0); opticalPhysicsList = new G4OpticalPhysics(0);
opticalPhysicsList->SetMaxNumPhotonsPerStep(100); opticalPhysicsList->SetMaxNumPhotonsPerStep(100);
opticalPhysicsList->SetScintillationYieldFactor(0.1); opticalPhysicsList->SetScintillationYieldFactor(0.1);
opticalPhysicsList->SetTrackSecondariesFirst(kScintillation,true); opticalPhysicsList->SetTrackSecondariesFirst(kScintillation,true);
opticalPhysicsList->SetTrackSecondariesFirst(kCerenkov,true); opticalPhysicsList->SetTrackSecondariesFirst(kCerenkov,true);
cout << "//// Using Optical Photon Physics List ////" << endl; cout << "//// Using Optical Photon Physics List ////" << endl;
} }
// Drift electron for gazeous detector simulation // Drift electron for gazeous detector simulation
if(m_DriftElectronPhysics){ if(m_DriftElectronPhysics){
driftElectronPhysicsList = new G4DriftElectronPhysics(0); driftElectronPhysicsList = new G4DriftElectronPhysics(0);
driftElectronPhysicsList->SetMaxNumDriftElectronPerStep(1e6); driftElectronPhysicsList->SetMaxNumDriftElectronPerStep(1e6);
cout << "//// Using Drift Electron Physics List ////" << endl; cout << "//// Using Drift Electron Physics List ////" << endl;
} }
if(m_IonGasModels){ if(m_IonGasModels){
AddIonGasModels(); AddIonGasModels();
cout << "//// Using Ion Gas Model Physics List ////" << endl; cout << "//// Using Ion Gas Model Physics List ////" << endl;
} }
if(m_pai){ if(m_pai){
AddPAIModel("pai"); AddPAIModel("pai");
cout << "//// Using PAI Model Physics List ////" << endl; cout << "//// Using PAI Model Physics List ////" << endl;
} }
if(m_pai_photon){ if(m_pai_photon){
AddPAIModel("pai_photon"); AddPAIModel("pai_photon");
cout << "//// Using PAI PHOTON Model Physics List ////" << endl; cout << "//// Using PAI PHOTON Model Physics List ////" << endl;
} }
// Decay physics // Decay physics
// Add Radioactive decay // Add Radioactive decay
// Gamma decay of known states // Gamma decay of known states
if(m_Decay){ if(m_Decay){
decay_List = new G4DecayPhysics(); decay_List = new G4DecayPhysics();
radioactiveDecay_List = new G4RadioactiveDecayPhysics() ; radioactiveDecay_List = new G4RadioactiveDecayPhysics() ;
m_PhysList["decay_list"]= decay_List; m_PhysList["decay_list"]= decay_List;
m_PhysList["radioactiveDecay"]=radioactiveDecay_List; m_PhysList["radioactiveDecay"]=radioactiveDecay_List;
} }
else{ else{
decay_List = 0; decay_List = 0;
radioactiveDecay_List = 0; radioactiveDecay_List = 0;
} }
} }
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
void PhysicsList::ReadConfiguration(std::string filename){ void PhysicsList::ReadConfiguration(std::string filename){
m_IonBinaryCascadePhysics = 0; m_IonBinaryCascadePhysics = 0;
m_EmExtraPhysics = 0; m_EmExtraPhysics = 0;
m_HadronElasticPhysics = 0; m_HadronElasticPhysics = 0;
m_NPIonInelasticPhysics = 0; m_NPIonInelasticPhysics = 0;
m_StoppingPhysics = 0; m_StoppingPhysics = 0;
m_OpticalPhysics = 0; m_OpticalPhysics = 0;
m_DriftElectronPhysics = 0; m_DriftElectronPhysics = 0;
m_HadronPhysicsQGSP_BIC_HP = 0; m_HadronPhysicsQGSP_BIC_HP = 0;
m_HadronPhysicsINCLXX = 0; m_HadronPhysicsINCLXX = 0;
m_Decay = 0; m_Decay = 0;
m_IonGasModels = 0; m_IonGasModels = 0;
m_pai= 0; m_pai= 0;
m_pai_photon= 0; m_pai_photon= 0;
m_Menate_R = 0; m_Menate_R = 0;
std::ifstream file(filename.c_str()); std::ifstream file(filename.c_str());
if(!file.is_open()){ if(!file.is_open()){
std::cout << "WARNING: Could not find Physics List option file " << filename << " | Using default Physics List" << std::endl; std::cout << "WARNING: Could not find Physics List option file " << filename << " | Using default Physics List" << std::endl;
return; return;
} }
std::cout << "Reading Physics list option file " << filename << std::endl; std::cout << "Reading Physics list option file " << filename << std::endl;
std::string name; std::string name;
std::string str_value; std::string str_value;
double value; double value;
while(file >> name >> str_value){ while(file >> name >> str_value){
value = std::atof(str_value.c_str()); value = std::atof(str_value.c_str());
if(name == "EmPhysicsList") if(name == "EmPhysicsList")
m_EmList = str_value; m_EmList = str_value;
else if(name == "DefaultCutOff") else if(name == "DefaultCutOff")
defaultCutValue = value*mm; defaultCutValue = value*mm;
else if(name == "IonBinaryCascadePhysics") else if(name == "IonBinaryCascadePhysics")
m_IonBinaryCascadePhysics = value; m_IonBinaryCascadePhysics = value;
else if(name == "NPIonInelasticPhysics") else if(name == "NPIonInelasticPhysics")
m_NPIonInelasticPhysics = value; m_NPIonInelasticPhysics = value;
else if (name == "EmExtraPhysics") else if (name == "EmExtraPhysics")
m_EmExtraPhysics= value; m_EmExtraPhysics= value;
else if (name == "HadronElasticPhysics") else if (name == "HadronElasticPhysics")
m_HadronElasticPhysics= value; m_HadronElasticPhysics= value;
else if (name == "StoppingPhysics") else if (name == "StoppingPhysics")
m_StoppingPhysics= value; m_StoppingPhysics= value;
else if (name == "OpticalPhysics") else if (name == "OpticalPhysics")
m_OpticalPhysics= value; m_OpticalPhysics= value;
else if (name == "DriftElectronPhysics") else if (name == "DriftElectronPhysics")
m_DriftElectronPhysics= value; m_DriftElectronPhysics= value;
else if (name == "HadronPhysicsQGSP_BIC_HP") else if (name == "HadronPhysicsQGSP_BIC_HP")
m_HadronPhysicsQGSP_BIC_HP= value; m_HadronPhysicsQGSP_BIC_HP= value;
else if (name == "HadronPhysicsINCLXX") else if (name == "HadronPhysicsINCLXX")
m_HadronPhysicsINCLXX= value; m_HadronPhysicsINCLXX= value;
else if (name == "Decay") else if (name == "Decay")
m_Decay = value; m_Decay = value;
else if (name == "IonGasModels") else if (name == "IonGasModels")
m_IonGasModels= value; m_IonGasModels= value;
else if (name == "pai") else if (name == "pai")
m_pai = value; m_pai = value;
else if (name == "pai_photon") else if (name == "pai_photon")
m_pai_photon = value; m_pai_photon = value;
else if (name == "Menate_R") else if (name == "Menate_R")
m_Menate_R = value; m_Menate_R = value;
else else
std::cout <<"WARNING: Physics List Token '" << name << "' unknown. Token is ignored." << std::endl; std::cout <<"WARNING: Physics List Token '" << name << "' unknown. Token is ignored." << std::endl;
} }
// Most special process need decay to be activated // Most special process need decay to be activated
if( (m_IonBinaryCascadePhysics || m_EmExtraPhysics || m_HadronElasticPhysics || m_NPIonInelasticPhysics if( (m_IonBinaryCascadePhysics || m_EmExtraPhysics || m_HadronElasticPhysics || m_NPIonInelasticPhysics
|| m_StoppingPhysics || m_HadronPhysicsQGSP_BIC_HP || m_HadronPhysicsINCLXX) && !m_Decay){ || m_StoppingPhysics || m_HadronPhysicsQGSP_BIC_HP || m_HadronPhysicsINCLXX) && !m_Decay){
m_Decay = true; m_Decay = true;
std::cout << "Information: Selected process require Decay to be activated." << std::endl; std::cout << "Information: Selected process require Decay to be activated." << std::endl;
} }
} }
///////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////
PhysicsList::~PhysicsList(){ PhysicsList::~PhysicsList(){
...@@ -284,199 +289,230 @@ PhysicsList::~PhysicsList(){ ...@@ -284,199 +289,230 @@ PhysicsList::~PhysicsList(){
///////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////
void PhysicsList::AddIonGasModels() void PhysicsList::AddIonGasModels()
{ {
G4ParticleTable::G4PTblDicIterator* particleIterator = G4ParticleTable::GetParticleTable()->GetIterator(); G4ParticleTable::G4PTblDicIterator* particleIterator = G4ParticleTable::GetParticleTable()->GetIterator();
particleIterator->reset(); particleIterator->reset();
while ((*particleIterator)()) while ((*particleIterator)())
{ {
G4ParticleDefinition* particle = particleIterator->value(); G4ParticleDefinition* particle = particleIterator->value();
G4String partname = particle->GetParticleName(); G4String partname = particle->GetParticleName();
if(partname == "p" || partname == "proton" || partname == "alpha" || partname == "He3" || partname == "GenericIon") { if(partname == "p" || partname == "proton" || partname == "alpha" || partname == "He3" || partname == "GenericIon") {
G4BraggIonGasModel* mod1 = new G4BraggIonGasModel(); G4BraggIonGasModel* mod1 = new G4BraggIonGasModel();
G4BetheBlochIonGasModel* mod2 = new G4BetheBlochIonGasModel(); G4BetheBlochIonGasModel* mod2 = new G4BetheBlochIonGasModel();
G4double eth = 2.*MeV*particle->GetPDGMass()/proton_mass_c2; G4double eth = 2.*MeV*particle->GetPDGMass()/proton_mass_c2;
emConfig->SetExtraEmModel(partname,"ionIoni",mod1,"",0.0,eth,new G4IonFluctuations()); emConfig->SetExtraEmModel(partname,"ionIoni",mod1,"",0.0,eth,new G4IonFluctuations());
emConfig->SetExtraEmModel(partname,"ionIoni",mod2,"",eth,100*TeV,new G4UniversalFluctuation()); emConfig->SetExtraEmModel(partname,"ionIoni",mod2,"",eth,100*TeV,new G4UniversalFluctuation());
} }
} }
} }
///////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////
void PhysicsList::AddPAIModel(const G4String& modname) void PhysicsList::AddPAIModel(const G4String& modname)
{ {
G4ParticleTable::G4PTblDicIterator* particleIterator = G4ParticleTable::GetParticleTable()->GetIterator(); G4ParticleTable::G4PTblDicIterator* particleIterator = G4ParticleTable::GetParticleTable()->GetIterator();
particleIterator->reset(); particleIterator->reset();
while ((*particleIterator)()) while ((*particleIterator)())
{ {
G4ParticleDefinition* particle = particleIterator->value(); G4ParticleDefinition* particle = particleIterator->value();
G4String partname = particle->GetParticleName(); G4String partname = particle->GetParticleName();
if(partname == "e-" || partname == "e+") { if(partname == "e-" || partname == "e+") {
NewPAIModel(particle, modname, "eIoni"); NewPAIModel(particle, modname, "eIoni");
} else if(partname == "mu-" || partname == "mu+") { } else if(partname == "mu-" || partname == "mu+") {
NewPAIModel(particle, modname, "muIoni"); NewPAIModel(particle, modname, "muIoni");
} else if(partname == "proton" || } else if(partname == "proton" ||
partname == "pi+" || partname == "pi+" ||
partname == "pi-" partname == "pi-"
) { ) {
NewPAIModel(particle, modname, "hIoni"); NewPAIModel(particle, modname, "hIoni");
}
} }
}
} }
///////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////
void PhysicsList::NewPAIModel(const G4ParticleDefinition* part, void PhysicsList::NewPAIModel(const G4ParticleDefinition* part,
const G4String& modname, const G4String& modname,
const G4String& procname) const G4String& procname)
{ {
G4String partname = part->GetParticleName(); G4String partname = part->GetParticleName();
if(modname == "pai"){ if(modname == "pai"){
G4PAIModel* pai = new G4PAIModel(part,"PAIModel"); G4PAIModel* pai = new G4PAIModel(part,"PAIModel");
emConfig->SetExtraEmModel(partname,procname,pai,"GasDetector",0.0,100.*TeV,pai); emConfig->SetExtraEmModel(partname,procname,pai,"GasDetector",0.0,100.*TeV,pai);
} }
else if(modname == "pai_photon") { else if(modname == "pai_photon") {
G4PAIPhotModel* pai = new G4PAIPhotModel(part,"PAIPhotModel"); G4PAIPhotModel* pai = new G4PAIPhotModel(part,"PAIPhotModel");
emConfig->SetExtraEmModel(partname,procname,pai,"GasDetector",0.0,100.*TeV,pai); emConfig->SetExtraEmModel(partname,procname,pai,"GasDetector",0.0,100.*TeV,pai);
} }
} }
///////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////
void PhysicsList::ConstructParticle(){ void PhysicsList::ConstructParticle(){
if(m_OpticalPhysics){ if(m_OpticalPhysics){
//G4OpticalPhoton::OpticalPhotonDefinition(); //G4OpticalPhoton::OpticalPhotonDefinition();
((G4VPhysicsConstructor*) opticalPhysicsList)->ConstructParticle(); ((G4VPhysicsConstructor*) opticalPhysicsList)->ConstructParticle();
} }
if(m_DriftElectronPhysics){ if(m_DriftElectronPhysics){
((G4VPhysicsConstructor*) driftElectronPhysicsList)->ConstructParticle(); ((G4VPhysicsConstructor*) driftElectronPhysicsList)->ConstructParticle();
} }
if(decay_List){ if(decay_List){
decay_List -> ConstructParticle(); decay_List -> ConstructParticle();
radioactiveDecay_List->ConstructParticle(); radioactiveDecay_List->ConstructParticle();
} }
else{ else{
// If decay is not activated we have to declare the particle ourself // If decay is not activated we have to declare the particle ourself
G4He3::He3Definition(); G4He3::He3Definition();
G4Deuteron::DeuteronDefinition(); G4Deuteron::DeuteronDefinition();
G4Triton::TritonDefinition(); G4Triton::TritonDefinition();
G4Alpha::AlphaDefinition(); G4Alpha::AlphaDefinition();
G4Proton::ProtonDefinition(); G4Proton::ProtonDefinition();
G4AntiProton::AntiProtonDefinition(); G4AntiProton::AntiProtonDefinition();
G4Neutron::NeutronDefinition(); G4Neutron::NeutronDefinition();
G4Electron::ElectronDefinition(); G4Electron::ElectronDefinition();
G4Positron::PositronDefinition(); G4Positron::PositronDefinition();
G4NeutrinoE::NeutrinoEDefinition(); G4NeutrinoE::NeutrinoEDefinition();
G4MuonPlus::MuonPlusDefinition(); G4MuonPlus::MuonPlusDefinition();
G4MuonMinus::MuonMinusDefinition(); G4MuonMinus::MuonMinusDefinition();
G4AntiNeutrinoE::AntiNeutrinoEDefinition(); G4AntiNeutrinoE::AntiNeutrinoEDefinition();
G4Geantino::GeantinoDefinition(); G4Geantino::GeantinoDefinition();
G4ChargedGeantino::ChargedGeantinoDefinition(); G4ChargedGeantino::ChargedGeantinoDefinition();
G4Gamma::GammaDefinition(); G4Gamma::GammaDefinition();
// mesons // mesons
G4PionPlus ::PionPlusDefinition(); G4PionPlus ::PionPlusDefinition();
G4PionMinus ::PionMinusDefinition(); G4PionMinus ::PionMinusDefinition();
G4PionZero ::PionZeroDefinition(); G4PionZero ::PionZeroDefinition();
G4Eta ::EtaDefinition(); G4Eta ::EtaDefinition();
G4EtaPrime ::EtaPrimeDefinition(); G4EtaPrime ::EtaPrimeDefinition();
// G4RhoZero ::RhoZeroDefinition(); // G4RhoZero ::RhoZeroDefinition();
G4KaonPlus ::KaonPlusDefinition(); G4KaonPlus ::KaonPlusDefinition();
G4KaonMinus ::KaonMinusDefinition(); G4KaonMinus ::KaonMinusDefinition();
G4KaonZero ::KaonZeroDefinition(); G4KaonZero ::KaonZeroDefinition();
G4AntiKaonZero ::AntiKaonZeroDefinition(); G4AntiKaonZero ::AntiKaonZeroDefinition();
G4KaonZeroLong ::KaonZeroLongDefinition(); G4KaonZeroLong ::KaonZeroLongDefinition();
G4KaonZeroShort::KaonZeroShortDefinition(); G4KaonZeroShort::KaonZeroShortDefinition();
// Ion // Ion
G4IonConstructor ionConstructor ; G4IonConstructor ionConstructor ;
ionConstructor.ConstructParticle() ; ionConstructor.ConstructParticle() ;
} }
} }
///////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////
void PhysicsList::ConstructProcess(){ void PhysicsList::ConstructProcess(){
// Transportation // Transportation
AddTransportation(); AddTransportation();
// Electromagnetic physics // Electromagnetic physics
emPhysicsList -> ConstructProcess(); emPhysicsList -> ConstructProcess();
if(opticalPhysicsList){ if(opticalPhysicsList){
((G4VPhysicsConstructor*) opticalPhysicsList)->ConstructProcess(); ((G4VPhysicsConstructor*) opticalPhysicsList)->ConstructProcess();
} }
if(driftElectronPhysicsList){ if(driftElectronPhysicsList){
((G4VPhysicsConstructor*) driftElectronPhysicsList)->ConstructProcess(); ((G4VPhysicsConstructor*) driftElectronPhysicsList)->ConstructProcess();
} }
// Hadronic physics // Hadronic physics
std::map<std::string,G4VPhysicsConstructor*>::iterator it; std::map<std::string,G4VPhysicsConstructor*>::iterator it;
for(it = m_PhysList.begin(); it!= m_PhysList.end(); it++){ for(it = m_PhysList.begin(); it!= m_PhysList.end(); it++){
it->second -> ConstructProcess(); it->second -> ConstructProcess();
} }
BiasCrossSectionByFactor(m_IonBinaryCascadePhysics); BiasCrossSectionByFactor(m_IonBinaryCascadePhysics);
em_option.SetBuildCSDARange(true); em_option.SetBuildCSDARange(true);
em_option.SetDEDXBinningForCSDARange(10*10); em_option.SetDEDXBinningForCSDARange(10*10);
em_option.SetFluo(true); em_option.SetFluo(true);
em_option.SetAuger(true); em_option.SetAuger(true);
AddParametrisation(); AddParametrisation();
AddLevelData();
return; return;
} }
///////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////
void PhysicsList::AddStepMax(){ void PhysicsList::AddStepMax(){
} }
///////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////
void PhysicsList::AddParametrisation(){ void PhysicsList::AddParametrisation(){
G4FastSimulationManagerProcess* BeamReaction = G4FastSimulationManagerProcess* BeamReaction =
new G4FastSimulationManagerProcess("NPSimulationProcess"); new G4FastSimulationManagerProcess("NPSimulationProcess");
// For 10.3 and higher // For 10.3 and higher
#ifndef theParticleIterator #ifndef theParticleIterator
G4ParticleTable::G4PTblDicIterator* theParticleIterator = GetParticleIterator(); G4ParticleTable::G4PTblDicIterator* theParticleIterator = GetParticleIterator();
#endif #endif
theParticleIterator->reset(); theParticleIterator->reset();
while ((*theParticleIterator)()){ while ((*theParticleIterator)()){
G4ParticleDefinition* particle = theParticleIterator->value(); G4ParticleDefinition* particle = theParticleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager(); G4ProcessManager* pmanager = particle->GetProcessManager();
std::string name = particle->GetParticleName(); std::string name = particle->GetParticleName();
pmanager->AddDiscreteProcess(BeamReaction); pmanager->AddDiscreteProcess(BeamReaction);
// Add a Step limiter to the beam particle. // Add a Step limiter to the beam particle.
// This will be used to limit the step of the beam in the target // This will be used to limit the step of the beam in the target
pmanager->AddProcess(new G4StepLimiter,-1,-1,5); pmanager->AddProcess(new G4StepLimiter,-1,-1,5);
// Add menate R to the process manager for neutrons // Add menate R to the process manager for neutrons
// (if selected as a option) // (if selected as a option)
if(m_Menate_R > 0 && name == "neutron") { if(m_Menate_R > 0 && name == "neutron") {
menate_R* theMenate = new menate_R("menateR_neutron"); menate_R* theMenate = new menate_R("menateR_neutron");
theMenate->SetMeanFreePathCalcMethod("ORIGINAL"); theMenate->SetMeanFreePathCalcMethod("ORIGINAL");
pmanager->AddDiscreteProcess(theMenate); pmanager->AddDiscreteProcess(theMenate);
std::cout <<"||--------------------------------------------------||" << std::endl; std::cout <<"||--------------------------------------------------||" << std::endl;
std::cout <<" MENATE_R Added to Process Manager " << std::endl; std::cout <<" MENATE_R Added to Process Manager " << std::endl;
std::cout <<"||--------------------------------------------------||" << std::endl; std::cout <<"||--------------------------------------------------||" << std::endl;
} }
} }
} }
////////////////////////////////////////////////////////////////////////////////
void PhysicsList::AddLevelData(){
NPL::InputParser parser(NPOptionManager::getInstance()->GetReactionFile());
vector<NPL::InputBlock*> blocks;
blocks = parser.GetAllBlocksWithToken("LevelData");
// Loop over all blocks
if (blocks.size()>0 && NPOptionManager::getInstance()->GetVerboseLevel())
cout << endl << "//// Reading user level data " << endl;
for (unsigned int i = 0; i < blocks.size(); i++) {
NPL::Nucleus N(blocks[i]->GetMainValue());
int Z =N.GetZ() ; int A=N.GetA();
cout << " Nuclei: " << N.GetName() << endl;
string level_path = blocks[i]->GetString("Path");
G4IonTable* ionTable = G4IonTable::GetIonTable();
G4NuclearLevelData* levelData = G4NuclearLevelData::GetInstance();
levelData->AddPrivateData(Z, A, level_path);
const G4LevelManager* levelManager = levelData->GetLevelManager(Z, A);
G4int Nentries = levelManager->NumberOfTransitions()+1;
for(G4int j = 1; j < Nentries; j++){ // Excited states
cout << " - Level " << j
<< " energy = " << levelManager->LevelEnergy(j)
<< endl;
G4ParticleDefinition* excitedState
= ionTable->GetIon(Z,A,levelManager->LevelEnergy(j));
excitedState->SetPDGStable(false);
excitedState->SetPDGLifeTime(levelManager->LifeTime(j));
}
}
}
///////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////
void PhysicsList::SetCuts(){ void PhysicsList::SetCuts(){
if (verboseLevel >0){ if (verboseLevel >0){
G4cout << "PhysicsList::SetCuts:"; G4cout << "PhysicsList::SetCuts:";
G4cout << "CutLength : " << G4BestUnit(defaultCutValue,"Length") << G4endl; G4cout << "CutLength : " << G4BestUnit(defaultCutValue,"Length") << G4endl;
} }
// Special Cut for optical photon to be emmitted // Special Cut for optical photon to be emmitted
SetCutsWithDefault(); SetCutsWithDefault();
} }
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
/////// Friend Method for CS biasing /////// Friend Method for CS biasing
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
void PhysicsList::BiasCrossSectionByFactor(double factor){ void PhysicsList::BiasCrossSectionByFactor(double factor){
factor++; factor++;
} }
NPSimulation/Process/PhysicsList.hh
+ 1
− 1
View file @ ad528b71
...@@ -89,7 +89,7 @@ class PhysicsList: public G4VUserPhysicsList{ ...@@ -89,7 +89,7 @@ class PhysicsList: public G4VUserPhysicsList{
void AddIonGasModels(); void AddIonGasModels();
void AddPAIModel(const G4String& modname); void AddPAIModel(const G4String& modname);
void NewPAIModel(const G4ParticleDefinition* part, const G4String& modname,const G4String& procname); void NewPAIModel(const G4ParticleDefinition* part, const G4String& modname,const G4String& procname);
void AddLevelData();
private: private:
std::map<std::string,G4VPhysicsConstructor*> m_PhysList; std::map<std::string,G4VPhysicsConstructor*> m_PhysList;
G4EmConfigurator em_config; G4EmConfigurator em_config;
......
Projects/e748/Analysis.cxx
+ 194
− 194
View file @ ad528b71
This diff is collapsed.
Projects/e748/Analysis.h
+ 4
− 4
View file @ ad528b71
...@@ -74,10 +74,10 @@ class Analysis: public NPL::VAnalysis{ ...@@ -74,10 +74,10 @@ class Analysis: public NPL::VAnalysis{
double BeamEnergy; double BeamEnergy;
int run_major; int run_major;
int run_minor; int run_minor;
NPL::EnergyLoss LightTarget; NPL::EnergyLoss He3CD2 ;
NPL::EnergyLoss LightAl ; NPL::EnergyLoss He3Al ;
NPL::EnergyLoss LightSi ; NPL::EnergyLoss He3Si ;
NPL::EnergyLoss BeamTarget ; NPL::EnergyLoss BeamCD2 ;
NPL::EnergyLoss BeamMylar ; NPL::EnergyLoss BeamMylar ;
NPL::EnergyLoss BeamIsobutane ; NPL::EnergyLoss BeamIsobutane ;
......
Projects/e748/configs/ConfigMust2.dat 100644 → 100755
+ 0
− 0
View file @ ad528b71
File mode changed from 100644 to 100755
Projects/e748/e748.detector 100644 → 100755
+ 1
− 1
View file @ ad528b71
...@@ -116,7 +116,7 @@ CATSDetector ...@@ -116,7 +116,7 @@ CATSDetector
ModularLeaf ModularLeaf
DefaultValue= -1000 DefaultValue= -1000
Leafs= QPlast QCaviar CHIO T_CATS1_CAV T_CATS1_2 T_MUVI_CATS1 T_PL_CATS2 EXL_HF T_PL_CHIO T_PL_CATS E1D6 Leafs= QPlast QCaviar CHIO T_CATS1_CAV T_CATS1_2 T_MUVI_CATS1 T_PL_CATS2 EXL_HF T_PL_CHIO T_PL_CATS1
Chio_dig Chio_dig
Pos= -30. -30. 600. mm Pos= -30. -30. 600. mm
......
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