| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27
28 #include "G4DNADingfelderChargeDecreaseModel.h
29 #include "G4PhysicalConstants.hh"
30 #include "G4SystemOfUnits.hh"
31 #include "G4DNAMolecularMaterial.hh"
32 #include "G4DNAChemistryManager.hh"
33 #include "G4Log.hh"
34 #include "G4Pow.hh"
35 #include "G4Alpha.hh"
36
37
38 static G4Pow * gpow = G4Pow::GetInstance();
39
40 //....oooOO0OOooo........oooOO0OOooo........oo
41
42 using namespace std;
43
44 //....oooOO0OOooo........oooOO0OOooo........oo
45
46 G4DNADingfelderChargeDecreaseModel::G4DNADingf
47
48 G4VEmModel(nam)
49 {
50 numberOfPartialCrossSections[0] = 0;
51 numberOfPartialCrossSections[1] = 0;
52 numberOfPartialCrossSections[2] = 0;
53
54 verboseLevel = 0;
55 // Verbosity scale:
56 // 0 = nothing
57 // 1 = warning for energy non-conservation
58 // 2 = details of energy budget
59 // 3 = calculation of cross sections, file o
60 // 4 = entering in methods
61
62 if (verboseLevel > 0)
63 {
64 G4cout << "Dingfelder charge decrease mode
65 }
66 // Selection of stationary mode
67
68 statCode = false;
69 }
70
71 //....oooOO0OOooo........oooOO0OOooo........oo
72
73 void G4DNADingfelderChargeDecreaseModel::Initi
74
75 {
76
77 if (verboseLevel > 3)
78 {
79 G4cout << "Calling G4DNADingfelderChargeDe
80 << G4endl;
81 }
82
83 // Energy limits
84
85 G4DNAGenericIonsManager *instance;
86 instance = G4DNAGenericIonsManager::Instance
87 protonDef = G4Proton::ProtonDefinition();
88 alphaPlusPlusDef = G4Alpha::Alpha();
89 alphaPlusDef = instance->GetIon("alpha+");
90 hydrogenDef = instance->GetIon("hydrogen");
91 heliumDef = instance->GetIon("helium");
92
93 G4String proton;
94 G4String alphaPlusPlus;
95 G4String alphaPlus;
96
97 // Limits
98
99 proton = protonDef->GetParticleName();
100 lowEnergyLimit[proton] = 100. * eV;
101 highEnergyLimit[proton] = 100. * MeV;
102
103 alphaPlusPlus = alphaPlusPlusDef->GetParticl
104 lowEnergyLimit[alphaPlusPlus] = 1. * keV;
105 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
106
107 alphaPlus = alphaPlusDef->GetParticleName();
108 lowEnergyLimit[alphaPlus] = 1. * keV;
109 highEnergyLimit[alphaPlus] = 400. * MeV;
110
111 //
112
113 if (particle==protonDef)
114 {
115 SetLowEnergyLimit(lowEnergyLimit[proton]);
116 SetHighEnergyLimit(highEnergyLimit[proton]
117 }
118
119 if (particle==alphaPlusPlusDef)
120 {
121 SetLowEnergyLimit(lowEnergyLimit[alphaPlus
122 SetHighEnergyLimit(highEnergyLimit[alphaPl
123 }
124
125 if (particle==alphaPlusDef)
126 {
127 SetLowEnergyLimit(lowEnergyLimit[alphaPlus
128 SetHighEnergyLimit(highEnergyLimit[alphaPl
129 }
130
131 // Final state
132
133 //PROTON
134 f0[0][0]=1.;
135 a0[0][0]=-0.180;
136 a1[0][0]=-3.600;
137 b0[0][0]=-18.22;
138 b1[0][0]=-1.997;
139 c0[0][0]=0.215;
140 d0[0][0]=3.550;
141 x0[0][0]=3.450;
142 x1[0][0]=5.251;
143
144 numberOfPartialCrossSections[0] = 1;
145
146 //ALPHA++
147 f0[0][1]=1.; a0[0][1]=0.95;
148 a1[0][1]=-2.75;
149 b0[0][1]=-23.00;
150 c0[0][1]=0.215;
151 d0[0][1]=2.95;
152 x0[0][1]=3.50;
153
154 f0[1][1]=1.;
155 a0[1][1]=0.95;
156 a1[1][1]=-2.75;
157 b0[1][1]=-23.73;
158 c0[1][1]=0.250;
159 d0[1][1]=3.55;
160 x0[1][1]=3.72;
161
162 x1[0][1]=-1.;
163 b1[0][1]=-1.;
164
165 x1[1][1]=-1.;
166 b1[1][1]=-1.;
167
168 numberOfPartialCrossSections[1] = 2;
169
170 // ALPHA+
171 f0[0][2]=1.;
172 a0[0][2]=0.65;
173 a1[0][2]=-2.75;
174 b0[0][2]=-21.81;
175 c0[0][2]=0.232;
176 d0[0][2]=2.95;
177 x0[0][2]=3.53;
178
179 x1[0][2]=-1.;
180 b1[0][2]=-1.;
181
182 numberOfPartialCrossSections[2] = 1;
183
184 //
185
186 if( verboseLevel>0 )
187 {
188 G4cout << "Dingfelder charge decrease mode
189 << "Energy range: "
190 << LowEnergyLimit() / keV << " keV - "
191 << HighEnergyLimit() / MeV << " MeV for "
192 << particle->GetParticleName()
193 << G4endl;
194 }
195
196 // Initialize water density pointer
197 fpMolWaterDensity = G4DNAMolecularMaterial::
198
199 if (isInitialised)
200 { return;}
201 fParticleChangeForGamma = GetParticleChangeF
202 isInitialised = true;
203
204 }
205
206 //....oooOO0OOooo........oooOO0OOooo........oo
207
208 G4double G4DNADingfelderChargeDecreaseModel::C
209
210
211
212
213 {
214 if (verboseLevel > 3)
215 {
216 G4cout
217 << "Calling CrossSectionPerVolume() of
218 << G4endl;
219 }
220
221 // Calculate total cross section for model
222 if (
223 particleDefinition != protonDef
224 &&
225 particleDefinition != alphaPlusPlusDef
226 &&
227 particleDefinition != alphaPlusDef
228 )
229
230 return 0;
231
232 G4double lowLim = 0;
233 G4double highLim = 0;
234 G4double crossSection = 0.;
235
236 G4double waterDensity = (*fpMolWaterDensity)
237
238 const G4String& particleName = particleDefin
239
240 std::map< G4String,G4double,std::less<G4Stri
241 pos1 = lowEnergyLimit.find(particleName);
242
243 if (pos1 != lowEnergyLimit.end())
244 {
245 lowLim = pos1->second;
246 }
247
248 std::map< G4String,G4double,std::less<G4Stri
249 pos2 = highEnergyLimit.find(particleName);
250
251 if (pos2 != highEnergyLimit.end())
252 {
253 highLim = pos2->second;
254 }
255
256 if (k >= lowLim && k <= highLim)
257 {
258 crossSection = Sum(k,particleDefinition);
259 }
260
261 if (verboseLevel > 2)
262 {
263 G4cout << "_______________________________
264 G4cout << "G4DNADingfelderChargeDecreaeMod
265 G4cout << "Kinetic energy(eV)=" << k/eV <<
266 G4cout << "Cross section per water molecul
267 G4cout << "Cross section per water molecul
268 waterDensity/(1./cm) << G4endl;
269 // material->GetAtomicNumDensityVector()[1
270 }
271
272 return crossSection*waterDensity;
273 //return crossSection*material->GetAtomicNum
274
275 }
276
277 //....oooOO0OOooo........oooOO0OOooo........oo
278
279 void G4DNADingfelderChargeDecreaseModel::Sampl
280
281
282
283
284
285 {
286 if (verboseLevel > 3)
287 {
288 G4cout
289 << "Calling SampleSecondaries() of G4D
290 << G4endl;
291 }
292
293 G4double inK = aDynamicParticle->GetKineticE
294
295 G4ParticleDefinition* definition = aDynamicP
296
297 G4double particleMass = definition->GetPDGMa
298
299 G4int finalStateIndex = RandomSelect(inK,def
300
301 G4int n = NumberOfFinalStates(definition, fi
302 G4double waterBindingEnergy = WaterBindingEn
303 G4double outgoingParticleBindingEnergy = Out
304
305 G4double outK = 0.;
306
307 if (definition==G4Proton::Proton())
308 {
309 if (!statCode) outK = inK - n*(inK*electron
310 - waterBindingEne
311 else outK = inK;
312 }
313
314 else
315 {
316 if (!statCode) outK = inK - n*(inK*electron
317 - waterBindingEne
318 else outK = inK;
319 }
320
321 if (outK<0)
322 {
323 G4Exception("G4DNADingfelderChargeDecrease
324 FatalException,"Final kinetic energy i
325 }
326
327 fParticleChangeForGamma->ProposeTrackStatus(
328
329 if (!statCode) fParticleChangeForGamma->Prop
330
331 else
332
333 {
334 if (definition==G4Proton::Proton())
335 fParticleChangeForGamma->ProposeLocalEner
336 + waterBindingEnergy - outgoingParticleBi
337 else
338 fParticleChangeForGamma->ProposeLocalEner
339 + waterBindingEnergy - outgoingParticleBi
340 }
341
342 auto dp = new G4DynamicParticle (OutgoingPa
343 aDynamicParticle->GetMomentumDirection()
344 outK);
345 fvect->push_back(dp);
346
347 const G4Track * theIncomingTrack = fParticle
348 G4DNAChemistryManager::Instance()->Creat
349 1,
350 theIncomingTrack);
351
352 }
353
354 //....oooOO0OOooo........oooOO0OOooo........oo
355
356 G4int G4DNADingfelderChargeDecreaseModel::Numb
357
358
359 {
360 if (particleDefinition == G4Proton::Proton()
361 return 1;
362
363 if (particleDefinition == alphaPlusPlusDef)
364 {
365 if (finalStateIndex == 0)
366 return 1;
367 return 2;
368 }
369
370 if (particleDefinition == alphaPlusDef)
371 return 1;
372
373 return 0;
374 }
375
376 //....oooOO0OOooo........oooOO0OOooo........oo
377
378 G4ParticleDefinition* G4DNADingfelderChargeDec
379
380 {
381 if (particleDefinition == G4Proton::Proton()
382 return hydrogenDef;
383
384 if (particleDefinition == alphaPlusPlusDef)
385 {
386 if (finalStateIndex == 0)
387 return alphaPlusDef;
388 return heliumDef;
389 }
390
391 if (particleDefinition == alphaPlusDef)
392 return heliumDef;
393
394 return nullptr;
395 }
396
397 //....oooOO0OOooo........oooOO0OOooo........oo
398
399 G4double G4DNADingfelderChargeDecreaseModel::W
400
401 {
402 // Ionization energy of first water shell
403 // Rad. Phys. Chem. 59 p.267 by Dingf. et al
404 // W + 10.79 eV -> W+ + e-
405
406 if (particleDefinition == G4Proton::Proton()
407 return 10.79 * eV;
408
409 if (particleDefinition == alphaPlusPlusDef)
410 {
411 // Binding energy for W+ -> W++ + e-
412 // Binding energy for W -> W+ + e-
413 //
414 // Ionization energy of first water shell
415 // Rad. Phys. Chem. 59 p.267 by Dingf. et
416
417 if (finalStateIndex == 0)
418 return 10.79 * eV;
419
420 return 10.79 * 2 * eV;
421 }
422
423 if (particleDefinition == alphaPlusDef)
424 {
425 // Binding energy for W+ -> W++ + e-
426 // Binding energy for W -> W+ + e-
427 //
428 // Ionization energy of first water shell
429 // Rad. Phys. Chem. 59 p.267 by Dingf. et
430
431 return 10.79 * eV;
432 }
433
434 return 0;
435 }
436
437 //....oooOO0OOooo........oooOO0OOooo........oo
438
439 G4double G4DNADingfelderChargeDecreaseModel::O
440
441 {
442 if (particleDefinition == G4Proton::Proton()
443 return 13.6 * eV;
444
445 if (particleDefinition == alphaPlusPlusDef)
446 {
447 // Binding energy for He+ -> He++ + e-
448 // Binding energy for He -> He+ + e-
449
450 if (finalStateIndex == 0)
451 return 54.509 * eV;
452
453 return (54.509 + 24.587) * eV;
454 }
455
456 if (particleDefinition == alphaPlusDef)
457 {
458 // Binding energy for He+ -> He++ + e-
459 // Binding energy for He -> He+ + e-
460
461 return 24.587 * eV;
462 }
463
464 return 0;
465 }
466
467 //....oooOO0OOooo........oooOO0OOooo........oo
468
469 G4double G4DNADingfelderChargeDecreaseModel::P
470
471
472 {
473 G4int particleTypeIndex = 0;
474
475 if (particleDefinition == protonDef)
476 particleTypeIndex = 0;
477
478 if (particleDefinition == alphaPlusPlusDef)
479 particleTypeIndex = 1;
480
481 if (particleDefinition == alphaPlusDef)
482 particleTypeIndex = 2;
483
484 //
485 // sigma(T) = f0 10 ^ y(log10(T/eV))
486 //
487 // / a0 x + b0 x
488 // |
489 // y(x) = < a0 x + b0 - c0 (x - x0)^d0 x
490 // |
491 // \ a1 x + b1 x
492 //
493 //
494 // f0, a0, a1, b0, b1, c0, d0, x0, x1 are pa
495 //
496 // f0 has been added to the code in order to
497 //
498 // From Rad. Phys. and Chem. 59 (2000) 255-2
499 // Inelastic-collision cross sections of liq
500 //
501
502 if (x1[index][particleTypeIndex] < x0[index]
503 {
504 //
505 // if x1 < x0 means that x1 and b1 will be
506 //
507 // x1 = x0 + ((a0 - a1)/(c0 * d0)) ^ (1 /
508 //
509 // b1 = (a0 - a1) * x1 + b0 - c0 * (x1 - x
510 //
511
512 x1[index][particleTypeIndex] = x0[index][p
513 + gpow->powA((a0[index][particleTypeIn
514 / (c0[index][particleTy
515 * d0[index][particl
516 1. / (d0[index][particleTyp
517 b1[index][particleTypeIndex] = (a0[index][
518 - a1[index][particleTypeIndex]) * x1[i
519 + b0[index][particleTypeIndex]
520 - c0[index][particleTypeIndex]
521 * gpow->powA(x1[index][particleTyp
522 - x0[index][particl
523 d0[index][particleTypeI
524 }
525
526 G4double x(G4Log(k / eV)/gpow->logZ(10));
527 G4double y;
528
529 if (x < x0[index][particleTypeIndex])
530 y = a0[index][particleTypeIndex] * x + b0[
531 else if (x < x1[index][particleTypeIndex])
532 y = a0[index][particleTypeIndex] * x + b0[
533 - c0[index][particleTypeIndex]
534 * gpow->powA(x - x0[index][particl
535 d0[index][particleTypeI
536 else
537 y = a1[index][particleTypeIndex] * x + b1[
538
539 return f0[index][particleTypeIndex] * gpow->
540
541 }
542
543 G4int G4DNADingfelderChargeDecreaseModel::Rand
544
545 {
546 G4int particleTypeIndex = 0;
547
548 if (particleDefinition == protonDef)
549 particleTypeIndex = 0;
550
551 if (particleDefinition == alphaPlusPlusDef)
552 particleTypeIndex = 1;
553
554 if (particleDefinition == alphaPlusDef)
555 particleTypeIndex = 2;
556
557 const G4int n = numberOfPartialCrossSections
558 auto values(new G4double[n]);
559 G4double value(0);
560 G4int i = n;
561
562 while (i > 0)
563 {
564 i--;
565 values[i] = PartialCrossSection(k, i, part
566 value += values[i];
567 }
568
569 value *= G4UniformRand();
570
571 i = n;
572 while (i > 0)
573 {
574 i--;
575
576 if (values[i] > value)
577 break;
578
579 value -= values[i];
580 }
581
582 delete[] values;
583
584 return i;
585 }
586
587 //....oooOO0OOooo........oooOO0OOooo........oo
588
589 G4double G4DNADingfelderChargeDecreaseModel::S
590
591 {
592 G4int particleTypeIndex = 0;
593
594 if (particleDefinition == protonDef)
595 particleTypeIndex = 0;
596
597 if (particleDefinition == alphaPlusPlusDef)
598 particleTypeIndex = 1;
599
600 if (particleDefinition == alphaPlusDef)
601 particleTypeIndex = 2;
602
603 G4double totalCrossSection = 0.;
604
605 for (G4int i = 0; i < numberOfPartialCrossSe
606 {
607 totalCrossSection += PartialCrossSection(k
608 }
609
610 return totalCrossSection;
611 }
612
613