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 //
0030 //    ******      SUGGESTED PHYSICS FOR ACCURATE SIMULATIONS    *********
0031 //    ******            IN MEDICAL PHYSICS APPLICATIONS         *********
0032 //
0033 // 'HADRONTHERAPY_1' is more suited for protons only
0034 // 'HADRONTHERAPY_2' is suggested for better precision with ions
0035 // 'HADRONTHERAPY_3' test that uses Bertini cascade
0036 // It can be activated inside any macro file using the command:
0037 // /Physics/addPhysics HADRONTHERAPY_1 (HADRONTHERAPY_2) (HADRONTHERAPY_3)
0038
0039 #include "G4SystemOfUnits.hh"
0040 #include "G4RunManager.hh"
0041 #include "G4Region.hh"
0042 #include "G4RegionStore.hh"
0043 #include "HadrontherapyPhysicsList.hh"
0044 #include "HadrontherapyPhysicsListMessenger.hh"
0045 #include "HadrontherapyStepMax.hh"
0046 #include "G4PhysListFactory.hh"
0047 #include "G4VPhysicsConstructor.hh"
0048 #include "G4HadronPhysicsQGSP_BIC_HP.hh"
0049 #include "G4HadronPhysicsQGSP_BIC.hh"
0050 #include "G4EmStandardPhysics_option4.hh"
0051 #include "G4EmStandardPhysics.hh"
0052 #include "G4EmExtraPhysics.hh"
0053 #include "G4StoppingPhysics.hh"
0054 #include "G4DecayPhysics.hh"
0055 #include "G4HadronElasticPhysics.hh"
0056 #include "G4HadronElasticPhysicsHP.hh"
0057 #include "G4RadioactiveDecayPhysics.hh"
0058 #include "G4IonBinaryCascadePhysics.hh"
0059 #include "G4DecayPhysics.hh"
0060 #include "G4NeutronTrackingCut.hh"
0061 #include "G4LossTableManager.hh"
0062 #include "G4UnitsTable.hh"
0063 #include "G4ProcessManager.hh"
0064 #include "G4IonFluctuations.hh"
0065 #include "G4IonParametrisedLossModel.hh"
0066 #include "G4EmParameters.hh"
0067 #include "G4ParallelWorldPhysics.hh"
0068 #include "G4EmLivermorePhysics.hh"
0069 #include "G4AutoDelete.hh"
0070 #include "G4HadronPhysicsQGSP_BIC_AllHP.hh"
0071 #include "QGSP_BIC_HP.hh"
0072 #include "G4HadronPhysicsQGSP_BERT.hh"
0073 #include "G4HadronPhysicsQGSP_BERT_HP.hh"
0074
0075 /////////////////////////////////////////////////////////////////////////////
0076 HadrontherapyPhysicsList::HadrontherapyPhysicsList() : G4VModularPhysicsList()
0077 {
0078     G4LossTableManager::Instance();
0079     defaultCutValue = 1.*mm;
0080     cutForGamma     = defaultCutValue;
0081     cutForElectron  = defaultCutValue;
0082     cutForPositron  = defaultCutValue;
0083
0084     pMessenger = new HadrontherapyPhysicsListMessenger(this);
0085     SetVerboseLevel(1);
0086     decay_List = new G4DecayPhysics();
0087     // Elecromagnetic physics
0088     //
0089     emPhysicsList = new G4EmStandardPhysics_option4();
0090
0091 }
0092
0093 /////////////////////////////////////////////////////////////////////////////
0094 HadrontherapyPhysicsList::~HadrontherapyPhysicsList()
0095 {
0096     delete pMessenger;
0097     delete emPhysicsList;
0098     delete decay_List;
0099     //delete radioactiveDecay_List;
0100     hadronPhys.clear();
0101     for(size_t i=0; i<hadronPhys.size(); i++)
0102     {
0103         delete hadronPhys[i];
0104     }
0105 }
0106
0107 /////////////////////////////////////////////////////////////////////////////
0108 void HadrontherapyPhysicsList::ConstructParticle()
0109 {
0110     decay_List -> ConstructParticle();
0111
0112 }
0113
0114 /////////////////////////////////////////////////////////////////////////////
0115 void HadrontherapyPhysicsList::ConstructProcess()
0116 {
0117     // Transportation
0118     //
0119     AddTransportation();
0120
0121     decay_List -> ConstructProcess();
0122     emPhysicsList -> ConstructProcess();
0123
0124
0125     //em_config.AddModels();
0126
0127     // Hadronic physics
0128     //
0129     for(size_t i=0; i < hadronPhys.size(); i++)
0130     {
0131         hadronPhys[i] -> ConstructProcess();
0132     }
0133
0134     // step limitation (as a full process)
0135     //
0136     AddStepMax();
0137
0138     //Parallel world sensitivity
0139     //
0140     G4ParallelWorldPhysics* pWorld = new G4ParallelWorldPhysics("DetectorROGeometry");
0141     pWorld->ConstructProcess();
0142
0143     return;
0144 }
0145
0146 /////////////////////////////////////////////////////////////////////////////
0147 void HadrontherapyPhysicsList::AddPhysicsList(const G4String& name)
0148 {
0149     if (verboseLevel>1) {
0150         G4cout << "PhysicsList::AddPhysicsList: <" << name << ">" << G4endl;
0151     }
0152     if (name == emName) return;
0153
0154     ///////////////////////////////////
0155     //   ELECTROMAGNETIC MODELS
0156     ///////////////////////////////////
0157     if (name == "standard_opt4") {
0158         emName = name;
0159         delete emPhysicsList;
0160         hadronPhys.clear();
0161         emPhysicsList = new G4EmStandardPhysics_option4();
0162         G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
0163         G4cout << "THE FOLLOWING ELECTROMAGNETIC PHYSICS LIST HAS BEEN ACTIVATED: G4EmStandardPhysics_option4" << G4endl;
0164
0165         ////////////////////////////////////////
0166         //   ELECTROMAGNETIC + HADRONIC MODELS
0167         ////////////////////////////////////////
0168
0169     }  else if (name == "HADRONTHERAPY_1") {
0170
0171         AddPhysicsList("standard_opt4");
0172         hadronPhys.push_back( new G4DecayPhysics());
0173         hadronPhys.push_back( new G4RadioactiveDecayPhysics());
0174         hadronPhys.push_back( new G4IonBinaryCascadePhysics());
0175         hadronPhys.push_back( new G4EmExtraPhysics());
0176         hadronPhys.push_back( new G4HadronElasticPhysicsHP());
0177         hadronPhys.push_back( new G4StoppingPhysics());
0178         hadronPhys.push_back( new G4HadronPhysicsQGSP_BIC_HP());
0179         hadronPhys.push_back( new G4NeutronTrackingCut());
0180
0181         G4cout << "HADRONTHERAPY_1 PHYSICS LIST has been activated" << G4endl;
0182     }
0183
0184     else if (name == "HADRONTHERAPY_2") {
0185
0186         AddPhysicsList("standard_opt4");
0187         hadronPhys.push_back( new G4DecayPhysics());
0188         hadronPhys.push_back( new G4RadioactiveDecayPhysics());
0189         hadronPhys.push_back( new G4IonBinaryCascadePhysics());
0190         hadronPhys.push_back( new G4EmExtraPhysics());
0191         hadronPhys.push_back( new G4HadronElasticPhysics());
0192         hadronPhys.push_back( new G4StoppingPhysics());
0193         hadronPhys.push_back( new G4HadronPhysicsQGSP_BIC_AllHP());
0194         hadronPhys.push_back( new G4NeutronTrackingCut());
0195
0196         G4cout << "HADRONTHERAPY_2 PHYSICS LIST has been activated" << G4endl;
0197     }
0198
0199     else if (name == "HADRONTHERAPY_3"){
0200         AddPhysicsList("standard_opt4");
0201         hadronPhys.push_back( new G4DecayPhysics());
0202         hadronPhys.push_back( new G4RadioactiveDecayPhysics());
0203         hadronPhys.push_back( new G4IonBinaryCascadePhysics());
0204         hadronPhys.push_back( new G4EmExtraPhysics());
0205         hadronPhys.push_back( new G4HadronElasticPhysics());
0206         hadronPhys.push_back( new G4StoppingPhysics());
0207         hadronPhys.push_back( new G4HadronPhysicsQGSP_BERT_HP());
0208         hadronPhys.push_back( new G4NeutronTrackingCut());
0209
0210         G4cout << "HADRONTHERAPY_3 PHYSICS LIST has been activated" << G4endl;
0211     }
0212
0213     else {
0214         G4cout << "PhysicsList::AddPhysicsList: <" << name << ">"
0215         << " is not defined"
0216         << G4endl;
0217     }
0218
0219 }
0220
0221 /////////////////////////////////////////////////////////////////////////////
0222 void HadrontherapyPhysicsList::AddStepMax()
0223 {
0224     // Step limitation seen as a process
0225     // This process must exist in all threads.
0226     //
0227     HadrontherapyStepMax* stepMaxProcess  = new HadrontherapyStepMax();
0228
0229
0230     auto particleIterator = GetParticleIterator();
0231     particleIterator->reset();
0232     while ((*particleIterator)()){
0233         G4ParticleDefinition* particle = particleIterator->value();
0234         G4ProcessManager* pmanager = particle->GetProcessManager();
0235
0236         if (stepMaxProcess->IsApplicable(*particle) && pmanager)
0237         {
0238             pmanager ->AddDiscreteProcess(stepMaxProcess);
0239         }
0240     }
0241 }