| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 /// \file PhysicsList.cc
28 /// \brief Implementation of the PhysicsList class.
29
30 #include "PhysicsList.hh"
31
32 #include "G4SystemOfUnits.hh"
33 #include "G4UIcmdWithAnInteger.hh"
34
35 // Geant4-DNA MODELS
36
37 #include "G4DNAAttachment.hh"
38 #include "G4DNABornExcitationModel.hh"
39 #include "G4DNABornIonisationModel.hh"
40 #include "G4DNAChampionElasticModel.hh"
41 #include "G4DNAChargeDecrease.hh"
42 #include "G4DNAChargeIncrease.hh"
43 #include "G4DNADingfelderChargeDecreaseModel.hh"
44 #include "G4DNADingfelderChargeIncreaseModel.hh"
45 #include "G4DNAElastic.hh"
46 #include "G4DNAElectronSolvation.hh"
47 #include "G4DNAExcitation.hh"
48 #include "G4DNAIonisation.hh"
49 #include "G4DNAMeltonAttachmentModel.hh"
50 #include "G4DNAMillerGreenExcitationModel.hh"
51 #include "G4DNARuddIonisationModel.hh"
52 #include "G4DNASancheExcitationModel.hh"
53 #include "G4DNAScreenedRutherfordElasticModel.hh"
54 #include "G4DNAVibExcitation.hh"
55
56 //
57
58 #include "SplitProcess.hh"
59
60 #include "G4BetheBlochIonGasModel.hh"
61 #include "G4BraggIonGasModel.hh"
62 #include "G4DummyModel.hh"
63 #include "G4ElectronCapture.hh"
64 #include "G4EmConfigurator.hh"
65 #include "G4IonFluctuations.hh"
66 #include "G4LossTableManager.hh"
67 #include "G4MollerBhabhaModel.hh"
68 #include "G4UniversalFluctuation.hh"
69 #include "G4UrbanMscModel.hh"
70 #include "G4VEmModel.hh"
71 #include "G4eIonisation.hh"
72 #include "G4eMultipleScattering.hh"
73 #include "G4hIonisation.hh"
74 #include "G4hMultipleScattering.hh"
75 #include "G4ionIonisation.hh"
76
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78
79 PhysicsList::PhysicsList() : G4VUserPhysicsList(), G4UImessenger()
80
81 {
82 fpNumberOfSplit = new G4UIcmdWithAnInteger("/vrt/numberOfSplit", this);
83 fpNumberOfSplit->SetGuidance("Number of split for particle splitting");
84 fpNumberOfSplit->SetParameterName("NumberOfSplit", false);
85 fpNumberOfSplit->SetDefaultValue(1);
86
87 defaultCutValue = 1 * micrometer;
88 fCutForGamma = defaultCutValue;
89 fCutForElectron = defaultCutValue;
90 fCutForPositron = defaultCutValue;
91 fCutForProton = defaultCutValue;
92
93 SetVerboseLevel(1);
94
95 // fNumberOfSplit = 1;
96 }
97
98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
99
100 PhysicsList::~PhysicsList() {}
101
102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
103 void PhysicsList::SetNewValue(G4UIcommand* command, G4String newValue)
104 {
105 if (command == fpNumberOfSplit) fNumberOfSplit = fpNumberOfSplit->GetNewIntValue(newValue);
106 }
107
108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
109
110 void PhysicsList::ConstructParticle()
111 {
112 ConstructBosons();
113 ConstructLeptons();
114 ConstructBarions();
115 }
116
117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
118
119 void PhysicsList::ConstructBosons()
120 {
121 // gamma
122 G4Gamma::GammaDefinition();
123 }
124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
125
126 void PhysicsList::ConstructLeptons()
127 {
128 // leptons
129 G4Electron::ElectronDefinition();
130 G4Positron::PositronDefinition();
131 }
132
133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
134
135 // DNA
136 #include "G4DNAGenericIonsManager.hh"
137 // ENDDNA
138
139 void PhysicsList::ConstructBarions()
140 {
141 // baryons
142 G4Proton::ProtonDefinition();
143 G4GenericIon::GenericIonDefinition();
144
145 // Geant4 DNA new particles
146 G4DNAGenericIonsManager* genericIonsManager;
147 genericIonsManager = G4DNAGenericIonsManager::Instance();
148 genericIonsManager->GetIon("alpha++");
149 genericIonsManager->GetIon("alpha+");
150 genericIonsManager->GetIon("helium");
151 genericIonsManager->GetIon("hydrogen");
152 }
153
154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
155
156 void PhysicsList::ConstructProcess()
157 {
158 AddTransportation();
159 ConstructEM();
160 ConstructGeneral();
161 }
162
163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
164
165 void PhysicsList::ConstructEM()
166 {
167 auto particleIterator = GetParticleIterator();
168 particleIterator->reset();
169 G4bool splitProcessActive = false;
170 // Only if particle split is larger than 1, VRT is set
171 if (fNumberOfSplit > 1) splitProcessActive = true;
172
173 while ((*particleIterator)()) {
174 G4ParticleDefinition* particle = particleIterator->value();
175 G4ProcessManager* pmanager = particle->GetProcessManager();
176 G4String particleName = particle->GetParticleName();
177
178 // *********************************
179 // 1) Processes for the World region
180 // *********************************
181
182 if (particleName == "e-") {
183 // STANDARD msc is active in the world
184 G4eMultipleScattering* msc = new G4eMultipleScattering();
185 msc->AddEmModel(1, new G4UrbanMscModel());
186 pmanager->AddProcess(msc, -1, 1, -1);
187
188 // STANDARD ionisation is active in the world
189 G4eIonisation* eion = new G4eIonisation();
190 eion->SetEmModel(new G4MollerBhabhaModel());
191 pmanager->AddProcess(eion, -1, 2, 2);
192
193 // DNA elastic is not active in the world
194 G4DNAElastic* theDNAElasticProcess = new G4DNAElastic("e-_G4DNAElastic");
195 theDNAElasticProcess->SetEmModel(new G4DummyModel());
196 pmanager->AddDiscreteProcess(theDNAElasticProcess);
197
198 // DNA excitation is not active in the world
199 G4DNAExcitation* dnaex = new G4DNAExcitation("e-_G4DNAExcitation");
200 dnaex->SetEmModel(new G4DummyModel());
201 pmanager->AddDiscreteProcess(dnaex);
202
203 // DNA ionisation is not active in the world
204 G4DNAIonisation* dnaioni = new G4DNAIonisation("e-_G4DNAIonisation");
205 dnaioni->SetEmModel(new G4DummyModel());
206
207 // Variance reduction is set here
208 if (splitProcessActive) {
209 G4String splitRegion = "Target";
210 G4cout << "- Split for e-e created in e-_G4DNAIonisation process activated"
211 << "with split number " << fNumberOfSplit << " in region " << splitRegion << "-- "
212 << G4endl;
213 SplitProcess* splitProcess = new SplitProcess(splitRegion, fNumberOfSplit);
214 splitProcess->RegisterProcess(dnaioni);
215
216 pmanager->AddDiscreteProcess(splitProcess);
217 }
218 else {
219 pmanager->AddDiscreteProcess(dnaioni);
220 }
221
222 // DNA attachment is not active in the world
223 G4DNAAttachment* dnaatt = new G4DNAAttachment("e-_G4DNAAttachment");
224 dnaatt->SetEmModel(new G4DummyModel());
225 pmanager->AddDiscreteProcess(dnaatt);
226
227 // DNA vib. excitation is not active in the world
228 G4DNAVibExcitation* dnavib = new G4DNAVibExcitation("e-_G4DNAVibExcitation");
229 dnavib->SetEmModel(new G4DummyModel());
230 pmanager->AddDiscreteProcess(dnavib);
231
232 // THE FOLLOWING PROCESS WILL KILL ALL ELECTRONS BELOW A
233 // SELECTED ENERY THRESHOLD
234 // Capture of low-energy e-
235 G4ElectronCapture* ecap = new G4ElectronCapture("Target", 7.4 * eV);
236 pmanager->AddDiscreteProcess(ecap);
237 // 7.4 eV is compatible with the validity range of the Champion's model
238 }
239 else if (particleName == "proton") {
240 // STANDARD msc is active in the world
241 G4hMultipleScattering* msc = new G4hMultipleScattering();
242 msc->AddEmModel(1, new G4UrbanMscModel());
243 pmanager->AddProcess(msc, -1, 1, -1);
244
245 // STANDARD ionisation is active in the world
246 G4hIonisation* hion = new G4hIonisation();
247 hion->SetEmModel(new G4BraggIonGasModel());
248 hion->SetEmModel(new G4BetheBlochIonGasModel());
249 pmanager->AddProcess(hion, -1, 2, 2);
250
251 // DNA excitation is not active in the world
252 G4DNAExcitation* dnaex = new G4DNAExcitation("proton_G4DNAExcitation");
253 dnaex->SetEmModel(new G4DummyModel());
254 dnaex->SetEmModel(new G4DummyModel());
255 pmanager->AddDiscreteProcess(dnaex);
256
257 // DNA ionisation is not active in the world
258 G4DNAIonisation* dnaioni = new G4DNAIonisation("proton_G4DNAIonisation");
259 dnaioni->SetEmModel(new G4DummyModel());
260 dnaioni->SetEmModel(new G4DummyModel());
261 pmanager->AddDiscreteProcess(dnaioni);
262
263 // DNA charge decrease is ACTIVE in the world since
264 // no corresponding STANDARD process exist
265 pmanager->AddDiscreteProcess(new G4DNAChargeDecrease("proton_G4DNAChargeDecrease"));
266 }
267 else if (particleName == "hydrogen") {
268 // DNA processes are ACTIVE in the world since
269 // no corresponding STANDARD processes exist
270 pmanager->AddDiscreteProcess(new G4DNAIonisation("hydrogen_G4DNAIonisation"));
271 pmanager->AddDiscreteProcess(new G4DNAExcitation("hydrogen_G4DNAExcitation"));
272 pmanager->AddDiscreteProcess(new G4DNAChargeIncrease("hydrogen_G4DNAChargeIncrease"));
273 }
274 else if (particleName == "GenericIon") {
275 // THIS IS NEEDED FOR STANDARD ALPHA G4ionIonisation PROCESS
276
277 // STANDARD msc is active in the world
278 G4hMultipleScattering* msc = new G4hMultipleScattering();
279 msc->AddEmModel(1, new G4UrbanMscModel());
280 pmanager->AddProcess(msc, -1, 1, -1);
281
282 // STANDARD ionisation is active in the world
283 G4ionIonisation* hion = new G4ionIonisation();
284 hion->SetEmModel(new G4BraggIonGasModel());
285 hion->SetEmModel(new G4BetheBlochIonGasModel());
286 pmanager->AddProcess(hion, -1, 2, 2);
287 }
288 else if (particleName == "alpha") {
289 // STANDARD msc is active in the world
290 G4hMultipleScattering* msc = new G4hMultipleScattering();
291 msc->AddEmModel(1, new G4UrbanMscModel());
292 pmanager->AddProcess(msc, -1, 1, -1);
293
294 // STANDARD ionisation is active in the world
295 G4ionIonisation* hion = new G4ionIonisation();
296 hion->SetEmModel(new G4BraggIonGasModel());
297 hion->SetEmModel(new G4BetheBlochIonGasModel());
298 pmanager->AddProcess(hion, -1, 2, 2);
299
300 // DNA excitation is not active in the world
301 G4DNAExcitation* dnaex = new G4DNAExcitation("alpha_G4DNAExcitation");
302 dnaex->SetEmModel(new G4DummyModel());
303 pmanager->AddDiscreteProcess(dnaex);
304
305 // DNA ionisation is not active in the world
306 G4DNAIonisation* dnaioni = new G4DNAIonisation("alpha_G4DNAIonisation");
307 dnaioni->SetEmModel(new G4DummyModel());
308 pmanager->AddDiscreteProcess(dnaioni);
309
310 // DNA charge decrease is ACTIVE in the world since no
311 // corresponding STANDARD process exist
312 pmanager->AddDiscreteProcess(new G4DNAChargeDecrease("alpha_G4DNAChargeDecrease"));
313 }
314 else if (particleName == "alpha+") {
315 // DNA processes are ACTIVE in the world since no
316 // corresponding STANDARD processes exist
317 pmanager->AddDiscreteProcess(new G4DNAExcitation("alpha+_G4DNAExcitation"));
318 pmanager->AddDiscreteProcess(new G4DNAIonisation("alpha+_G4DNAIonisation"));
319 pmanager->AddDiscreteProcess(new G4DNAChargeDecrease("alpha+_G4DNAChargeDecrease"));
320 pmanager->AddDiscreteProcess(new G4DNAChargeIncrease("alpha+_G4DNAChargeIncrease"));
321 }
322 else if (particleName == "helium") {
323 // DNA processes are ACTIVE in the world since no
324 // corresponding STANDARD processes exist
325 pmanager->AddDiscreteProcess(new G4DNAExcitation("helium_G4DNAExcitation"));
326 pmanager->AddDiscreteProcess(new G4DNAIonisation("helium_G4DNAIonisation"));
327 pmanager->AddDiscreteProcess(new G4DNAChargeIncrease("helium_G4DNAChargeIncrease"));
328 }
329 }
330
331 // **************************************
332 // 2) Define processes for Target region
333 // **************************************
334
335 // STANDARD EM processes should be inactivated when
336 // corresponding DNA processes are used
337 // - STANDARD EM e- processes are inactivated below 1 MeV
338 // - STANDARD EM proton & alpha processes are inactivated below
339 // standEnergyLimit
340 G4double standEnergyLimit = 9.9 * MeV;
341 //
342
343 G4double massFactor = 1.0079 / 4.0026;
344 G4EmConfigurator* em_config = G4LossTableManager::Instance()->EmConfigurator();
345
346 G4VEmModel* mod;
347
348 // *** e-
349
350 // ---> STANDARD EM processes are inactivated below 1 MeV
351
352 mod = new G4UrbanMscModel();
353 mod->SetActivationLowEnergyLimit(1. * MeV);
354 em_config->SetExtraEmModel("e-", "msc", mod, "Target");
355
356 mod = new G4MollerBhabhaModel();
357 mod->SetActivationLowEnergyLimit(1 * MeV);
358 em_config->SetExtraEmModel("e-", "eIoni", mod, "Target", 0.0, 100 * TeV,
359 new G4UniversalFluctuation());
360
361 // ---> DNA processes activated
362
363 mod = new G4DNAChampionElasticModel();
364 em_config->SetExtraEmModel("e-", "e-_G4DNAElastic", mod, "Target", 7.4 * eV, 1. * MeV);
365
366 mod = new G4DNABornIonisationModel();
367 em_config->SetExtraEmModel("e-", "e-_G4DNAIonisation", mod, "Target", 11. * eV, 1. * MeV);
368 // Note: valid from 11 eV to 0.999.. MeV then switch to std models at
369 // higher energies ; same for other models
370
371 mod = new G4DNABornExcitationModel();
372 em_config->SetExtraEmModel("e-", "e-_G4DNAExcitation", mod, "Target", 9. * eV, 1. * MeV);
373
374 mod = new G4DNAMeltonAttachmentModel();
375 em_config->SetExtraEmModel("e-", "e-_G4DNAAttachment", mod, "Target", 4. * eV, 13. * eV);
376
377 mod = new G4DNASancheExcitationModel();
378 em_config->SetExtraEmModel("e-", "e-_G4DNAVibExcitation", mod, "Target", 2. * eV, 100. * eV);
379
380 // *** proton
381
382 // ---> STANDARD EM processes inactivated below standEnergyLimit
383
384 // STANDARD msc is still active
385 // Inactivate following STANDARD processes
386
387 mod = new G4BraggIonGasModel();
388 mod->SetActivationLowEnergyLimit(standEnergyLimit);
389 em_config->SetExtraEmModel("proton", "hIoni", mod, "Target", 0.0, 2 * MeV,
390 new G4IonFluctuations());
391
392 mod = new G4BetheBlochIonGasModel();
393 mod->SetActivationLowEnergyLimit(standEnergyLimit);
394 em_config->SetExtraEmModel("proton", "hIoni", mod, "Target", 2 * MeV, 100 * TeV,
395 new G4UniversalFluctuation());
396
397 // ---> DNA processes activated
398
399 mod = new G4DNARuddIonisationModel();
400 em_config->SetExtraEmModel("proton", "proton_G4DNAIonisation", mod, "Target", 0.0, 0.5 * MeV);
401
402 mod = new G4DNABornIonisationModel();
403 em_config->SetExtraEmModel("proton", "proton_G4DNAIonisation", mod, "Target", 0.5 * MeV,
404 10 * MeV);
405
406 mod = new G4DNAMillerGreenExcitationModel();
407 em_config->SetExtraEmModel("proton", "proton_G4DNAExcitation", mod, "Target", 10 * eV, 0.5 * MeV);
408
409 mod = new G4DNABornExcitationModel();
410 em_config->SetExtraEmModel("proton", "proton_G4DNAExcitation", mod, "Target", 0.5 * MeV,
411 10 * MeV);
412
413 // *** alpha
414
415 // ---> STANDARD EM processes inactivated below standEnergyLimit
416
417 // STANDARD msc is still active
418 // Inactivate following STANDARD processes
419
420 mod = new G4BraggIonGasModel();
421 mod->SetActivationLowEnergyLimit(standEnergyLimit);
422 em_config->SetExtraEmModel("alpha", "ionIoni", mod, "Target", 0.0, 2 * MeV / massFactor,
423 new G4IonFluctuations());
424
425 mod = new G4BetheBlochIonGasModel();
426 mod->SetActivationLowEnergyLimit(standEnergyLimit);
427 em_config->SetExtraEmModel("alpha", "ionIoni", mod, "Target", 2 * MeV / massFactor, 100 * TeV,
428 new G4UniversalFluctuation());
429
430 // ---> DNA processes activated
431
432 mod = new G4DNARuddIonisationModel();
433 em_config->SetExtraEmModel("alpha", "alpha_G4DNAIonisation", mod, "Target", 0.0, 10 * MeV);
434
435 mod = new G4DNAMillerGreenExcitationModel();
436 em_config->SetExtraEmModel("alpha", "alpha_G4DNAExcitation", mod, "Target", 1 * keV, 10 * MeV);
437 }
438
439 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
440
441 void PhysicsList::ConstructGeneral() {}
442
443 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
444
445 void PhysicsList::SetCuts()
446 {
447 if (verboseLevel > 0) {
448 G4cout << "PhysicsList::SetCuts:";
449 G4cout << "CutLength : " << G4BestUnit(defaultCutValue, "Length") << G4endl;
450 }
451
452 // set cut values for gamma at first and for e- second and next for e+,
453 // because some processes for e+/e- need cut values for gamma
454 SetCutValue(fCutForGamma, "gamma");
455 SetCutValue(fCutForElectron, "e-");
456 SetCutValue(fCutForPositron, "e+");
457 SetCutValue(fCutForProton, "proton");
458
459 if (verboseLevel > 0) {
460 DumpCutValuesTable();
461 }
462 }
463
464 void PhysicsList::SetNumberOfSplit(G4int nb)
465 {
466 fNumberOfSplit = nb;
467 }
468
469 G4int PhysicsList::GetNumberOfSplit()
470 {
471 return fNumberOfSplit;
472 }
473