| 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 // Author V.Ivanchenko
27 //
28
29 #include "G4EmDNAPhysicsActivator.hh"
30
31 #include "G4EmParameters.hh"
32 #include "G4ParticleDefinition.hh"
33 #include "G4ParticleTable.hh"
34 #include "G4Region.hh"
35 #include "G4VEnergyLossProcess.hh"
36
37 #include "G4Gamma.hh"
38 #include "G4Electron.hh"
39 #include "G4Positron.hh"
40 #include "G4Proton.hh"
41 #include "G4GenericIon.hh"
42 #include "G4Alpha.hh"
43
44 #include "G4ProcessManager.hh"
45 #include "G4DummyModel.hh"
46 #include "G4EmProcessSubType.hh"
47 #include "G4PhysicsListHelper.hh"
48
49 #include "G4BraggModel.hh"
50 #include "G4BraggIonModel.hh"
51 #include "G4BetheBlochModel.hh"
52 #include "G4UrbanMscModel.hh"
53 #include "G4GoudsmitSaundersonMscModel.hh"
54 #include "G4MollerBhabhaModel.hh"
55 #include "G4SeltzerBergerModel.hh"
56 #include "G4IonFluctuations.hh"
57 #include "G4UniversalFluctuation.hh"
58 #include "G4IonFluctuations.hh"
59 #include "G4LowECapture.hh"
60 #include "G4eMultipleScattering.hh"
61 #include "G4hMultipleScattering.hh"
62 #include "G4eCoulombScatteringModel.hh"
63 #include "G4IonCoulombScatteringModel.hh"
64 #include "G4eIonisation.hh"
65 #include "G4eBremsstrahlung.hh"
66 #include "G4hIonisation.hh"
67 #include "G4ionIonisation.hh"
68 #include "G4NuclearStopping.hh"
69 #include "G4ICRU49NuclearStoppingModel.hh"
70 #include "G4Generator2BS.hh"
71
72 #include "G4Threading.hh"
73 #include "G4EmDNABuilder.hh"
74 #include "G4EmUtility.hh"
75 #include "G4PhysListUtil.hh"
76 #include "G4SystemOfUnits.hh"
77 #include <vector>
78
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
80
81 G4EmDNAPhysicsActivator::G4EmDNAPhysicsActivator(G4int ver)
82 : G4VPhysicsConstructor("G4EmDNAPhysicsActivator"), verbose(ver)
83 {
84 theParameters = G4EmParameters::Instance();
85 theParameters->ActivateDNA();
86 theParameters->SetFluo(true);
87 }
88
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90
91 G4bool G4EmDNAPhysicsActivator::IsVerbose() const
92 {
93 return (0 < verbose && G4Threading::IsMasterThread());
94 }
95
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
97
98 void G4EmDNAPhysicsActivator::ConstructParticle()
99 {
100 G4EmDNABuilder::ConstructDNAParticles();
101 }
102
103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
104
105 void G4EmDNAPhysicsActivator::ConstructProcess()
106 {
107 const std::vector<G4String>& regnamesDNA = theParameters->RegionsDNA();
108 std::size_t nreg = regnamesDNA.size();
109 if(0 == nreg)
110 {
111 return;
112 }
113
114 const std::vector<G4String>& typesDNA = theParameters->TypesDNA();
115 G4bool fast = theParameters->DNAFast();
116 G4bool st = theParameters->DNAStationary();
117
118 const G4double emaxDNA = 1.*CLHEP::MeV;
119 const G4double emaxIonDNA = 300.*CLHEP::MeV;
120 const G4double emaxLightIonDNA = 400.*CLHEP::MeV;
121 const G4double eminBorn = 500.*CLHEP::keV;
122 const G4double emax = theParameters->MaxKinEnergy();
123
124 if(IsVerbose()) {
125 G4cout << "### G4EmDNAPhysicsActivator::ConstructProcess for " << nreg
126 << " regions; DNA physics type " << G4endl;
127 }
128
129 // list of particles
130 G4ParticleDefinition* prot = G4Proton::Proton();
131 G4ParticleDefinition* gion = G4GenericIon::GenericIon();
132
133 G4DNAGenericIonsManager * genericIonsManager =
134 G4DNAGenericIonsManager::Instance();
135 G4ParticleDefinition* alpha2 = G4Alpha::Alpha();
136 G4ParticleDefinition* alpha1 = genericIonsManager->GetIon("alpha+");
137 G4ParticleDefinition* alpha0 = genericIonsManager->GetIon("helium");
138 G4ParticleDefinition* h0 = genericIonsManager->GetIon("hydrogen");
139
140 // loop over regions
141 for(std::size_t i = 0; i < nreg; ++i)
142 {
143 if(IsVerbose())
144 {
145 G4cout << "### DNA models type " << typesDNA[i]
146 << " are activated for G4Region " << regnamesDNA[i] << G4endl;
147 }
148 const G4Region* reg = G4EmUtility::FindRegion(regnamesDNA[i], verbose);
149 if(nullptr == reg) { continue; }
150 G4int opt = 0;
151 if(typesDNA[i] == "DNA_Opt1") {
152 opt = 1;
153 } else if(typesDNA[i] == "DNA_Opt2") {
154 opt = 2;
155 } else if(typesDNA[i] == "DNA_Opt3") {
156 opt = 3;
157 } else if(typesDNA[i] == "DNA_Opt4") {
158 opt = 4;
159 } else if(typesDNA[i] == "DNA_Opt5") {
160 opt = 4;
161 } else if(typesDNA[i] == "DNA_Opt6") {
162 opt = 6;
163 } else if(typesDNA[i] == "DNA_Opt7") {
164 opt = 6;
165 } else if(typesDNA[i] == "DNA_Opt8") {
166 opt = 8;
167 }
168 DeactivateElectronProcesses(emaxDNA, emax, reg);
169 G4EmDNABuilder::ConstructDNAElectronPhysics(emaxDNA, opt, fast, st, reg);
170 DeactivateHadronProcesses(prot, emaxIonDNA, emax, reg);
171 G4EmDNABuilder::ConstructDNAProtonPhysics(eminBorn, emaxIonDNA, opt, fast, st, reg);
172 DeactivateIonProcesses(gion, emaxIonDNA, emax, reg);
173 G4EmDNABuilder::ConstructDNAIonPhysics(emaxIonDNA, st, reg);
174 DeactivateIonProcesses(alpha2, emaxLightIonDNA, emax, reg);
175 G4EmDNABuilder::ConstructDNALightIonPhysics(alpha2, 2, opt, emaxLightIonDNA, fast, st, reg);
176 DeactivateHadronProcesses(alpha1, emaxLightIonDNA, emax, reg);
177 G4EmDNABuilder::ConstructDNALightIonPhysics(alpha1, 1, opt, emaxLightIonDNA, fast, st, reg);
178 G4EmDNABuilder::ConstructDNALightIonPhysics(alpha0, 0, opt, emaxLightIonDNA, fast, st, reg);
179 G4EmDNABuilder::ConstructDNALightIonPhysics(h0, 0, opt, emaxIonDNA, fast, st, reg);
180 }
181 }
182
183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
184
185 void G4EmDNAPhysicsActivator::DeactivateElectronProcesses(const G4double emaxDNA,
186 const G4double emax,
187 const G4Region* reg)
188 {
189 if(emaxDNA >= emax) { return; }
190 const G4double msclimit = 100.*CLHEP::MeV;
191 G4ParticleDefinition* elec = G4Electron::Electron();
192 G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
193
194 G4VProcess* p;
195 if(emaxDNA < msclimit) {
196 p = G4PhysListUtil::FindProcess(elec, fMultipleScattering);
197 G4VMultipleScattering* msc = dynamic_cast<G4VMultipleScattering*>(p);
198 G4double elim = std::min(msclimit, emax);
199 if(nullptr == msc) {
200 msc = new G4eMultipleScattering();
201 ph->RegisterProcess(msc, elec);
202 }
203 auto mod = new G4GoudsmitSaundersonMscModel();
204 mod->SetActivationLowEnergyLimit(emaxDNA);
205 mod->SetHighEnergyLimit(elim);
206 msc->AddEmModel(-2, mod, reg);
207 }
208
209 p = G4PhysListUtil::FindProcess(elec, fIonisation);
210 G4VEnergyLossProcess* ptr = dynamic_cast<G4VEnergyLossProcess*>(p);
211 G4VEmFluctuationModel* fluc = nullptr;
212 if(nullptr == ptr) {
213 ptr = new G4eIonisation();
214 ph->RegisterProcess(ptr, elec);
215 }
216 auto modi = new G4MollerBhabhaModel();
217 modi->SetActivationLowEnergyLimit(emaxDNA);
218 modi->SetHighEnergyLimit(emax);
219 fluc = new G4UniversalFluctuation();
220 ptr->AddEmModel(-2, modi, fluc, reg);
221
222 p = G4PhysListUtil::FindProcess(elec, fBremsstrahlung);
223 ptr = dynamic_cast<G4VEnergyLossProcess*>(p);
224 if(nullptr == ptr) {
225 ptr = new G4eBremsstrahlung();
226 ph->RegisterProcess(ptr, elec);
227 }
228 auto modb = new G4SeltzerBergerModel();
229 modb->SetAngularDistribution(new G4Generator2BS());
230 modb->SetActivationLowEnergyLimit(emaxDNA);
231 modb->SetHighEnergyLimit(emax);
232 fluc = nullptr;
233 ptr->AddEmModel(-2, modb, fluc, reg);
234 }
235
236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
237
238 void
239 G4EmDNAPhysicsActivator::DeactivateHadronProcesses(G4ParticleDefinition* part,
240 const G4double emaxDNA,
241 const G4double emax,
242 const G4Region* reg)
243 {
244 if(emaxDNA >= emax) { return; }
245 G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
246 G4VProcess* p = G4PhysListUtil::FindProcess(part, fMultipleScattering);
247 G4VMultipleScattering* msc = dynamic_cast<G4VMultipleScattering*>(p);
248 if(nullptr == msc) {
249 msc = new G4hMultipleScattering();
250 ph->RegisterProcess(msc, part);
251 }
252 G4VMscModel* mod = new G4UrbanMscModel();
253 mod->SetActivationLowEnergyLimit(emaxDNA);
254 mod->SetHighEnergyLimit(emax);
255 msc->AddEmModel(-2, mod, reg);
256
257 const G4double braggmax = 2*CLHEP::MeV;
258 p = G4PhysListUtil::FindProcess(part, fIonisation);
259 G4VEnergyLossProcess* ptr = dynamic_cast<G4VEnergyLossProcess*>(p);
260 G4VEmFluctuationModel* fluc;
261 G4VEmModel* br;
262 if(part == G4GenericIon::GenericIon() || part == G4Alpha::Alpha()) {
263 br = new G4BraggIonModel();
264 fluc = new G4IonFluctuations();
265 } else {
266 br = new G4BraggModel();
267 fluc = new G4UniversalFluctuation();
268 }
269 if(nullptr == ptr) {
270 if(part == G4GenericIon::GenericIon() || part == G4Alpha::Alpha()) {
271 ptr = new G4ionIonisation();
272 } else {
273 ptr = new G4hIonisation();
274 }
275 ptr->SetFluctModel(fluc);
276 ph->RegisterProcess(ptr, part);
277 }
278 br->SetActivationLowEnergyLimit(emaxDNA);
279 br->SetHighEnergyLimit(braggmax);
280 ptr->AddEmModel(-2, br, fluc, reg);
281
282 auto be = new G4BetheBlochModel();
283 be->SetLowEnergyLimit(braggmax);
284 be->SetActivationLowEnergyLimit(braggmax);
285 be->SetHighEnergyLimit(emax);
286 ptr->AddEmModel(-3, be, fluc, reg);
287
288 DeactivateNuclearStopping(part, emaxDNA, reg);
289 }
290
291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
292
293 void
294 G4EmDNAPhysicsActivator::DeactivateIonProcesses(G4ParticleDefinition* part,
295 const G4double emaxDNA,
296 const G4double emax,
297 const G4Region* reg)
298 {
299 if(emaxDNA >= emax) { return; }
300 G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
301 G4VProcess* p = G4PhysListUtil::FindProcess(part, fMultipleScattering);
302 G4VMultipleScattering* msc = dynamic_cast<G4VMultipleScattering*>(p);
303 if(nullptr == msc) {
304 msc = new G4hMultipleScattering();
305 ph->RegisterProcess(msc, part);
306 }
307 auto mod = new G4UrbanMscModel();
308 mod->SetActivationLowEnergyLimit(emaxDNA);
309 mod->SetHighEnergyLimit(emax);
310 msc->AddEmModel(-2, mod, reg);
311
312 const G4double braggmax = 2*CLHEP::MeV;
313 p = G4PhysListUtil::FindProcess(part, fIonisation);
314 G4VEnergyLossProcess* ptr = dynamic_cast<G4VEnergyLossProcess*>(p);
315 G4VEmFluctuationModel* fluc = new G4IonFluctuations();
316 if(nullptr == ptr) {
317 ptr = new G4ionIonisation();
318 ptr->SetFluctModel(fluc);
319 ph->RegisterProcess(ptr, part);
320 }
321 auto br = new G4BraggIonModel();
322 br->SetActivationLowEnergyLimit(emaxDNA);
323 br->SetHighEnergyLimit(braggmax);
324 ptr->AddEmModel(-2, br, fluc, reg);
325
326 auto be = new G4BetheBlochModel();
327 be->SetLowEnergyLimit(braggmax);
328 be->SetActivationLowEnergyLimit(braggmax);
329 be->SetHighEnergyLimit(emax);
330 ptr->AddEmModel(-3, be, fluc, reg);
331
332 DeactivateNuclearStopping(part, emaxDNA, reg);
333 }
334
335 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
336
337 void
338 G4EmDNAPhysicsActivator::DeactivateNuclearStopping(const G4ParticleDefinition* part,
339 const G4double emax,
340 const G4Region* reg)
341 {
342 G4VProcess* p = G4PhysListUtil::FindProcess(part, fNuclearStopping);
343 G4NuclearStopping* ptr = dynamic_cast<G4NuclearStopping*>(p);
344 if(nullptr != ptr) {
345 auto mod = new G4ICRU49NuclearStoppingModel();
346 mod->SetActivationLowEnergyLimit(emax);
347 ptr->AddEmModel(-2, mod, reg);
348 }
349 }
350
351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
352