| 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):
>> 46 // export PHYSLIST=QGSP_BIC_EMY
>> 47 // User must considered that, in this second case, radioactive processes are not activated
>> 48
49 #include "G4SystemOfUnits.hh" 49 #include "G4SystemOfUnits.hh"
50 #include "G4RunManager.hh" 50 #include "G4RunManager.hh"
51 #include "G4Region.hh" 51 #include "G4Region.hh"
52 #include "G4RegionStore.hh" 52 #include "G4RegionStore.hh"
53 #include "HadrontherapyPhysicsList.hh" 53 #include "HadrontherapyPhysicsList.hh"
54 #include "HadrontherapyPhysicsListMessenger.hh 54 #include "HadrontherapyPhysicsListMessenger.hh"
55 #include "HadrontherapyStepMax.hh" 55 #include "HadrontherapyStepMax.hh"
56 #include "G4PhysListFactory.hh" 56 #include "G4PhysListFactory.hh"
57 #include "G4VPhysicsConstructor.hh" 57 #include "G4VPhysicsConstructor.hh"
58 #include "G4HadronPhysicsQGSP_BIC_HP.hh" << 58
59 #include "G4HadronPhysicsQGSP_BIC.hh" << 59 // Local physic directly implemented in the Hadronthrapy directory
60 #include "G4EmStandardPhysics_option4.hh" << 60 #include "LocalIonIonInelasticPhysic.hh" // Physic dedicated to the ion-ion inelastic processes
61 #include "G4EmStandardPhysics.hh" << 61
>> 62 // Physic lists (contained inside the Geant4 source code, in the 'physicslists folder')
>> 63 #include "HadronPhysicsQGSP_BIC.hh"
>> 64 #include "G4EmStandardPhysics_option3.hh"
>> 65 #include "G4EmLivermorePhysics.hh"
>> 66 #include "G4EmPenelopePhysics.hh"
62 #include "G4EmExtraPhysics.hh" 67 #include "G4EmExtraPhysics.hh"
63 #include "G4StoppingPhysics.hh" 68 #include "G4StoppingPhysics.hh"
64 #include "G4DecayPhysics.hh" 69 #include "G4DecayPhysics.hh"
65 #include "G4HadronElasticPhysics.hh" 70 #include "G4HadronElasticPhysics.hh"
66 #include "G4HadronElasticPhysicsHP.hh" 71 #include "G4HadronElasticPhysicsHP.hh"
67 #include "G4RadioactiveDecayPhysics.hh" 72 #include "G4RadioactiveDecayPhysics.hh"
68 #include "G4IonBinaryCascadePhysics.hh" 73 #include "G4IonBinaryCascadePhysics.hh"
69 #include "G4DecayPhysics.hh" 74 #include "G4DecayPhysics.hh"
70 #include "G4NeutronTrackingCut.hh" 75 #include "G4NeutronTrackingCut.hh"
71 #include "G4LossTableManager.hh" 76 #include "G4LossTableManager.hh"
72 #include "G4UnitsTable.hh" 77 #include "G4UnitsTable.hh"
73 #include "G4ProcessManager.hh" 78 #include "G4ProcessManager.hh"
74 #include "G4IonFluctuations.hh" 79 #include "G4IonFluctuations.hh"
75 #include "G4IonParametrisedLossModel.hh" 80 #include "G4IonParametrisedLossModel.hh"
76 #include "G4EmParameters.hh" << 81 #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 82
96 ////////////////////////////////////////////// 83 /////////////////////////////////////////////////////////////////////////////
97 HadrontherapyPhysicsList::HadrontherapyPhysics 84 HadrontherapyPhysicsList::HadrontherapyPhysicsList() : G4VModularPhysicsList()
98 { 85 {
99 G4LossTableManager::Instance(); << 86 G4LossTableManager::Instance();
100 defaultCutValue = 1.*mm; << 87 defaultCutValue = 1.*mm;
101 cutForGamma = defaultCutValue; << 88 cutForGamma = defaultCutValue;
102 cutForElectron = defaultCutValue; << 89 cutForElectron = defaultCutValue;
103 cutForPositron = defaultCutValue; << 90 cutForPositron = defaultCutValue;
104 << 91
105 pMessenger = new HadrontherapyPhysicsListM << 92 helIsRegistered = false;
106 SetVerboseLevel(1); << 93 bicIsRegistered = false;
107 decay_List = new G4DecayPhysics(); << 94 biciIsRegistered = false;
108 // Elecromagnetic physics << 95 locIonIonInelasticIsRegistered = false;
109 // << 96 radioactiveDecayIsRegistered = false;
110 emPhysicsList = new G4EmStandardPhysics_op << 97
>> 98 stepMaxProcess = 0;
>> 99
>> 100 pMessenger = new HadrontherapyPhysicsListMessenger(this);
>> 101
>> 102 SetVerboseLevel(1);
>> 103
>> 104 // ****** Definition of some defaults for the physics *****
>> 105 // ****** in case no physics is called by the macro file *****
>> 106 // EM physics
>> 107 emPhysicsList = new G4EmStandardPhysics_option3(1);
>> 108 emName = G4String("emstandard_opt3");
>> 109
>> 110 // Decay physics and all particles
>> 111 decPhysicsList = new G4DecayPhysics();
>> 112 raddecayList = new G4RadioactiveDecayPhysics();
111 } 113 }
112 114
113 ////////////////////////////////////////////// 115 /////////////////////////////////////////////////////////////////////////////
114 HadrontherapyPhysicsList::~HadrontherapyPhysic 116 HadrontherapyPhysicsList::~HadrontherapyPhysicsList()
115 { 117 {
116 delete pMessenger; << 118 delete pMessenger;
117 delete emPhysicsList; << 119 delete emPhysicsList;
118 delete decay_List; << 120 delete decPhysicsList;
119 //delete radioactiveDecay_List; << 121 delete raddecayList;
120 hadronPhys.clear(); << 122
121 for(size_t i=0; i<hadronPhys.size(); i++) << 123 for(size_t i=0; i<hadronPhys.size(); i++) {delete hadronPhys[i];}
122 { <<
123 delete hadronPhys[i]; <<
124 } <<
125 } 124 }
126 125
127 ////////////////////////////////////////////// 126 /////////////////////////////////////////////////////////////////////////////
128 void HadrontherapyPhysicsList::ConstructPartic 127 void HadrontherapyPhysicsList::ConstructParticle()
129 { 128 {
130 decay_List -> ConstructParticle(); << 129 decPhysicsList->ConstructParticle();
131 <<
132 } 130 }
133 131
134 ////////////////////////////////////////////// 132 /////////////////////////////////////////////////////////////////////////////
135 void HadrontherapyPhysicsList::ConstructProces 133 void HadrontherapyPhysicsList::ConstructProcess()
136 { 134 {
137 // Transportation << 135 // transportation
138 // << 136 AddTransportation();
139 AddTransportation(); << 137
140 << 138 // electromagnetic physics list
141 decay_List -> ConstructProcess(); << 139 emPhysicsList->ConstructProcess();
142 emPhysicsList -> ConstructProcess(); << 140 em_config.AddModels();
143 << 141
144 << 142 // decay physics list
145 //em_config.AddModels(); << 143 decPhysicsList->ConstructProcess();
146 << 144 raddecayList->ConstructProcess();
147 // Hadronic physics << 145
148 // << 146 // hadronic physics lists
149 for(size_t i=0; i < hadronPhys.size(); i++ << 147 for(size_t i=0; i<hadronPhys.size(); i++) {
150 { << 148 hadronPhys[i] -> ConstructProcess();
151 hadronPhys[i] -> ConstructProcess(); << 149 }
152 } << 150
153 << 151 // step limitation (as a full process)
154 // step limitation (as a full process) << 152 //
155 // << 153 AddStepMax();
156 AddStepMax(); <<
157 <<
158 //Parallel world sensitivity <<
159 // <<
160 G4ParallelWorldPhysics* pWorld = new G4Par <<
161 pWorld->ConstructProcess(); <<
162 <<
163 return; <<
164 } 154 }
165 155
166 ////////////////////////////////////////////// 156 /////////////////////////////////////////////////////////////////////////////
167 void HadrontherapyPhysicsList::AddPhysicsList( 157 void HadrontherapyPhysicsList::AddPhysicsList(const G4String& name)
168 { 158 {
169 if (verboseLevel>1) { << 159
170 G4cout << "PhysicsList::AddPhysicsList << 160 if (verboseLevel>1) {
171 } << 161 G4cout << "PhysicsList::AddPhysicsList: <" << name << ">" << G4endl;
172 if (name == emName) return; << 162 }
173 << 163 if (name == emName) return;
174 /////////////////////////////////// << 164
175 // ELECTROMAGNETIC MODELS << 165 /////////////////////////////////////////////////////////////////////////////
176 /////////////////////////////////// << 166 // ELECTROMAGNETIC MODELS
177 if (name == "standard_opt4") { << 167 /////////////////////////////////////////////////////////////////////////////
178 emName = name; << 168 if (name == "standard_opt3") {
179 delete emPhysicsList; << 169 emName = name;
180 hadronPhys.clear(); << 170 delete emPhysicsList;
181 emPhysicsList = new G4EmStandardPhysic << 171 emPhysicsList = new G4EmStandardPhysics_option3();
182 G4RunManager::GetRunManager() -> Physi << 172 G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
183 G4cout << "THE FOLLOWING ELECTROMAGNET << 173 G4cout << "THE FOLLOWING ELECTROMAGNETIC PHYSICS LIST HAS BEEN ACTIVATED: G4EmStandardPhysics_option3" << G4endl;
184 << 174
185 ////////////////////////////////////// << 175 } else if (name == "LowE_Livermore") {
186 // ELECTROMAGNETIC + HADRONIC MODELS << 176 emName = name;
187 ////////////////////////////////////// << 177 delete emPhysicsList;
188 << 178 emPhysicsList = new G4EmLivermorePhysics();
189 } else if (name == "HADRONTHERAPY_1") { << 179 G4RunManager::GetRunManager()-> PhysicsHasBeenModified();
190 << 180 G4cout << "THE FOLLOWING ELECTROMAGNETIC PHYSICS LIST HAS BEEN ACTIVATED: G4EmLivermorePhysics" << G4endl;
191 AddPhysicsList("standard_opt4"); << 181
192 hadronPhys.push_back( new G4DecayPhysi << 182 } else if (name == "LowE_Penelope") {
193 hadronPhys.push_back( new G4Radioactiv << 183 emName = name;
194 hadronPhys.push_back( new G4IonBinaryC << 184 delete emPhysicsList;
195 hadronPhys.push_back( new G4EmExtraPhy << 185 emPhysicsList = new G4EmPenelopePhysics();
196 hadronPhys.push_back( new G4HadronElas << 186 G4RunManager::GetRunManager()-> PhysicsHasBeenModified();
197 hadronPhys.push_back( new G4StoppingPh << 187 G4cout << "THE FOLLOWING ELECTROMAGNETIC PHYSICS LIST HAS BEEN ACTIVATED: G4EmPenelopePhysics" << G4endl;
198 hadronPhys.push_back( new G4HadronPhys << 188
199 hadronPhys.push_back( new G4NeutronTra << 189 } else if (name == "local_ion_ion_inelastic") {
200 << 190 hadronPhys.push_back(new LocalIonIonInelasticPhysic());
201 G4cout << "HADRONTHERAPY_1 PHYSICS LIS << 191 locIonIonInelasticIsRegistered = true;
202 } << 192
203 << 193 } else if (name == "QGSP_BIC_EMY") {
204 else if (name == "HADRONTHERAPY_2") { << 194 AddPhysicsList("emstandard_opt3");
205 << 195 hadronPhys.push_back( new HadronPhysicsQGSP_BIC());
206 AddPhysicsList("standard_opt4"); << 196 hadronPhys.push_back( new G4EmExtraPhysics());
207 hadronPhys.push_back( new G4DecayPhysi << 197 hadronPhys.push_back( new G4HadronElasticPhysics());
208 hadronPhys.push_back( new G4Radioactiv << 198 hadronPhys.push_back( new G4StoppingPhysics());
209 hadronPhys.push_back( new G4IonBinaryC << 199 hadronPhys.push_back( new G4IonBinaryCascadePhysics());
210 hadronPhys.push_back( new G4EmExtraPhy << 200 hadronPhys.push_back( new G4NeutronTrackingCut());
211 hadronPhys.push_back( new G4HadronElas << 201
212 hadronPhys.push_back( new G4StoppingPh << 202 } else {
213 hadronPhys.push_back( new G4HadronPhys << 203 G4cout << "PhysicsList::AddPhysicsList: <" << name << ">"
214 hadronPhys.push_back( new G4NeutronTra << 204 << " is not defined"
215 << 205 << G4endl;
216 G4cout << "HADRONTHERAPY_2 PHYSICS LIS << 206 }
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 } 207 }
269 208
270 ////////////////////////////////////////////// 209 /////////////////////////////////////////////////////////////////////////////
271 void HadrontherapyPhysicsList::AddStepMax() 210 void HadrontherapyPhysicsList::AddStepMax()
272 { 211 {
273 // Step limitation seen as a process << 212 // Step limitation seen as a process
274 // This process must exist in all threads. << 213 stepMaxProcess = new HadrontherapyStepMax();
275 // << 214
276 HadrontherapyStepMax* stepMaxProcess = ne << 215 theParticleIterator->reset();
277 << 216 while ((*theParticleIterator)()){
278 << 217 G4ParticleDefinition* particle = theParticleIterator->value();
279 auto particleIterator = GetParticleIterato << 218 G4ProcessManager* pmanager = particle->GetProcessManager();
280 particleIterator->reset(); << 219
281 while ((*particleIterator)()){ << 220 if (stepMaxProcess->IsApplicable(*particle) && pmanager)
282 G4ParticleDefinition* particle = parti << 221 {
283 G4ProcessManager* pmanager = particle- << 222 pmanager ->AddDiscreteProcess(stepMaxProcess);
284 << 223 }
285 if (stepMaxProcess->IsApplicable(*part << 224 }
286 { << 225 }
287 pmanager ->AddDiscreteProcess(step << 226
288 } << 227 /////////////////////////////////////////////////////////////////////////////
289 } << 228 void HadrontherapyPhysicsList::SetCuts()
>> 229 {
>> 230
>> 231 if (verboseLevel >0){
>> 232 G4cout << "PhysicsList::SetCuts:";
>> 233 G4cout << "CutLength : " << G4BestUnit(defaultCutValue,"Length") << G4endl;
>> 234 }
>> 235
>> 236 // set cut values for gamma at first and for e- second and next for e+,
>> 237 // because some processes for e+/e- need cut values for gamma
>> 238 SetCutValue(cutForGamma, "gamma");
>> 239 SetCutValue(cutForElectron, "e-");
>> 240 SetCutValue(cutForPositron, "e+");
>> 241
>> 242 // Set cuts for detector
>> 243 SetDetectorCut(defaultCutValue);
>> 244 if (verboseLevel>0) DumpCutValuesTable();
>> 245 }
>> 246
>> 247 /////////////////////////////////////////////////////////////////////////////
>> 248 void HadrontherapyPhysicsList::SetCutForGamma(G4double cut)
>> 249 {
>> 250 cutForGamma = cut;
>> 251 SetParticleCuts(cutForGamma, G4Gamma::Gamma());
>> 252 }
>> 253
>> 254 /////////////////////////////////////////////////////////////////////////////
>> 255 void HadrontherapyPhysicsList::SetCutForElectron(G4double cut)
>> 256 {
>> 257 cutForElectron = cut;
>> 258 SetParticleCuts(cutForElectron, G4Electron::Electron());
290 } 259 }
>> 260
>> 261 /////////////////////////////////////////////////////////////////////////////
>> 262 void HadrontherapyPhysicsList::SetCutForPositron(G4double cut)
>> 263 {
>> 264 cutForPositron = cut;
>> 265 SetParticleCuts(cutForPositron, G4Positron::Positron());
>> 266 }
>> 267
>> 268 void HadrontherapyPhysicsList::SetDetectorCut(G4double cut)
>> 269 {
>> 270 G4String regionName = "DetectorLog";
>> 271 G4Region* region = G4RegionStore::GetInstance()->GetRegion(regionName);
>> 272
>> 273 G4ProductionCuts* cuts = new G4ProductionCuts ;
>> 274 cuts -> SetProductionCut(cut,G4ProductionCuts::GetIndex("gamma"));
>> 275 cuts -> SetProductionCut(cut,G4ProductionCuts::GetIndex("e-"));
>> 276 cuts -> SetProductionCut(cut,G4ProductionCuts::GetIndex("e+"));
>> 277 region -> SetProductionCuts(cuts);
>> 278 }
>> 279
291 280