iort_therapy/src/IORTPhysicsList.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 // This is the *BASIC* version of IORT, a Geant4-based application
0027 //
0028 // Main Authors: G.Russo(a,b), C.Casarino*(c), G.C. Candiano(c), G.A.P. Cirrone(d), F.Romano(d)
0029 // Contributor Authors: S.Guatelli(e)
0030 // Past Authors: G.Arnetta(c), S.E.Mazzaglia(d)
0031 //
0032 //   (a) Fondazione Istituto San Raffaele G.Giglio, Cefalù, Italy
0033 //   (b) IBFM-CNR , Segrate (Milano), Italy
0034 //   (c) LATO (Laboratorio di Tecnologie Oncologiche), Cefalù, Italy
0035 //   (d) Laboratori Nazionali del Sud of the INFN, Catania, Italy
0036 //   (e) University of Wollongong, Australia
0037 //
0038 //   *Corresponding author, email to carlo.casarino@polooncologicocefalu.it
0039 //////////////////////////////////////////////////////////////////////////////////////////////
0040 //
0041
0042 #include "G4SystemOfUnits.hh"
0043 #include "G4RunManager.hh"
0044 #include "G4Region.hh"
0045 #include "G4RegionStore.hh"
0046 #include "IORTPhysicsList.hh"
0047 #include "IORTPhysicsListMessenger.hh"
0048 #include "IORTStepMax.hh"
0049 #include "G4VPhysicsConstructor.hh"
0050
0051 // Physic lists (contained inside the Geant4 source code, in the 'physicslists folder')
0052 #include "G4EmStandardPhysics_option3.hh"
0053 #include "G4EmStandardPhysics_option4.hh"
0054 #include "G4EmLivermorePhysics.hh"
0055 #include "G4EmPenelopePhysics.hh"
0056 #include "G4EmExtraPhysics.hh"
0057 #include "G4ParticleDefinition.hh"
0058 #include "G4ProductionCutsTable.hh"
0059 #include "G4ProcessManager.hh"
0060 #include "globals.hh"
0061 #include "G4Electron.hh"
0062 #include "G4Gamma.hh"
0063 #include "G4Positron.hh"
0064 #include "G4UnitsTable.hh"
0065 #include "G4DecayPhysics.hh"
0066
0067 IORTPhysicsList::IORTPhysicsList() : G4VModularPhysicsList()
0068 {
0069   defaultCutValue = 0.1 *mm;
0070   cutForGamma     = defaultCutValue;
0071   cutForElectron  = defaultCutValue;
0072   cutForPositron  = defaultCutValue;
0073
0074   // set cut values for gamma at first and for e- second and next for e+,
0075   // because some processes for e+/e- need cut values for gamma
0076   SetCutValue(cutForGamma, "gamma");
0077   SetCutValue(cutForElectron, "e-");
0078   SetCutValue(cutForPositron, "e+");
0079
0080   DumpCutValuesTable();
0081
0082   stepMaxProcess  = 0;
0083
0084   pMessenger = new IORTPhysicsListMessenger(this);
0085
0086   SetVerboseLevel(1);
0087
0088   // EM physics
0089   emPhysicsList = new G4EmStandardPhysics_option4(1);
0090   emName = G4String("emstandard_opt4");
0091   decPhysicsList = new G4DecayPhysics();
0092 }
0093
0094 IORTPhysicsList::~IORTPhysicsList()
0095 {
0096   delete pMessenger;
0097   delete emPhysicsList;
0098   delete decPhysicsList;
0099 }
0100
0101 void IORTPhysicsList::ConstructParticle()
0102 {
0103   decPhysicsList->ConstructParticle();
0104 }
0105
0106 void IORTPhysicsList::ConstructProcess()
0107 {
0108   // transportation
0109   AddTransportation();
0110
0111   // electromagnetic physics list
0112   emPhysicsList->ConstructProcess();
0113   em_config.AddModels();
0114
0115   // Set cuts for detector
0116   SetDetectorCut(defaultCutValue);
0117   // step limitation (as a full process)
0118   //
0119   AddStepMax();
0120 }
0121
0122 void IORTPhysicsList::AddPhysicsList(const G4String& name)
0123 {
0124
0125     if (verboseLevel>1) {
0126     G4cout << "PhysicsList::AddPhysicsList: <" << name << ">" << G4endl;
0127     }
0128     if (name == emName) return;
0129
0130   if (name == "standard_opt3") {
0131     emName = name;
0132     delete emPhysicsList;
0133     emPhysicsList = new G4EmStandardPhysics_option3();
0134     G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
0135     G4cout << "THE FOLLOWING ELECTROMAGNETIC PHYSICS LIST HAS BEEN ACTIVATED: G4EmStandardPhysics_option3" << G4endl;
0136
0137
0138     } else if (name == "livermore") {
0139     emName = name;
0140     delete emPhysicsList;
0141     emPhysicsList = new G4EmLivermorePhysics();
0142     G4RunManager::GetRunManager()-> PhysicsHasBeenModified();
0143     G4cout << "THE FOLLOWING ELECTROMAGNETIC PHYSICS LIST HAS BEEN ACTIVATED: G4EmLivermorePhysics" << G4endl;
0144
0145     } else if (name == "penelope") {
0146     emName = name;
0147     delete emPhysicsList;
0148     emPhysicsList = new G4EmPenelopePhysics();
0149     G4RunManager::GetRunManager()-> PhysicsHasBeenModified();
0150     G4cout << "THE FOLLOWING ELECTROMAGNETIC PHYSICS LIST HAS BEEN ACTIVATED: G4EmPenelopePhysics" << G4endl;}
0151
0152  else if (name == "standard_opt4") {
0153     emName = name;
0154     delete emPhysicsList;
0155     emPhysicsList = new G4EmStandardPhysics_option4();
0156     G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
0157     G4cout << "THE FOLLOWING ELECTROMAGNETIC PHYSICS LIST HAS BEEN ACTIVATED: G4EmStandardPhysics_option4" << G4endl;
0158         }
0159
0160     else {
0161
0162     G4cout << "PhysicsList::AddPhysicsList: <" << name << ">"
0163         << " is not defined"
0164         << G4endl;
0165     }
0166 }
0167
0168 void IORTPhysicsList::AddStepMax()
0169 {
0170   // Step limitation seen as a process
0171   stepMaxProcess = new IORTStepMax();
0172
0173   auto particleIterator=GetParticleIterator();
0174   particleIterator->reset();
0175   while ((*particleIterator)()){
0176     G4ParticleDefinition* particle = particleIterator->value();
0177     G4ProcessManager* pmanager = particle->GetProcessManager();
0178
0179     if (stepMaxProcess->IsApplicable(*particle) && pmanager)
0180       {
0181     pmanager ->AddDiscreteProcess(stepMaxProcess);
0182       }
0183   }
0184 }
0185
0186 void IORTPhysicsList::SetCutForGamma(G4double cut)
0187 {
0188   cutForGamma = cut;
0189   SetParticleCuts(cutForGamma, G4Gamma::Gamma());
0190 }
0191
0192 void IORTPhysicsList::SetCutForElectron(G4double cut)
0193 {
0194   cutForElectron = cut;
0195   SetParticleCuts(cutForElectron, G4Electron::Electron());
0196 }
0197
0198 void IORTPhysicsList::SetCutForPositron(G4double cut)
0199 {
0200   cutForPositron = cut;
0201   SetParticleCuts(cutForPositron, G4Positron::Positron());
0202 }
0203
0204 void IORTPhysicsList::SetDetectorCut(G4double cut)
0205 {
0206   G4String regionName = "DetectorLog";
0207   G4Region* region = G4RegionStore::GetInstance()->GetRegion(regionName);
0208
0209   G4ProductionCuts* cuts = new G4ProductionCuts ;
0210   cuts -> SetProductionCut(cut,"gamma");
0211   cuts -> SetProductionCut(cut,"e-");
0212   cuts -> SetProductionCut(cut,"e+");
0213   region -> SetProductionCuts(cuts);
0214 }
0215