hadrontherapy/src/HadrontherapyPhysicsList.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // Hadrontherapy advanced example for Geant4
0027 // See more at: https://twiki.cern.ch/twiki/bin/view/Geant4/AdvancedExamplesHadrontherapy
0028 //
0029 // Using the builder concepts of Geant4 we assembled (and tested) two different
0030 // Physics Lists that are particuilarly suited for Hadronterapy applications:
0031 //
0032 // 'HADRONTHERAPY_1' is more suited for protons only
0033 // 'HADRONTHERAPY_2' is suggested for better precision with ions
0034 //
0035 // The Reference physics lists (already present in the Geant4 kernel) can
0036 // be used as well. In this case the more suitable "Reference physics lists" are:
0037 // "QBBC", "QGSP_BIC", "Shielding", "QGSP_BERT",
0038 // "QGSP_BIC_AllHP" and "QGSP_BIC_HP"
0039 //
0040 // NOTE: to activate the "_HP" physics you have to set the G4PARTICLEHPDATA environment
0041 // variable pointing to the external dataset named "G4TENDL".
0042 //
0043 // All the lists can be  activated inside any macro file using the command:
0044 // /Physics/addPhysics
0045 //
0046 // Examples of usage are:
0047 // /Physics/addPhysics HADRONTHERAPY_1 or /Physics/addPhysics QGSP_BIC_HP
0048
0049 #include "G4SystemOfUnits.hh"
0050 #include "G4RunManager.hh"
0051 #include "G4Region.hh"
0052 #include "G4RegionStore.hh"
0053 #include "HadrontherapyPhysicsList.hh"
0054 #include "HadrontherapyPhysicsListMessenger.hh"
0055 #include "HadrontherapyStepMax.hh"
0056 #include "G4PhysListFactory.hh"
0057 #include "G4VPhysicsConstructor.hh"
0058 #include "G4HadronPhysicsQGSP_BIC_HP.hh"
0059 #include "G4HadronPhysicsQGSP_BIC.hh"
0060 #include "G4EmStandardPhysics_option4.hh"
0061 #include "G4EmStandardPhysics.hh"
0062 #include "G4EmExtraPhysics.hh"
0063 #include "G4StoppingPhysics.hh"
0064 #include "G4DecayPhysics.hh"
0065 #include "G4HadronElasticPhysics.hh"
0066 #include "G4HadronElasticPhysicsHP.hh"
0067 #include "G4RadioactiveDecayPhysics.hh"
0068 #include "G4IonBinaryCascadePhysics.hh"
0069 #include "G4DecayPhysics.hh"
0070 #include "G4NeutronTrackingCut.hh"
0071 #include "G4LossTableManager.hh"
0072 #include "G4UnitsTable.hh"
0073 #include "G4ProcessManager.hh"
0074 #include "G4IonFluctuations.hh"
0075 #include "G4IonParametrisedLossModel.hh"
0076 #include "G4EmParameters.hh"
0077 #include "G4ParallelWorldPhysics.hh"
0078 #include "G4EmLivermorePhysics.hh"
0079 #include "G4AutoDelete.hh"
0080 #include "G4HadronPhysicsQGSP_BIC_AllHP.hh"
0081 #include "QGSP_BIC_HP.hh"
0082 #include "QGSP_BIC.hh"
0083 #include "G4HadronPhysicsQGSP_BERT.hh"
0084 #include "G4HadronPhysicsQGSP_BERT_HP.hh"
0085 #include "G4ParallelWorldPhysics.hh"
0086 // Physics List
0087 #include "QBBC.hh"
0088 #include "QGSP_BIC.hh"
0089 #include "Shielding.hh"
0090 #include "QGSP_BERT.hh"
0091 #include "QGSP_BIC_AllHP.hh"
0092 #include "QGSP_BIC_HP.hh"
0093
0094
0095
0096 /////////////////////////////////////////////////////////////////////////////
0097 HadrontherapyPhysicsList::HadrontherapyPhysicsList() : G4VModularPhysicsList()
0098 {
0099     G4LossTableManager::Instance();
0100     defaultCutValue = 1.*mm;
0101     cutForGamma     = defaultCutValue;
0102     cutForElectron  = defaultCutValue;
0103     cutForPositron  = defaultCutValue;
0104
0105     pMessenger = new HadrontherapyPhysicsListMessenger(this);
0106     SetVerboseLevel(1);
0107     decay_List = new G4DecayPhysics();
0108     // Elecromagnetic physics
0109     //
0110     emPhysicsList = new G4EmStandardPhysics_option4();
0111 }
0112
0113 /////////////////////////////////////////////////////////////////////////////
0114 HadrontherapyPhysicsList::~HadrontherapyPhysicsList()
0115 {
0116     delete pMessenger;
0117     delete emPhysicsList;
0118     delete decay_List;
0119     //delete radioactiveDecay_List;
0120     hadronPhys.clear();
0121     for(size_t i=0; i<hadronPhys.size(); i++)
0122     {
0123         delete hadronPhys[i];
0124     }
0125 }
0126
0127 /////////////////////////////////////////////////////////////////////////////
0128 void HadrontherapyPhysicsList::ConstructParticle()
0129 {
0130     decay_List -> ConstructParticle();
0131
0132 }
0133
0134 /////////////////////////////////////////////////////////////////////////////
0135 void HadrontherapyPhysicsList::ConstructProcess()
0136 {
0137     // Transportation
0138     //
0139     AddTransportation();
0140
0141     decay_List -> ConstructProcess();
0142     emPhysicsList -> ConstructProcess();
0143
0144
0145     //em_config.AddModels();
0146
0147     // Hadronic physics
0148     //
0149     for(size_t i=0; i < hadronPhys.size(); i++)
0150     {
0151         hadronPhys[i] -> ConstructProcess();
0152     }
0153
0154     // step limitation (as a full process)
0155     //
0156     AddStepMax();
0157
0158     //Parallel world sensitivity
0159     //
0160     G4ParallelWorldPhysics* pWorld = new G4ParallelWorldPhysics("DetectorROGeometry");
0161     pWorld->ConstructProcess();
0162
0163     return;
0164 }
0165
0166 /////////////////////////////////////////////////////////////////////////////
0167 void HadrontherapyPhysicsList::AddPhysicsList(const G4String& name)
0168 {
0169     if (verboseLevel>1) {
0170         G4cout << "PhysicsList::AddPhysicsList: <" << name << ">" << G4endl;
0171     }
0172     if (name == emName) return;
0173
0174     ///////////////////////////////////
0175     //   ELECTROMAGNETIC MODELS
0176     ///////////////////////////////////
0177     if (name == "standard_opt4") {
0178         emName = name;
0179         delete emPhysicsList;
0180         hadronPhys.clear();
0181         emPhysicsList = new G4EmStandardPhysics_option4();
0182         G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
0183         G4cout << "THE FOLLOWING ELECTROMAGNETIC PHYSICS LIST HAS BEEN ACTIVATED: G4EmStandardPhysics_option4" << G4endl;
0184
0185         ////////////////////////////////////////
0186         //   ELECTROMAGNETIC + HADRONIC MODELS
0187         ////////////////////////////////////////
0188
0189     }  else if (name == "HADRONTHERAPY_1") {
0190
0191         AddPhysicsList("standard_opt4");
0192         hadronPhys.push_back( new G4DecayPhysics());
0193         hadronPhys.push_back( new G4RadioactiveDecayPhysics());
0194         hadronPhys.push_back( new G4IonBinaryCascadePhysics());
0195         hadronPhys.push_back( new G4EmExtraPhysics());
0196         hadronPhys.push_back( new G4HadronElasticPhysicsHP());
0197         hadronPhys.push_back( new G4StoppingPhysics());
0198         hadronPhys.push_back( new G4HadronPhysicsQGSP_BIC_HP());
0199         hadronPhys.push_back( new G4NeutronTrackingCut());
0200
0201         G4cout << "HADRONTHERAPY_1 PHYSICS LIST has been activated" << G4endl;
0202     }
0203
0204     else if (name == "HADRONTHERAPY_2") {
0205
0206         AddPhysicsList("standard_opt4");
0207         hadronPhys.push_back( new G4DecayPhysics());
0208         hadronPhys.push_back( new G4RadioactiveDecayPhysics());
0209         hadronPhys.push_back( new G4IonBinaryCascadePhysics());
0210         hadronPhys.push_back( new G4EmExtraPhysics());
0211         hadronPhys.push_back( new G4HadronElasticPhysics());
0212         hadronPhys.push_back( new G4StoppingPhysics());
0213         hadronPhys.push_back( new G4HadronPhysicsQGSP_BIC_AllHP());
0214         hadronPhys.push_back( new G4NeutronTrackingCut());
0215
0216         G4cout << "HADRONTHERAPY_2 PHYSICS LIST has been activated" << G4endl;
0217
0218     }
0219
0220     else if (name == "QGSP_BIC"){
0221         auto physicsList = new QGSP_BIC;
0222         G4RunManager::GetRunManager() -> SetUserInitialization(physicsList);
0223         G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
0224         physicsList -> RegisterPhysics(new G4ParallelWorldPhysics("DetectorROGeometry"));
0225     }
0226
0227     else if (name == "QGSP_BERT"){
0228         auto physicsList = new QGSP_BERT;
0229         G4RunManager::GetRunManager() -> SetUserInitialization(physicsList);
0230         G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
0231         physicsList -> RegisterPhysics(new G4ParallelWorldPhysics("DetectorROGeometry"));
0232     }
0233
0234     else if (name == "QGSP_BIC_AllHP"){
0235         auto physicsList = new QGSP_BIC_AllHP;
0236         G4RunManager::GetRunManager() -> SetUserInitialization(physicsList);
0237         G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
0238         physicsList -> RegisterPhysics(new G4ParallelWorldPhysics("DetectorROGeometry"));
0239     }
0240
0241     else if (name == "QGSP_BIC_HP"){
0242         auto physicsList = new QGSP_BIC_HP;
0243         G4RunManager::GetRunManager() -> SetUserInitialization(physicsList);
0244         G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
0245         physicsList -> RegisterPhysics(new G4ParallelWorldPhysics("DetectorROGeometry"));
0246     }
0247
0248     else if (name == "Shielding"){
0249         auto physicsList = new Shielding;
0250         G4RunManager::GetRunManager() -> SetUserInitialization(physicsList);
0251         G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
0252         physicsList -> RegisterPhysics(new G4ParallelWorldPhysics("DetectorROGeometry"));
0253     }
0254
0255     else if (name == "QBBC"){
0256         auto physicsList = new QBBC;
0257         G4RunManager::GetRunManager() -> SetUserInitialization(physicsList);
0258         G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
0259         physicsList -> RegisterPhysics(new G4ParallelWorldPhysics("DetectorROGeometry"));
0260     }
0261
0262     else {
0263         G4cout << "PhysicsList::AddPhysicsList: <" << name << ">"
0264         << " is not defined"
0265         << G4endl;
0266     }
0267
0268 }
0269
0270 /////////////////////////////////////////////////////////////////////////////
0271 void HadrontherapyPhysicsList::AddStepMax()
0272 {
0273     // Step limitation seen as a process
0274     // This process must exist in all threads.
0275     //
0276     HadrontherapyStepMax* stepMaxProcess  = new HadrontherapyStepMax();
0277
0278
0279     auto particleIterator = GetParticleIterator();
0280     particleIterator->reset();
0281     while ((*particleIterator)()){
0282         G4ParticleDefinition* particle = particleIterator->value();
0283         G4ProcessManager* pmanager = particle->GetProcessManager();
0284
0285         if (stepMaxProcess->IsApplicable(*particle) && pmanager)
0286         {
0287             pmanager ->AddDiscreteProcess(stepMaxProcess);
0288         }
0289     }
0290 }