| 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 //---------------------------------------------------------------------------
28 //
29 // ClassName: G4EmModelActivator
30 //
31 // Author: V.Ivanchenko 01.06.2015
32 //
33 // Organisation: G4AI
34 // Contract: ESA contract 4000107387/12/NL/AT
35 // Customer: ESA/ESTEC
36 //
37 // Modified:
38 //
39 //----------------------------------------------------------------------------
40 //
41
42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
44
45 #include "G4EmModelActivator.hh"
46 #include "G4EmParameters.hh"
47 #include "G4ParticleDefinition.hh"
48 #include "G4ParticleTable.hh"
49 #include "G4RegionStore.hh"
50 #include "G4Region.hh"
51 #include "G4VEnergyLossProcess.hh"
52 #include "G4VMscModel.hh"
53 #include "G4LossTableManager.hh"
54 #include "G4EmConfigurator.hh"
55 #include "G4PAIModel.hh"
56 #include "G4PAIPhotModel.hh"
57 #include "G4Gamma.hh"
58 #include "G4Electron.hh"
59 #include "G4Positron.hh"
60 #include "G4MuonPlus.hh"
61 #include "G4MuonMinus.hh"
62 #include "G4Proton.hh"
63 #include "G4GenericIon.hh"
64 #include "G4Alpha.hh"
65 #include "G4ProcessManager.hh"
66 #include "G4DummyModel.hh"
67 #include "G4EmProcessSubType.hh"
68 #include "G4PhysicsListHelper.hh"
69 #include "G4EmParticleList.hh"
70 #include "G4EmUtility.hh"
71
72 #include "G4MicroElecElastic.hh"
73 #include "G4MicroElecElasticModel.hh"
74
75 #include "G4MicroElecInelastic.hh"
76 #include "G4MicroElecInelasticModel.hh"
77
78 #include "G4BraggModel.hh"
79 #include "G4BraggIonModel.hh"
80 #include "G4BetheBlochModel.hh"
81 #include "G4UrbanMscModel.hh"
82 #include "G4MollerBhabhaModel.hh"
83 #include "G4IonFluctuations.hh"
84 #include "G4UniversalFluctuation.hh"
85 #include "G4LowECapture.hh"
86 #include "G4hMultipleScattering.hh"
87 #include "G4eCoulombScatteringModel.hh"
88 #include "G4IonCoulombScatteringModel.hh"
89 #include "G4ionIonisation.hh"
90 #include "G4KleinNishinaModel.hh"
91
92 #include "G4GammaGeneralProcess.hh"
93 #include "G4CoulombScattering.hh"
94 #include "G4eCoulombScatteringModel.hh"
95 #include "G4WentzelVIModel.hh"
96 #include "G4UniversalFluctuation.hh"
97 #include "G4RayleighScattering.hh"
98 #include "G4UrbanMscModel.hh"
99 #include "G4GoudsmitSaundersonMscModel.hh"
100 #include "G4LowEPComptonModel.hh"
101 #include "G4BetheHeitler5DModel.hh"
102 #include "G4LindhardSorensenIonModel.hh"
103
104 #include "G4LivermorePhotoElectricModel.hh"
105 #include "G4LivermoreComptonModel.hh"
106 #include "G4LivermoreGammaConversionModel.hh"
107 #include "G4LivermoreRayleighModel.hh"
108 #include "G4LivermoreIonisationModel.hh"
109 #include "G4LivermoreBremsstrahlungModel.hh"
110
111 #include "G4PenelopePhotoElectricModel.hh"
112 #include "G4PenelopeComptonModel.hh"
113 #include "G4PenelopeGammaConversionModel.hh"
114 #include "G4PenelopeRayleighModel.hh"
115 #include "G4PenelopeIonisationModel.hh"
116 #include "G4PenelopeBremsstrahlungModel.hh"
117 #include "G4PenelopeAnnihilationModel.hh"
118
119 #include "G4SystemOfUnits.hh"
120 #include "G4Threading.hh"
121 #include <vector>
122
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
124
125 G4EmModelActivator::G4EmModelActivator(const G4String& emphys)
126 : baseName(emphys)
127 {
128 theParameters = G4EmParameters::Instance();
129
130 const std::vector<G4String>& regnamesPAI = theParameters->RegionsPAI();
131 if(regnamesPAI.size() > 0)
132 {
133 ActivatePAI();
134 }
135 const std::vector<G4String>& regnamesME = theParameters->RegionsMicroElec();
136 if(regnamesME.size() > 0)
137 {
138 ActivateMicroElec();
139 }
140 const std::vector<G4String>& regnamesMSC = theParameters->RegionsPhysics();
141 if(regnamesMSC.size() > 0)
142 {
143 ActivateEmOptions();
144 }
145 }
146
147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
148
149 void G4EmModelActivator::ActivateEmOptions()
150 {
151 const std::vector<G4String>& regnamesPhys = theParameters->RegionsPhysics();
152 std::size_t nreg = regnamesPhys.size();
153 if(0 == nreg) { return; }
154 G4int verbose = theParameters->Verbose() - 1;
155 if(verbose > 0) {
156 G4cout << "### G4EmModelActivator::ActivateEmOptions for " << nreg << " regions"
157 << G4endl;
158 }
159 const std::vector<G4String>& typesPhys = theParameters->TypesPhysics();
160
161 // start configuration of models
162 const G4ParticleDefinition* elec = G4Electron::Electron();
163 const G4ParticleDefinition* posi = G4Positron::Positron();
164 const G4ParticleDefinition* phot = G4Gamma::Gamma();
165 const G4ParticleDefinition* prot = G4Proton::Proton();
166 G4LossTableManager* man = G4LossTableManager::Instance();
167 G4EmConfigurator* em_config = man->EmConfigurator();
168 G4VAtomDeexcitation* adeexc = man->AtomDeexcitation();
169 G4ParticleTable* table = G4ParticleTable::GetParticleTable();
170 G4VEmModel* mod = nullptr;
171 G4VEmProcess* proc = nullptr;
172
173 // high energy limit for low-energy e+- model of msc
174 G4double mscEnergyLimit = theParameters->MscEnergyLimit();
175
176 // high energy limit for Livermore and Penelope models
177 G4double highEnergyLimit = 1*GeV;
178
179 // general high energy limit
180 G4double highEnergy = theParameters->MaxKinEnergy();
181
182 for(std::size_t i=0; i<nreg; ++i) {
183 const G4String reg = regnamesPhys[i];
184 auto region = G4EmUtility::FindRegion(reg);
185 if(verbose > 0) {
186 G4cout << i << ". region <" << reg << ">; physics type <"
187 << typesPhys[i] << "> region ptr: " << region << G4endl;
188 }
189
190 if(baseName == typesPhys[i]) { continue; }
191
192 if("G4EmStandard" == typesPhys[i]) {
193 G4UrbanMscModel* msc = new G4UrbanMscModel();
194 AddStandardScattering(elec, em_config, msc, reg, mscEnergyLimit, highEnergy, typesPhys[i]);
195
196 msc = new G4UrbanMscModel();
197 AddStandardScattering(posi, em_config, msc, reg, mscEnergyLimit, highEnergy, typesPhys[i]);
198
199 } else if("G4EmStandard_opt1" == typesPhys[i] || "G4EmStandard_opt2" == typesPhys[i]) {
200 G4UrbanMscModel* msc = new G4UrbanMscModel();
201 AddStandardScattering(elec, em_config, msc, reg, mscEnergyLimit, highEnergy, typesPhys[i]);
202
203 msc = new G4UrbanMscModel();
204 AddStandardScattering(posi, em_config, msc, reg, mscEnergyLimit, highEnergy, typesPhys[i]);
205
206 } else if("G4EmStandard_opt3" == typesPhys[i]) {
207
208 G4DummyModel* dummy = new G4DummyModel();
209 G4UrbanMscModel* msc = new G4UrbanMscModel();
210 SetMscParameters(elec, msc, typesPhys[i]);
211 em_config->SetExtraEmModel("e-", "msc", msc, reg);
212 proc = FindOrAddProcess(elec, "CoulombScat");
213 proc->AddEmModel(-1, dummy, region);
214
215 msc = new G4UrbanMscModel();
216 SetMscParameters(posi, msc, typesPhys[i]);
217 em_config->SetExtraEmModel("e+", "msc", msc, reg);
218 proc = FindOrAddProcess(posi, "CoulombScat");
219 proc->AddEmModel(-1, dummy, region);
220
221 msc = new G4UrbanMscModel();
222 SetMscParameters(prot, msc, typesPhys[i]);
223 em_config->SetExtraEmModel("proton", "msc", msc, reg);
224 proc = FindOrAddProcess(prot, "CoulombScat");
225 proc->AddEmModel(-1, dummy, region);
226
227 theParameters->SetNumberOfBinsPerDecade(20);
228 if(G4Threading::IsMasterThread()) {
229 theParameters->SetDeexActiveRegion(reg, true, false, false);
230 }
231 theParameters->DefineRegParamForDeex(adeexc);
232 proc = FindOrAddProcess(phot, "Rayl");
233 proc->AddEmModel(-1, new G4LivermoreRayleighModel(), region);
234 proc = FindOrAddProcess(phot, "compt");
235 proc->AddEmModel(-1, new G4KleinNishinaModel(), region);
236
237 } else if("G4EmStandard_opt4" == typesPhys[i]) {
238 G4VMscModel* msc = new G4GoudsmitSaundersonMscModel();
239 AddStandardScattering(elec, em_config, msc, reg, mscEnergyLimit, highEnergy, typesPhys[i]);
240
241 msc = new G4GoudsmitSaundersonMscModel();
242 AddStandardScattering(posi, em_config, msc, reg, mscEnergyLimit, highEnergy, typesPhys[i]);
243
244 theParameters->SetNumberOfBinsPerDecade(20);
245 theParameters->SetUseMottCorrection(true);
246 if(G4Threading::IsMasterThread()) {
247 theParameters->SetDeexActiveRegion(reg, true, false, false);
248 }
249 theParameters->DefineRegParamForDeex(adeexc);
250
251 proc = FindOrAddProcess(phot, "Rayl");
252 proc->AddEmModel(-1, new G4LivermoreRayleighModel(), region);
253 proc = FindOrAddProcess(phot, "compt");
254 proc->AddEmModel(-1, new G4KleinNishinaModel(), region);
255 mod = new G4LowEPComptonModel();
256 mod->SetHighEnergyLimit(20*MeV);
257 proc->AddEmModel(-2, mod, region);
258 proc = FindOrAddProcess(phot, "conv");
259 proc->AddEmModel(-1, new G4BetheHeitler5DModel(), region);
260
261 } else if("G4EmStandardGS" == typesPhys[i]) {
262 G4GoudsmitSaundersonMscModel* msc = new G4GoudsmitSaundersonMscModel();
263 AddStandardScattering(elec, em_config, msc, reg, mscEnergyLimit, highEnergy, typesPhys[i]);
264
265 msc = new G4GoudsmitSaundersonMscModel();
266 AddStandardScattering(posi, em_config, msc, reg, mscEnergyLimit, highEnergy, typesPhys[i]);
267
268 } else if("G4EmStandardWVI" == typesPhys[i]) {
269 G4WentzelVIModel* msc = new G4WentzelVIModel();
270 AddStandardScattering(elec, em_config, msc, reg, mscEnergyLimit, highEnergy, typesPhys[i]);
271
272 msc = new G4WentzelVIModel();
273 AddStandardScattering(posi, em_config, msc, reg, mscEnergyLimit, highEnergy, typesPhys[i]);
274
275 theParameters->SetMscThetaLimit(0.15);
276
277 if(G4Threading::IsMasterThread()) {
278 theParameters->SetDeexActiveRegion(regnamesPhys[i], true, false, false);
279 }
280 theParameters->DefineRegParamForDeex(adeexc);
281
282 } else if("G4EmStandardSS" == typesPhys[i]) {
283 G4EmParticleList emList;
284 for(const auto& particleName : emList.EmChargedPartNames()) {
285 G4ParticleDefinition* particle = table->FindParticle(particleName);
286 if(nullptr != particle && 0.0 != particle->GetPDGCharge()) {
287 proc = FindOrAddProcess(particle, "CoulombScat");
288 if(nullptr != proc) {
289 proc->AddEmModel(-1, new G4DummyModel(), region);
290 }
291 auto pm = particle->GetProcessManager();
292 proc = new G4CoulombScattering("SingleCoulombScat", false);
293 pm->AddDiscreteProcess(proc);
294 proc->SetEmModel(new G4DummyModel());
295 G4VEmModel* scmod = nullptr;
296 if(particle->GetPDGMass() > CLHEP::GeV ||
297 particle->GetParticleType() == "nucleus") {
298 scmod = new G4IonCoulombScatteringModel();
299 } else {
300 scmod = new G4eCoulombScatteringModel(false);
301 }
302 scmod->SetPolarAngleLimit(0.0);
303 scmod->SetLocked(true);
304 proc->AddEmModel(-1, scmod, region);
305
306 // multiple scattering should be disabled
307 if(particleName == "mu+" || particleName == "mu-") {
308 em_config->SetExtraEmModel(particleName, "muMsc",
309 new G4DummyModel(), reg);
310 } else {
311 em_config->SetExtraEmModel(particleName, "msc",
312 new G4DummyModel(), reg);
313 }
314 }
315 }
316 if(G4Threading::IsMasterThread()) {
317 theParameters->SetDeexActiveRegion(regnamesPhys[i], true, true, true);
318 }
319 theParameters->DefineRegParamForDeex(adeexc);
320
321 } else if("G4EmLivermore" == typesPhys[i]) {
322
323 G4VMscModel* msc = new G4GoudsmitSaundersonMscModel();
324 AddStandardScattering(elec, em_config, msc, reg, mscEnergyLimit, highEnergy, typesPhys[i]);
325
326 msc = new G4GoudsmitSaundersonMscModel();
327 AddStandardScattering(posi, em_config, msc, reg, mscEnergyLimit, highEnergy, typesPhys[i]);
328
329 mod = new G4LivermorePhotoElectricModel();
330 em_config->SetExtraEmModel("gamma", "phot", mod, reg);
331 mod = new G4LivermoreComptonModel();
332 em_config->SetExtraEmModel("gamma", "compt", mod, reg);
333 mod = new G4LivermoreGammaConversionModel();
334 em_config->SetExtraEmModel("gamma", "conv", mod, reg);
335
336 FindOrAddProcess(phot, "Rayl");
337 mod = new G4LivermoreRayleighModel();
338 em_config->SetExtraEmModel("gamma", "Rayl", mod, reg);
339
340 mod = new G4LivermoreIonisationModel();
341 G4UniversalFluctuation* uf = new G4UniversalFluctuation();
342 em_config->SetExtraEmModel("e-", "eIoni", mod, reg, 0.0, 0.1*MeV, uf);
343 mod = new G4LivermoreBremsstrahlungModel();
344 em_config->SetExtraEmModel("e-", "eBrem", mod, reg, 0.0, highEnergyLimit);
345
346 theParameters->SetNumberOfBinsPerDecade(20);
347 theParameters->SetUseMottCorrection(true);
348 if(G4Threading::IsMasterThread()) {
349 theParameters->SetDeexActiveRegion(regnamesPhys[i], true, false, false);
350 }
351 theParameters->DefineRegParamForDeex(adeexc);
352
353 } else if("G4EmPenelope" == typesPhys[i]) {
354
355 G4VMscModel* msc = new G4GoudsmitSaundersonMscModel();
356 AddStandardScattering(elec, em_config, msc, reg, mscEnergyLimit, highEnergy, typesPhys[i]);
357
358 msc = new G4GoudsmitSaundersonMscModel();
359 AddStandardScattering(posi, em_config, msc, reg, mscEnergyLimit, highEnergy, typesPhys[i]);
360
361 mod = new G4PenelopePhotoElectricModel();
362 em_config->SetExtraEmModel("gamma", "phot", mod, reg);
363 mod = new G4PenelopeComptonModel();
364 em_config->SetExtraEmModel("gamma", "compt", mod, reg);
365 mod = new G4PenelopeGammaConversionModel();
366 em_config->SetExtraEmModel("gamma", "conv", mod, reg);
367
368 FindOrAddProcess(phot, "Rayl");
369 mod = new G4PenelopeRayleighModel();
370 em_config->SetExtraEmModel("gamma", "Rayl", mod, reg);
371
372 mod = new G4PenelopeIonisationModel();
373 G4UniversalFluctuation* uf = new G4UniversalFluctuation();
374 em_config->SetExtraEmModel("e-", "eIoni", mod, reg, 0.0, highEnergyLimit, uf);
375 mod = new G4PenelopeBremsstrahlungModel();
376 em_config->SetExtraEmModel("e-", "eBrem", mod, reg, 0.0, highEnergyLimit);
377
378 mod = new G4PenelopeIonisationModel();
379 uf = new G4UniversalFluctuation();
380 em_config->SetExtraEmModel("e+", "eIoni", mod, reg, 0.0, highEnergyLimit, uf);
381 mod = new G4PenelopeBremsstrahlungModel();
382 em_config->SetExtraEmModel("e+", "eBrem", mod, reg, 0.0, highEnergyLimit);
383 mod = new G4PenelopeAnnihilationModel();
384 em_config->SetExtraEmModel("e+", "annihil", mod, reg, 0.0, highEnergyLimit);
385
386 theParameters->SetNumberOfBinsPerDecade(20);
387 theParameters->SetUseMottCorrection(true);
388 if(G4Threading::IsMasterThread()) {
389 theParameters->SetDeexActiveRegion(regnamesPhys[i], true, false, false);
390 }
391 theParameters->DefineRegParamForDeex(adeexc);
392
393 } else {
394 if(verbose > 0 && G4Threading::IsMasterThread()) {
395 G4cout << "### G4EmModelActivator::ActivateEmOptions WARNING: \n"
396 << " EM Physics configuration name <" << typesPhys[i]
397 << "> is not known - ignored" << G4endl;
398 }
399 }
400 }
401 }
402
403 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
404
405 void G4EmModelActivator::ActivatePAI()
406 {
407 const std::vector<G4String> regnamesPAI = theParameters->RegionsPAI();
408 std::size_t nreg = regnamesPAI.size();
409 if(0 == nreg) { return; }
410 G4int verbose = theParameters->Verbose() - 1;
411 if(verbose > 0) {
412 G4cout << "### G4EmModelActivator::ActivatePAI for " << nreg << " regions"
413 << G4endl;
414 }
415 const std::vector<G4String> particlesPAI = theParameters->ParticlesPAI();
416 const std::vector<G4String> typesPAI = theParameters->TypesPAI();
417
418 const std::vector<G4VEnergyLossProcess*>& v = G4LossTableManager::Instance()
419 ->GetEnergyLossProcessVector();
420
421 G4RegionStore* regionStore = G4RegionStore::GetInstance();
422
423 const G4ParticleDefinition* elec = G4Electron::Electron();
424 const G4ParticleDefinition* posi = G4Positron::Positron();
425 const G4ParticleDefinition* mupl = G4MuonPlus::MuonPlus();
426 const G4ParticleDefinition* mumi = G4MuonMinus::MuonMinus();
427 const G4ParticleDefinition* gion = G4GenericIon::GenericIon();
428
429 for(std::size_t i = 0; i < nreg; ++i) {
430 const G4ParticleDefinition* p = nullptr;
431 if(particlesPAI[i] != "all") {
432 p = G4ParticleTable::GetParticleTable()->FindParticle(particlesPAI[i]);
433 if(!p) {
434 G4cout << "### WARNING: ActivatePAI::FindParticle fails to find "
435 << particlesPAI[i] << G4endl;
436 continue;
437 }
438 }
439 const G4Region* r = regionStore->GetRegion(regnamesPAI[i], false);
440 if(!r) {
441 G4cout << "### WARNING: ActivatePAI::GetRegion fails to find "
442 << regnamesPAI[i] << G4endl;
443 continue;
444 }
445
446 G4String name = "hIoni";
447 if(p == elec || p == posi)
448 { name = "eIoni"; }
449 else if (p == mupl || p == mumi)
450 { name = "muIoni"; }
451 else if (p == gion)
452 { name = "ionIoni"; }
453
454 for(auto proc : v) {
455
456 if(!proc->IsIonisationProcess()) { continue; }
457
458 G4String namep = proc->GetProcessName();
459 if(p) {
460 if(name != namep) { continue; }
461 } else {
462 if(namep != "hIoni" && namep != "muIoni" &&
463 namep != "eIoni" && namep != "ionIoni")
464 { continue; }
465 }
466 G4double emin = 50*CLHEP::keV;
467 if(namep == "eIoni") emin = 110*CLHEP::eV;
468 else if(namep == "muIoni") emin = 5*CLHEP::keV;
469
470 G4VEmModel* em = nullptr;
471 G4VEmFluctuationModel* fm = nullptr;
472 if(typesPAI[i] == "PAIphoton" || typesPAI[i] == "pai_photon") {
473 G4PAIPhotModel* mod = new G4PAIPhotModel(p,"PAIPhotModel");
474 em = mod;
475 fm = mod;
476 } else {
477 G4PAIModel* mod = new G4PAIModel(p,"PAIModel");
478 em = mod;
479 fm = mod;
480 }
481 // first added the default model for world
482 // second added low energy model below PAI threshold in the region
483 // finally added PAI for the region
484 G4VEmModel* em0 = nullptr;
485 G4VEmFluctuationModel* fm0 = nullptr;
486 if(namep == "eIoni") {
487 fm0 = new G4UniversalFluctuation();
488 proc->SetEmModel(new G4MollerBhabhaModel());
489 proc->SetFluctModel(fm0);
490 em0 = new G4MollerBhabhaModel();
491 } else if(namep == "ionIoni") {
492 fm0 = new G4IonFluctuations();
493 proc->SetEmModel(new G4LindhardSorensenIonModel());
494 proc->SetFluctModel(fm0);
495 em0 = new G4LindhardSorensenIonModel();
496 } else {
497 fm0 = new G4UniversalFluctuation();
498 proc->SetEmModel(new G4BraggModel());
499 proc->SetEmModel(new G4BetheBlochModel());
500 proc->SetFluctModel(fm0);
501 em0 = new G4BraggModel();
502 }
503 em0->SetHighEnergyLimit(emin);
504 proc->AddEmModel(-1, em0, fm0, r);
505 em->SetLowEnergyLimit(emin);
506 proc->AddEmModel(-1, em, fm, r);
507 if(verbose > 0) {
508 G4cout << "### G4EmModelActivator: add <" << typesPAI[i]
509 << "> model for " << particlesPAI[i]
510 << " in the " << regnamesPAI[i]
511 << " Emin(keV)= " << emin/CLHEP::keV << G4endl;
512 }
513 }
514 }
515 }
516
517 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
518
519 void G4EmModelActivator::ActivateMicroElec()
520 {
521 const std::vector<G4String> regnamesME = theParameters->RegionsMicroElec();
522 std::size_t nreg = regnamesME.size();
523 if(0 == nreg)
524 {
525 return;
526 }
527 G4int verbose = theParameters->Verbose() - 1;
528 if(verbose > 0)
529 {
530 G4cout << "### G4EmModelActivator::ActivateMicroElec for " << nreg
531 << " regions" << G4endl;
532 }
533 G4LossTableManager* man = G4LossTableManager::Instance();
534
535 const G4ParticleDefinition* elec = G4Electron::Electron();
536 const G4ParticleDefinition* prot = G4Proton::Proton();
537 const G4ParticleDefinition* gion = G4GenericIon::GenericIon();
538 G4ProcessManager* eman = elec->GetProcessManager();
539 G4ProcessManager* pman = prot->GetProcessManager();
540 G4ProcessManager* iman = gion->GetProcessManager();
541
542 G4bool emsc = HasMsc(eman);
543
544 // MicroElec elastic is not active in the world
545 G4MicroElecElastic* eElasProc =
546 new G4MicroElecElastic("e-G4MicroElecElastic");
547 eman->AddDiscreteProcess(eElasProc);
548
549 G4MicroElecInelastic* eInelProc =
550 new G4MicroElecInelastic("e-G4MicroElecInelastic");
551 eman->AddDiscreteProcess(eInelProc);
552
553 G4MicroElecInelastic* pInelProc =
554 new G4MicroElecInelastic("p_G4MicroElecInelastic");
555 pman->AddDiscreteProcess(pInelProc);
556
557 G4MicroElecInelastic* iInelProc =
558 new G4MicroElecInelastic("ion_G4MicroElecInelastic");
559 iman->AddDiscreteProcess(iInelProc);
560
561 // start configuration of models
562 G4EmConfigurator* em_config = man->EmConfigurator();
563 G4VEmModel* mod;
564
565 // limits for MicroElec applicability
566 G4double elowest = 16.7 * eV;
567 G4double elimel = 100 * MeV;
568 G4double elimin = 9 * MeV;
569 G4double pmin = 50 * keV;
570 G4double pmax = 99.9 * MeV;
571
572 G4LowECapture* ecap = new G4LowECapture(elowest);
573 eman->AddDiscreteProcess(ecap);
574
575 for(std::size_t i = 0; i < nreg; ++i)
576 {
577
578 G4String reg = regnamesME[i];
579 G4cout << "### MicroElec models are activated for G4Region " << reg
580 << G4endl
581 << " Energy limits for e- elastic: " << elowest/eV << " eV - "
582 << elimel/MeV << " MeV"
583 << G4endl
584 << " Energy limits for e- inelastic: " << elowest/eV << " eV - "
585 << elimin/MeV << " MeV"
586 << G4endl
587 << " Energy limits for hadrons/ions: " << pmin/MeV << " MeV - "
588 << pmax/MeV << " MeV"
589 << G4endl;
590
591 // e-
592 if(emsc)
593 {
594 G4UrbanMscModel* msc = new G4UrbanMscModel();
595 msc->SetActivationLowEnergyLimit(elimel);
596 em_config->SetExtraEmModel("e-", "msc", msc, reg);
597 }
598 else
599 {
600 mod = new G4DummyModel();
601 em_config->SetExtraEmModel("e-", "CoulombScat", mod, reg, 0.0, elimel);
602 }
603
604 mod = new G4MicroElecElasticModel();
605 em_config->SetExtraEmModel("e-",
606 "e-G4MicroElecElastic",
607 mod,
608 reg,
609 elowest,
610 elimin);
611
612 mod = new G4MollerBhabhaModel();
613 mod->SetActivationLowEnergyLimit(elimin);
614 em_config->SetExtraEmModel("e-",
615 "eIoni",
616 mod,
617 reg,
618 0.0,
619 10 * TeV,
620 new G4UniversalFluctuation());
621
622 mod = new G4MicroElecInelasticModel();
623 em_config->SetExtraEmModel("e-",
624 "e-G4MicroElecInelastic",
625 mod,
626 reg,
627 elowest,
628 elimin);
629
630 // proton
631 mod = new G4BraggModel();
632 mod->SetActivationHighEnergyLimit(pmin);
633 em_config->SetExtraEmModel("proton",
634 "hIoni",
635 mod,
636 reg,
637 0.0,
638 2 * MeV,
639 new G4UniversalFluctuation());
640
641 mod = new G4BetheBlochModel();
642 mod->SetActivationLowEnergyLimit(pmax);
643 em_config->SetExtraEmModel("proton",
644 "hIoni",
645 mod,
646 reg,
647 2 * MeV,
648 10 * TeV,
649 new G4UniversalFluctuation());
650
651 mod = new G4MicroElecInelasticModel();
652 em_config->SetExtraEmModel("proton",
653 "p_G4MicroElecInelastic",
654 mod,
655 reg,
656 pmin,
657 pmax);
658
659 // ions
660 mod = new G4BraggIonModel();
661 mod->SetActivationHighEnergyLimit(pmin);
662 em_config->SetExtraEmModel("GenericIon",
663 "ionIoni",
664 mod,
665 reg,
666 0.0,
667 2 * MeV,
668 new G4IonFluctuations());
669
670 mod = new G4BetheBlochModel();
671 mod->SetActivationLowEnergyLimit(pmax);
672 em_config->SetExtraEmModel("GenericIon",
673 "ionIoni",
674 mod,
675 reg,
676 2 * MeV,
677 10 * TeV,
678 new G4IonFluctuations());
679
680 mod = new G4MicroElecInelasticModel();
681 em_config->SetExtraEmModel("GenericIon",
682 "ion_G4MicroElecInelastic",
683 mod,
684 reg,
685 pmin,
686 pmax);
687 }
688 }
689
690 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
691
692 G4bool G4EmModelActivator::HasMsc(G4ProcessManager* pm) const
693 {
694 G4bool res = false;
695 G4ProcessVector* pv = pm->GetProcessList();
696 G4int nproc = pm->GetProcessListLength();
697 for(G4int i = 0; i < nproc; ++i)
698 {
699 if(((*pv)[i])->GetProcessSubType() == fMultipleScattering)
700 {
701 res = true;
702 break;
703 }
704 }
705 return res;
706 }
707
708 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
709
710 void G4EmModelActivator::AddStandardScattering(const G4ParticleDefinition* part,
711 G4EmConfigurator* em_config,
712 G4VMscModel* mscmod,
713 const G4String& reg,
714 G4double e1, G4double e2,
715 const G4String& type)
716 {
717 G4String pname = part->GetParticleName();
718
719 // low-energy msc model
720 SetMscParameters(part, mscmod, type);
721 em_config->SetExtraEmModel(pname, "msc", mscmod, reg, 0.0, e1);
722
723 // high energy msc model
724 G4WentzelVIModel* msc = new G4WentzelVIModel();
725 SetMscParameters(part, msc, type);
726 em_config->SetExtraEmModel(pname, "msc", msc, reg, e1, e2);
727
728 // high energy single scattering model
729 FindOrAddProcess(part, "CoulombScat");
730 G4eCoulombScatteringModel* mod = new G4eCoulombScatteringModel();
731 mod->SetActivationLowEnergyLimit(e1);
732 mod->SetLocked(true);
733 em_config->SetExtraEmModel(pname, "CoulombScat", mod, reg, 0.0, e2);
734 }
735
736 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
737
738 void G4EmModelActivator::SetMscParameters(const G4ParticleDefinition* part,
739 G4VMscModel* msc, const G4String& phys)
740 {
741 if(part == G4Electron::Electron() || part == G4Positron::Positron()) {
742 if(phys == "G4EmStandard_opt1" || phys == "G4EmStandard_opt2") {
743 msc->SetRangeFactor(0.2);
744 msc->SetStepLimitType(fMinimal);
745 } else if(phys == "G4EmStandard_opt3") {
746 msc->SetStepLimitType(fUseDistanceToBoundary);
747 } else if(phys == "G4EmStandard_opt4" || phys == "G4EmLivermore" || phys == "G4EmPenelope") {
748 msc->SetRangeFactor(0.08);
749 msc->SetStepLimitType(fUseSafetyPlus);
750 msc->SetSkin(3);
751 } else if(phys == "G4EmStandardGS") {
752 msc->SetRangeFactor(0.06);
753 }
754 } else {
755 if(phys != "G4EmStandard" && phys != "G4EmStandard_opt1" && phys != "G4EmStandard_opt2") {
756 msc->SetLateralDisplasmentFlag(true);
757 }
758 }
759 msc->SetLocked(true);
760 }
761
762 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
763
764 G4VEmProcess*
765 G4EmModelActivator::FindOrAddProcess(const G4ParticleDefinition* part,
766 const G4String& name)
767 {
768 G4VEmProcess* proc = nullptr;
769 auto pm = part->GetProcessManager();
770 auto emproc = pm->GetProcessList();
771 G4int n = (G4int)emproc->size();
772 for(auto i=0; i<n; ++i) {
773 auto ptr = (*emproc)[i];
774 if(part->GetPDGEncoding() == 22 &&
775 ptr->GetProcessSubType() == fGammaGeneralProcess) {
776 proc = (static_cast<G4GammaGeneralProcess*>(ptr))->GetEmProcess(name);
777 } else if(ptr->GetProcessName() == name) {
778 proc = dynamic_cast<G4VEmProcess*>(ptr);
779 }
780 if(nullptr != proc) { return proc; }
781 }
782 if(name == "Rayl") {
783 G4RayleighScattering* rs = new G4RayleighScattering();
784 rs->SetEmModel(new G4DummyModel());
785 pm->AddDiscreteProcess(rs);
786 return rs;
787 }
788 return proc;
789 }
790
791 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
792