| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // Hadrontherapy advanced example for Geant4 << 26 // This is the *BASIC* version of Hadrontherapy, a Geant4-based application
27 // See more at: https://twiki.cern.ch/twiki/bi << 27 // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
28 // 28 //
29 // Using the builder concepts of Geant4 we ass << 29 // Visit the Hadrontherapy web site (http://www.lns.infn.it/link/Hadrontherapy) to request
30 // Physics Lists that are particuilarly suited << 30 // the *COMPLETE* version of this program, together with its documentation;
>> 31 // Hadrontherapy (both basic and full version) are supported by the Italian INFN
>> 32 // Institute in the framework of the MC-INFN Group
31 // 33 //
32 // 'HADRONTHERAPY_1' is more suited for proton <<
33 // 'HADRONTHERAPY_2' is suggested for better p <<
34 // 34 //
35 // The Reference physics lists (already presen << 35 // Hadrontherapy Users are recommended to use the prepared macro files in order to activate the
36 // be used as well. In this case the more suit << 36 // most appropriate physics for hadrontherapy applications.
37 // "QBBC", "QGSP_BIC", "Shielding", "QGSP_BERT << 37 // As one can easily see the physics we suggest is contained in the
38 // "QGSP_BIC_AllHP" and "QGSP_BIC_HP" << 38 // 'QGSP_BIC_EMY' list adding to the decay processes (activated as default).
39 // 39 //
40 // NOTE: to activate the "_HP" physics you hav <<
41 // variable pointing to the external dataset n <<
42 // 40 //
43 // All the lists can be activated inside any << 41 // ****** SUGGESTED PHYSICS *********
44 // /Physics/addPhysics <<
45 // 42 //
46 // Examples of usage are: << 43 // AT MOMENT, IF ACCURATE RESULTS ARE NEEDED, WE STRONGLY RECOMMEND:
47 // /Physics/addPhysics HADRONTHERAPY_1 or /Phy << 44 // 1. The use of the macro 'hadron_therapy.mac', that is connected with the HadrontherapyPhysicsList.cc file.
48 << 45 // 2. the QGSP_BIC_EMY Reference Physics Lists (define the PHYSLIST eviroment variable):
49 #include "G4SystemOfUnits.hh" << 46 // export PHYSLIST=QGSP_BIC_EMY
>> 47 // User must considered that, in this second case, radioactive processes are not activated
>> 48
50 #include "G4RunManager.hh" 49 #include "G4RunManager.hh"
51 #include "G4Region.hh" 50 #include "G4Region.hh"
52 #include "G4RegionStore.hh" 51 #include "G4RegionStore.hh"
53 #include "HadrontherapyPhysicsList.hh" 52 #include "HadrontherapyPhysicsList.hh"
54 #include "HadrontherapyPhysicsListMessenger.hh 53 #include "HadrontherapyPhysicsListMessenger.hh"
55 #include "HadrontherapyStepMax.hh" 54 #include "HadrontherapyStepMax.hh"
56 #include "G4PhysListFactory.hh" 55 #include "G4PhysListFactory.hh"
57 #include "G4VPhysicsConstructor.hh" 56 #include "G4VPhysicsConstructor.hh"
58 #include "G4HadronPhysicsQGSP_BIC_HP.hh" << 57
59 #include "G4HadronPhysicsQGSP_BIC.hh" << 58 // Local physic directly implemented in the Hadronthrapy directory
60 #include "G4EmStandardPhysics_option4.hh" << 59 #include "LocalIonIonInelasticPhysic.hh" // Physic dedicated to the ion-ion inelastic processes
61 #include "G4EmStandardPhysics.hh" << 60
>> 61 // Physic lists (contained inside the Geant4 source code, in the 'physicslists folder')
>> 62 #include "HadronPhysicsQGSP_BIC.hh"
>> 63 #include "G4EmStandardPhysics_option3.hh"
>> 64 #include "G4EmLivermorePhysics.hh"
>> 65 #include "G4EmPenelopePhysics.hh"
62 #include "G4EmExtraPhysics.hh" 66 #include "G4EmExtraPhysics.hh"
63 #include "G4StoppingPhysics.hh" << 67 #include "G4QStoppingPhysics.hh"
64 #include "G4DecayPhysics.hh" 68 #include "G4DecayPhysics.hh"
65 #include "G4HadronElasticPhysics.hh" 69 #include "G4HadronElasticPhysics.hh"
66 #include "G4HadronElasticPhysicsHP.hh" 70 #include "G4HadronElasticPhysicsHP.hh"
67 #include "G4RadioactiveDecayPhysics.hh" 71 #include "G4RadioactiveDecayPhysics.hh"
68 #include "G4IonBinaryCascadePhysics.hh" 72 #include "G4IonBinaryCascadePhysics.hh"
69 #include "G4DecayPhysics.hh" 73 #include "G4DecayPhysics.hh"
70 #include "G4NeutronTrackingCut.hh" 74 #include "G4NeutronTrackingCut.hh"
71 #include "G4LossTableManager.hh" 75 #include "G4LossTableManager.hh"
72 #include "G4UnitsTable.hh" 76 #include "G4UnitsTable.hh"
73 #include "G4ProcessManager.hh" 77 #include "G4ProcessManager.hh"
74 #include "G4IonFluctuations.hh" 78 #include "G4IonFluctuations.hh"
75 #include "G4IonParametrisedLossModel.hh" 79 #include "G4IonParametrisedLossModel.hh"
76 #include "G4EmParameters.hh" << 80 #include "G4EmProcessOptions.hh"
77 #include "G4ParallelWorldPhysics.hh" <<
78 #include "G4EmLivermorePhysics.hh" <<
79 #include "G4AutoDelete.hh" <<
80 #include "G4HadronPhysicsQGSP_BIC_AllHP.hh" <<
81 #include "QGSP_BIC_HP.hh" <<
82 #include "QGSP_BIC.hh" <<
83 #include "G4HadronPhysicsQGSP_BERT.hh" <<
84 #include "G4HadronPhysicsQGSP_BERT_HP.hh" <<
85 #include "G4ParallelWorldPhysics.hh" <<
86 // Physics List <<
87 #include "QBBC.hh" <<
88 #include "QGSP_BIC.hh" <<
89 #include "Shielding.hh" <<
90 #include "QGSP_BERT.hh" <<
91 #include "QGSP_BIC_AllHP.hh" <<
92 #include "QGSP_BIC_HP.hh" <<
93 <<
94 <<
95 81
96 ////////////////////////////////////////////// 82 /////////////////////////////////////////////////////////////////////////////
97 HadrontherapyPhysicsList::HadrontherapyPhysics 83 HadrontherapyPhysicsList::HadrontherapyPhysicsList() : G4VModularPhysicsList()
98 { 84 {
99 G4LossTableManager::Instance(); << 85 G4LossTableManager::Instance();
100 defaultCutValue = 1.*mm; << 86 defaultCutValue = 1.*mm;
101 cutForGamma = defaultCutValue; << 87 cutForGamma = defaultCutValue;
102 cutForElectron = defaultCutValue; << 88 cutForElectron = defaultCutValue;
103 cutForPositron = defaultCutValue; << 89 cutForPositron = defaultCutValue;
104 << 90
105 pMessenger = new HadrontherapyPhysicsListM << 91 helIsRegistered = false;
106 SetVerboseLevel(1); << 92 bicIsRegistered = false;
107 decay_List = new G4DecayPhysics(); << 93 biciIsRegistered = false;
108 // Elecromagnetic physics << 94 locIonIonInelasticIsRegistered = false;
109 // << 95 radioactiveDecayIsRegistered = false;
110 emPhysicsList = new G4EmStandardPhysics_op << 96
>> 97 stepMaxProcess = 0;
>> 98
>> 99 pMessenger = new HadrontherapyPhysicsListMessenger(this);
>> 100
>> 101 SetVerboseLevel(1);
>> 102
>> 103 // ****** Definition of some defaults for the physics *****
>> 104 // ****** in case no physics is called by the macro file *****
>> 105 // EM physics
>> 106 emPhysicsList = new G4EmStandardPhysics_option3(1);
>> 107 emName = G4String("emstandard_opt3");
>> 108
>> 109 // Decay physics and all particles
>> 110 decPhysicsList = new G4DecayPhysics();
>> 111 raddecayList = new G4RadioactiveDecayPhysics();
111 } 112 }
112 113
113 ////////////////////////////////////////////// 114 /////////////////////////////////////////////////////////////////////////////
114 HadrontherapyPhysicsList::~HadrontherapyPhysic 115 HadrontherapyPhysicsList::~HadrontherapyPhysicsList()
115 { 116 {
116 delete pMessenger; << 117 delete pMessenger;
117 delete emPhysicsList; << 118 delete emPhysicsList;
118 delete decay_List; << 119 delete decPhysicsList;
119 //delete radioactiveDecay_List; << 120 delete raddecayList;
120 hadronPhys.clear(); << 121
121 for(size_t i=0; i<hadronPhys.size(); i++) << 122 for(size_t i=0; i<hadronPhys.size(); i++) {delete hadronPhys[i];}
122 { <<
123 delete hadronPhys[i]; <<
124 } <<
125 } 123 }
126 124
127 ////////////////////////////////////////////// 125 /////////////////////////////////////////////////////////////////////////////
128 void HadrontherapyPhysicsList::ConstructPartic 126 void HadrontherapyPhysicsList::ConstructParticle()
129 { 127 {
130 decay_List -> ConstructParticle(); << 128 decPhysicsList->ConstructParticle();
131 <<
132 } 129 }
133 130
134 ////////////////////////////////////////////// 131 /////////////////////////////////////////////////////////////////////////////
135 void HadrontherapyPhysicsList::ConstructProces 132 void HadrontherapyPhysicsList::ConstructProcess()
136 { 133 {
137 // Transportation << 134 // transportation
138 // << 135 AddTransportation();
139 AddTransportation(); << 136
140 << 137 // electromagnetic physics list
141 decay_List -> ConstructProcess(); << 138 emPhysicsList->ConstructProcess();
142 emPhysicsList -> ConstructProcess(); << 139 em_config.AddModels();
143 << 140
144 << 141 // decay physics list
145 //em_config.AddModels(); << 142 decPhysicsList->ConstructProcess();
146 << 143 raddecayList->ConstructProcess();
147 // Hadronic physics << 144
148 // << 145 // hadronic physics lists
149 for(size_t i=0; i < hadronPhys.size(); i++ << 146 for(size_t i=0; i<hadronPhys.size(); i++) {
150 { << 147 hadronPhys[i] -> ConstructProcess();
151 hadronPhys[i] -> ConstructProcess(); << 148 }
152 } << 149
153 << 150 // step limitation (as a full process)
154 // step limitation (as a full process) << 151 //
155 // << 152 AddStepMax();
156 AddStepMax(); <<
157 <<
158 //Parallel world sensitivity <<
159 // <<
160 G4ParallelWorldPhysics* pWorld = new G4Par <<
161 pWorld->ConstructProcess(); <<
162 <<
163 return; <<
164 } 153 }
165 154
166 ////////////////////////////////////////////// 155 /////////////////////////////////////////////////////////////////////////////
167 void HadrontherapyPhysicsList::AddPhysicsList( 156 void HadrontherapyPhysicsList::AddPhysicsList(const G4String& name)
168 { 157 {
169 if (verboseLevel>1) { << 158
170 G4cout << "PhysicsList::AddPhysicsList << 159 if (verboseLevel>1) {
171 } << 160 G4cout << "PhysicsList::AddPhysicsList: <" << name << ">" << G4endl;
172 if (name == emName) return; << 161 }
173 << 162 if (name == emName) return;
174 /////////////////////////////////// << 163
175 // ELECTROMAGNETIC MODELS << 164 /////////////////////////////////////////////////////////////////////////////
176 /////////////////////////////////// << 165 // ELECTROMAGNETIC MODELS
177 if (name == "standard_opt4") { << 166 /////////////////////////////////////////////////////////////////////////////
178 emName = name; << 167 if (name == "standard_opt3") {
179 delete emPhysicsList; << 168 emName = name;
180 hadronPhys.clear(); << 169 delete emPhysicsList;
181 emPhysicsList = new G4EmStandardPhysic << 170 emPhysicsList = new G4EmStandardPhysics_option3();
182 G4RunManager::GetRunManager() -> Physi << 171 G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
183 G4cout << "THE FOLLOWING ELECTROMAGNET << 172 G4cout << "THE FOLLOWING ELECTROMAGNETIC PHYSICS LIST HAS BEEN ACTIVATED: G4EmStandardPhysics_option3" << G4endl;
184 << 173
185 ////////////////////////////////////// << 174 } else if (name == "LowE_Livermore") {
186 // ELECTROMAGNETIC + HADRONIC MODELS << 175 emName = name;
187 ////////////////////////////////////// << 176 delete emPhysicsList;
188 << 177 emPhysicsList = new G4EmLivermorePhysics();
189 } else if (name == "HADRONTHERAPY_1") { << 178 G4RunManager::GetRunManager()-> PhysicsHasBeenModified();
190 << 179 G4cout << "THE FOLLOWING ELECTROMAGNETIC PHYSICS LIST HAS BEEN ACTIVATED: G4EmLivermorePhysics" << G4endl;
191 AddPhysicsList("standard_opt4"); << 180
192 hadronPhys.push_back( new G4DecayPhysi << 181 } else if (name == "LowE_Penelope") {
193 hadronPhys.push_back( new G4Radioactiv << 182 emName = name;
194 hadronPhys.push_back( new G4IonBinaryC << 183 delete emPhysicsList;
195 hadronPhys.push_back( new G4EmExtraPhy << 184 emPhysicsList = new G4EmPenelopePhysics();
196 hadronPhys.push_back( new G4HadronElas << 185 G4RunManager::GetRunManager()-> PhysicsHasBeenModified();
197 hadronPhys.push_back( new G4StoppingPh << 186 G4cout << "THE FOLLOWING ELECTROMAGNETIC PHYSICS LIST HAS BEEN ACTIVATED: G4EmPenelopePhysics" << G4endl;
198 hadronPhys.push_back( new G4HadronPhys << 187
199 hadronPhys.push_back( new G4NeutronTra << 188 } else if (name == "local_ion_ion_inelastic") {
200 << 189 hadronPhys.push_back(new LocalIonIonInelasticPhysic());
201 G4cout << "HADRONTHERAPY_1 PHYSICS LIS << 190 locIonIonInelasticIsRegistered = true;
202 } << 191
203 << 192 } else if (name == "QGSP_BIC_EMY") {
204 else if (name == "HADRONTHERAPY_2") { << 193 AddPhysicsList("emstandard_opt3");
205 << 194 hadronPhys.push_back( new HadronPhysicsQGSP_BIC());
206 AddPhysicsList("standard_opt4"); << 195 hadronPhys.push_back( new G4EmExtraPhysics());
207 hadronPhys.push_back( new G4DecayPhysi << 196 hadronPhys.push_back( new G4HadronElasticPhysics());
208 hadronPhys.push_back( new G4Radioactiv << 197 hadronPhys.push_back( new G4QStoppingPhysics());
209 hadronPhys.push_back( new G4IonBinaryC << 198 hadronPhys.push_back( new G4IonBinaryCascadePhysics());
210 hadronPhys.push_back( new G4EmExtraPhy << 199 hadronPhys.push_back( new G4NeutronTrackingCut());
211 hadronPhys.push_back( new G4HadronElas << 200
212 hadronPhys.push_back( new G4StoppingPh << 201 } else {
213 hadronPhys.push_back( new G4HadronPhys << 202 G4cout << "PhysicsList::AddPhysicsList: <" << name << ">"
214 hadronPhys.push_back( new G4NeutronTra << 203 << " is not defined"
215 << 204 << G4endl;
216 G4cout << "HADRONTHERAPY_2 PHYSICS LIS << 205 }
217 <<
218 } <<
219 <<
220 else if (name == "QGSP_BIC"){ <<
221 auto physicsList = new QGSP_BIC; <<
222 G4RunManager::GetRunManager() -> SetUs <<
223 G4RunManager::GetRunManager() -> Physi <<
224 physicsList -> RegisterPhysics(new G4P <<
225 } <<
226 <<
227 else if (name == "QGSP_BERT"){ <<
228 auto physicsList = new QGSP_BERT; <<
229 G4RunManager::GetRunManager() -> SetUs <<
230 G4RunManager::GetRunManager() -> Physi <<
231 physicsList -> RegisterPhysics(new G4P <<
232 } <<
233 <<
234 else if (name == "QGSP_BIC_AllHP"){ <<
235 auto physicsList = new QGSP_BIC_AllHP; <<
236 G4RunManager::GetRunManager() -> SetUs <<
237 G4RunManager::GetRunManager() -> Physi <<
238 physicsList -> RegisterPhysics(new G4P <<
239 } <<
240 <<
241 else if (name == "QGSP_BIC_HP"){ <<
242 auto physicsList = new QGSP_BIC_HP; <<
243 G4RunManager::GetRunManager() -> SetUs <<
244 G4RunManager::GetRunManager() -> Physi <<
245 physicsList -> RegisterPhysics(new G4P <<
246 } <<
247 <<
248 else if (name == "Shielding"){ <<
249 auto physicsList = new Shielding; <<
250 G4RunManager::GetRunManager() -> SetUs <<
251 G4RunManager::GetRunManager() -> Physi <<
252 physicsList -> RegisterPhysics(new G4P <<
253 } <<
254 <<
255 else if (name == "QBBC"){ <<
256 auto physicsList = new QBBC; <<
257 G4RunManager::GetRunManager() -> SetUs <<
258 G4RunManager::GetRunManager() -> Physi <<
259 physicsList -> RegisterPhysics(new G4P <<
260 } <<
261 <<
262 else { <<
263 G4cout << "PhysicsList::AddPhysicsList <<
264 << " is not defined" <<
265 << G4endl; <<
266 } <<
267 <<
268 } 206 }
269 207
270 ////////////////////////////////////////////// 208 /////////////////////////////////////////////////////////////////////////////
271 void HadrontherapyPhysicsList::AddStepMax() 209 void HadrontherapyPhysicsList::AddStepMax()
272 { 210 {
273 // Step limitation seen as a process << 211 // Step limitation seen as a process
274 // This process must exist in all threads. << 212 stepMaxProcess = new HadrontherapyStepMax();
275 // << 213
276 HadrontherapyStepMax* stepMaxProcess = ne << 214 theParticleIterator->reset();
277 << 215 while ((*theParticleIterator)()){
278 << 216 G4ParticleDefinition* particle = theParticleIterator->value();
279 auto particleIterator = GetParticleIterato << 217 G4ProcessManager* pmanager = particle->GetProcessManager();
280 particleIterator->reset(); << 218
281 while ((*particleIterator)()){ << 219 if (stepMaxProcess->IsApplicable(*particle) && pmanager)
282 G4ParticleDefinition* particle = parti << 220 {
283 G4ProcessManager* pmanager = particle- << 221 pmanager ->AddDiscreteProcess(stepMaxProcess);
284 << 222 }
285 if (stepMaxProcess->IsApplicable(*part << 223 }
286 { << 224 }
287 pmanager ->AddDiscreteProcess(step << 225
288 } << 226 /////////////////////////////////////////////////////////////////////////////
289 } << 227 void HadrontherapyPhysicsList::SetCuts()
>> 228 {
>> 229
>> 230 if (verboseLevel >0){
>> 231 G4cout << "PhysicsList::SetCuts:";
>> 232 G4cout << "CutLength : " << G4BestUnit(defaultCutValue,"Length") << G4endl;
>> 233 }
>> 234
>> 235 // set cut values for gamma at first and for e- second and next for e+,
>> 236 // because some processes for e+/e- need cut values for gamma
>> 237 SetCutValue(cutForGamma, "gamma");
>> 238 SetCutValue(cutForElectron, "e-");
>> 239 SetCutValue(cutForPositron, "e+");
>> 240
>> 241 // Set cuts for detector
>> 242 SetDetectorCut(defaultCutValue);
>> 243 if (verboseLevel>0) DumpCutValuesTable();
>> 244 }
>> 245
>> 246 /////////////////////////////////////////////////////////////////////////////
>> 247 void HadrontherapyPhysicsList::SetCutForGamma(G4double cut)
>> 248 {
>> 249 cutForGamma = cut;
>> 250 SetParticleCuts(cutForGamma, G4Gamma::Gamma());
>> 251 }
>> 252
>> 253 /////////////////////////////////////////////////////////////////////////////
>> 254 void HadrontherapyPhysicsList::SetCutForElectron(G4double cut)
>> 255 {
>> 256 cutForElectron = cut;
>> 257 SetParticleCuts(cutForElectron, G4Electron::Electron());
290 } 258 }
>> 259
>> 260 /////////////////////////////////////////////////////////////////////////////
>> 261 void HadrontherapyPhysicsList::SetCutForPositron(G4double cut)
>> 262 {
>> 263 cutForPositron = cut;
>> 264 SetParticleCuts(cutForPositron, G4Positron::Positron());
>> 265 }
>> 266
>> 267 void HadrontherapyPhysicsList::SetDetectorCut(G4double cut)
>> 268 {
>> 269 G4String regionName = "DetectorLog";
>> 270 G4Region* region = G4RegionStore::GetInstance()->GetRegion(regionName);
>> 271
>> 272 G4ProductionCuts* cuts = new G4ProductionCuts ;
>> 273 cuts -> SetProductionCut(cut,G4ProductionCuts::GetIndex("gamma"));
>> 274 cuts -> SetProductionCut(cut,G4ProductionCuts::GetIndex("e-"));
>> 275 cuts -> SetProductionCut(cut,G4ProductionCuts::GetIndex("e+"));
>> 276 region -> SetProductionCuts(cuts);
>> 277 }
>> 278
291 279