Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DNADingfelderChargeDecreaseModel.cc
Go to the documentation of this file.
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
30#include "G4SystemOfUnits.hh"
33
34//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
35
36using namespace std;
37
38//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
39
41 const G4String& nam) :
42G4VEmModel(nam), isInitialised(false)
43{
44 fpMolWaterDensity = 0;
45 numberOfPartialCrossSections[0] = 0;
46 numberOfPartialCrossSections[1] = 0;
47 numberOfPartialCrossSections[2] = 0;
48
49 verboseLevel = 0;
50 // Verbosity scale:
51 // 0 = nothing
52 // 1 = warning for energy non-conservation
53 // 2 = details of energy budget
54 // 3 = calculation of cross sections, file openings, sampling of atoms
55 // 4 = entering in methods
56
57 if (verboseLevel > 0)
58 {
59 G4cout << "Dingfelder charge decrease model is constructed " << G4endl;
60 }
62
63 // Selection of stationary mode
64
65 statCode = false;
66}
67
68//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69
71{
72}
73
74//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75
77 const G4DataVector& /*cuts*/)
78{
79
80 if (verboseLevel > 3)
81 {
82 G4cout << "Calling G4DNADingfelderChargeDecreaseModel::Initialise()"
83 << G4endl;
84 }
85
86 // Energy limits
87
91 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
92 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
93
94 G4String proton;
95 G4String alphaPlusPlus;
96 G4String alphaPlus;
97
98 // Limits
99
100 proton = protonDef->GetParticleName();
101 lowEnergyLimit[proton] = 100. * eV;
102 highEnergyLimit[proton] = 100. * MeV;
103
104 alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
105 lowEnergyLimit[alphaPlusPlus] = 1. * keV;
106 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
107
108 alphaPlus = alphaPlusDef->GetParticleName();
109 lowEnergyLimit[alphaPlus] = 1. * keV;
110 highEnergyLimit[alphaPlus] = 400. * MeV;
111
112 //
113
114 if (particle==protonDef)
115 {
116 SetLowEnergyLimit(lowEnergyLimit[proton]);
117 SetHighEnergyLimit(highEnergyLimit[proton]);
118 }
119
120 if (particle==alphaPlusPlusDef)
121 {
122 SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
123 SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
124 }
125
126 if (particle==alphaPlusDef)
127 {
128 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
129 SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
130 }
131
132 // Final state
133
134 //PROTON
135 f0[0][0]=1.;
136 a0[0][0]=-0.180;
137 a1[0][0]=-3.600;
138 b0[0][0]=-18.22;
139 b1[0][0]=-1.997;
140 c0[0][0]=0.215;
141 d0[0][0]=3.550;
142 x0[0][0]=3.450;
143 x1[0][0]=5.251;
144
145 numberOfPartialCrossSections[0] = 1;
146
147 //ALPHA++
148 f0[0][1]=1.; a0[0][1]=0.95;
149 a1[0][1]=-2.75;
150 b0[0][1]=-23.00;
151 c0[0][1]=0.215;
152 d0[0][1]=2.95;
153 x0[0][1]=3.50;
154
155 f0[1][1]=1.;
156 a0[1][1]=0.95;
157 a1[1][1]=-2.75;
158 b0[1][1]=-23.73;
159 c0[1][1]=0.250;
160 d0[1][1]=3.55;
161 x0[1][1]=3.72;
162
163 x1[0][1]=-1.;
164 b1[0][1]=-1.;
165
166 x1[1][1]=-1.;
167 b1[1][1]=-1.;
168
169 numberOfPartialCrossSections[1] = 2;
170
171 // ALPHA+
172 f0[0][2]=1.;
173 a0[0][2]=0.65;
174 a1[0][2]=-2.75;
175 b0[0][2]=-21.81;
176 c0[0][2]=0.232;
177 d0[0][2]=2.95;
178 x0[0][2]=3.53;
179
180 x1[0][2]=-1.;
181 b1[0][2]=-1.;
182
183 numberOfPartialCrossSections[2] = 1;
184
185 //
186
187 if( verboseLevel>0 )
188 {
189 G4cout << "Dingfelder charge decrease model is initialized " << G4endl
190 << "Energy range: "
191 << LowEnergyLimit() / keV << " keV - "
192 << HighEnergyLimit() / MeV << " MeV for "
193 << particle->GetParticleName()
194 << G4endl;
195 }
196
197 // Initialize water density pointer
199
200 if (isInitialised)
201 { return;}
203 isInitialised = true;
204
205}
206
207//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
208
210 const G4ParticleDefinition* particleDefinition,
211 G4double k,
212 G4double,
213 G4double)
214{
215 if (verboseLevel > 3)
216 {
217 G4cout
218 << "Calling CrossSectionPerVolume() of G4DNADingfelderChargeDecreaseModel"
219 << G4endl;
220 }
221
222 // Calculate total cross section for model
223
224 G4DNAGenericIonsManager *instance;
226
227 if (
228 particleDefinition != G4Proton::ProtonDefinition()
229 &&
230 particleDefinition != instance->GetIon("alpha++")
231 &&
232 particleDefinition != instance->GetIon("alpha+")
233 )
234
235 return 0;
236
237 G4double lowLim = 0;
238 G4double highLim = 0;
239 G4double crossSection = 0.;
240
241 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
242
243 const G4String& particleName = particleDefinition->GetParticleName();
244
245 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
246 pos1 = lowEnergyLimit.find(particleName);
247
248 if (pos1 != lowEnergyLimit.end())
249 {
250 lowLim = pos1->second;
251 }
252
253 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
254 pos2 = highEnergyLimit.find(particleName);
255
256 if (pos2 != highEnergyLimit.end())
257 {
258 highLim = pos2->second;
259 }
260
261 if (k >= lowLim && k <= highLim)
262 {
263 crossSection = Sum(k,particleDefinition);
264 }
265
266 if (verboseLevel > 2)
267 {
268 G4cout << "_______________________________________" << G4endl;
269 G4cout << "G4DNADingfelderChargeDecreaeModel" << G4endl;
270 G4cout << "Kinetic energy(eV)=" << k/eV << "particle :" << particleName << G4endl;
271 G4cout << "Cross section per water molecule (cm^2)=" << crossSection/cm/cm << G4endl;
272 G4cout << "Cross section per water molecule (cm^-1)=" << crossSection*
273 waterDensity/(1./cm) << G4endl;
274 // material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
275 }
276
277 return crossSection*waterDensity;
278 //return crossSection*material->GetAtomicNumDensityVector()[1];
279
280}
281
282//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
283
285 G4DynamicParticle*>* fvect,
286 const G4MaterialCutsCouple* /*couple*/,
287 const G4DynamicParticle* aDynamicParticle,
288 G4double,
289 G4double)
290{
291 if (verboseLevel > 3)
292 {
293 G4cout
294 << "Calling SampleSecondaries() of G4DNADingfelderChargeDecreaseModel"
295 << G4endl;
296 }
297
298 G4double inK = aDynamicParticle->GetKineticEnergy();
299
300 G4ParticleDefinition* definition = aDynamicParticle->GetDefinition();
301
302 G4double particleMass = definition->GetPDGMass();
303
304 G4int finalStateIndex = RandomSelect(inK,definition);
305
306 G4int n = NumberOfFinalStates(definition, finalStateIndex);
307 G4double waterBindingEnergy = WaterBindingEnergyConstant(definition, finalStateIndex);
308 G4double outgoingParticleBindingEnergy = OutgoingParticleBindingEnergyConstant(definition, finalStateIndex);
309
310 G4double outK = 0.;
311
312 if (definition==G4Proton::Proton())
313 {
314 if (!statCode) outK = inK - n*(inK*electron_mass_c2/proton_mass_c2)
315 - waterBindingEnergy + outgoingParticleBindingEnergy;
316 else outK = inK;
317 }
318
319 else
320 {
321 if (!statCode) outK = inK - n*(inK*electron_mass_c2/particleMass)
322 - waterBindingEnergy + outgoingParticleBindingEnergy;
323 else outK = inK;
324 }
325
326 if (outK<0)
327 {
328 G4Exception("G4DNADingfelderChargeDecreaseModel::SampleSecondaries","em0004",
329 FatalException,"Final kinetic energy is negative.");
330 }
331
333
334 if (!statCode) fParticleChangeForGamma->ProposeLocalEnergyDeposit(waterBindingEnergy);
335
336 else
337
338 {
339 if (definition==G4Proton::Proton())
340 fParticleChangeForGamma->ProposeLocalEnergyDeposit(n*(inK*electron_mass_c2/proton_mass_c2)
341 + waterBindingEnergy - outgoingParticleBindingEnergy);
342 else
343 fParticleChangeForGamma->ProposeLocalEnergyDeposit(n*(inK*electron_mass_c2/particleMass)
344 + waterBindingEnergy - outgoingParticleBindingEnergy);
345 }
346
347 G4DynamicParticle* dp = new G4DynamicParticle (OutgoingParticleDefinition(definition, finalStateIndex),
348 aDynamicParticle->GetMomentumDirection(),
349 outK);
350 fvect->push_back(dp);
351
352 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
354 1,
355 theIncomingTrack);
356
357}
358
359//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
360
361G4int G4DNADingfelderChargeDecreaseModel::NumberOfFinalStates(G4ParticleDefinition* particleDefinition,
362 G4int finalStateIndex)
363
364{
365 if (particleDefinition == G4Proton::Proton())
366 return 1;
367
370
371 if (particleDefinition == instance->GetIon("alpha++"))
372 {
373 if (finalStateIndex == 0)
374 return 1;
375 return 2;
376 }
377
378 if (particleDefinition == instance->GetIon("alpha+"))
379 return 1;
380
381 return 0;
382}
383
384//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
385
386G4ParticleDefinition* G4DNADingfelderChargeDecreaseModel::OutgoingParticleDefinition(G4ParticleDefinition* particleDefinition,
387 G4int finalStateIndex)
388{
390
391 if (particleDefinition == G4Proton::Proton())
392 return instance->GetIon("hydrogen");
393
394 if (particleDefinition == instance->GetIon("alpha++"))
395 {
396 if (finalStateIndex == 0)
397 return instance->GetIon("alpha+");
398 return instance->GetIon("helium");
399 }
400
401 if (particleDefinition == instance->GetIon("alpha+"))
402 return instance->GetIon("helium");
403
404 return 0;
405}
406
407//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
408
409G4double G4DNADingfelderChargeDecreaseModel::WaterBindingEnergyConstant(G4ParticleDefinition* particleDefinition,
410 G4int finalStateIndex)
411{
412 // Ionization energy of first water shell
413 // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
414 // W + 10.79 eV -> W+ + e-
415
417
418 if (particleDefinition == G4Proton::Proton())
419 return 10.79 * eV;
420
421 if (particleDefinition == instance->GetIon("alpha++"))
422 {
423 // Binding energy for W+ -> W++ + e- 10.79 eV
424 // Binding energy for W -> W+ + e- 10.79 eV
425 //
426 // Ionization energy of first water shell
427 // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
428
429 if (finalStateIndex == 0)
430 return 10.79 * eV;
431
432 return 10.79 * 2 * eV;
433 }
434
435 if (particleDefinition == instance->GetIon("alpha+"))
436 {
437 // Binding energy for W+ -> W++ + e- 10.79 eV
438 // Binding energy for W -> W+ + e- 10.79 eV
439 //
440 // Ionization energy of first water shell
441 // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
442
443 return 10.79 * eV;
444 }
445
446 return 0;
447}
448
449//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
450
451G4double G4DNADingfelderChargeDecreaseModel::OutgoingParticleBindingEnergyConstant(G4ParticleDefinition* particleDefinition,
452 G4int finalStateIndex)
453{
455
456 if (particleDefinition == G4Proton::Proton())
457 return 13.6 * eV;
458
459 if (particleDefinition == instance->GetIon("alpha++"))
460 {
461 // Binding energy for He+ -> He++ + e- 54.509 eV
462 // Binding energy for He -> He+ + e- 24.587 eV
463
464 if (finalStateIndex == 0)
465 return 54.509 * eV;
466
467 return (54.509 + 24.587) * eV;
468 }
469
470 if (particleDefinition == instance->GetIon("alpha+"))
471 {
472 // Binding energy for He+ -> He++ + e- 54.509 eV
473 // Binding energy for He -> He+ + e- 24.587 eV
474
475 return 24.587 * eV;
476 }
477
478 return 0;
479}
480
481//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
482
483G4double G4DNADingfelderChargeDecreaseModel::PartialCrossSection(G4double k,
484 G4int index,
485 const G4ParticleDefinition* particleDefinition)
486{
487 G4int particleTypeIndex = 0;
488 G4DNAGenericIonsManager* instance;
490
491 if (particleDefinition == G4Proton::ProtonDefinition())
492 particleTypeIndex = 0;
493
494 if (particleDefinition == instance->GetIon("alpha++"))
495 particleTypeIndex = 1;
496
497 if (particleDefinition == instance->GetIon("alpha+"))
498 particleTypeIndex = 2;
499
500 //
501 // sigma(T) = f0 10 ^ y(log10(T/eV))
502 //
503 // / a0 x + b0 x < x0
504 // |
505 // y(x) = < a0 x + b0 - c0 (x - x0)^d0 x0 <= x < x1
506 // |
507 // \ a1 x + b1 x >= x1
508 //
509 //
510 // f0, a0, a1, b0, b1, c0, d0, x0, x1 are parameters that change for protons and helium (0, +, ++)
511 //
512 // f0 has been added to the code in order to manage partial (shell-dependent) cross sections (if no shell dependence is present. f0=1. Sum of f0 over the considered shells should give 1)
513 //
514 // From Rad. Phys. and Chem. 59 (2000) 255-275, M. Dingfelder et al.
515 // Inelastic-collision cross sections of liquid water for interactions of energetic proton
516 //
517
518 if (x1[index][particleTypeIndex] < x0[index][particleTypeIndex])
519 {
520 //
521 // if x1 < x0 means that x1 and b1 will be calculated with the following formula (this piece of code is run on all alphas and not on protons)
522 //
523 // x1 = x0 + ((a0 - a1)/(c0 * d0)) ^ (1 / (d0 - 1))
524 //
525 // b1 = (a0 - a1) * x1 + b0 - c0 * (x1 - x0) ^ d0
526 //
527
528 x1[index][particleTypeIndex] = x0[index][particleTypeIndex]
529 + std::pow((a0[index][particleTypeIndex] - a1[index][particleTypeIndex])
530 / (c0[index][particleTypeIndex]
531 * d0[index][particleTypeIndex]),
532 1. / (d0[index][particleTypeIndex] - 1.));
533 b1[index][particleTypeIndex] = (a0[index][particleTypeIndex]
534 - a1[index][particleTypeIndex]) * x1[index][particleTypeIndex]
535 + b0[index][particleTypeIndex]
536 - c0[index][particleTypeIndex]
537 * std::pow(x1[index][particleTypeIndex]
538 - x0[index][particleTypeIndex],
539 d0[index][particleTypeIndex]);
540 }
541
542 G4double x(std::log10(k / eV));
543 G4double y;
544
545 if (x < x0[index][particleTypeIndex])
546 y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex];
547 else if (x < x1[index][particleTypeIndex])
548 y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex]
549 - c0[index][particleTypeIndex]
550 * std::pow(x - x0[index][particleTypeIndex],
551 d0[index][particleTypeIndex]);
552 else
553 y = a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex];
554
555 return f0[index][particleTypeIndex] * std::pow(10., y) * m * m;
556
557}
558
559G4int G4DNADingfelderChargeDecreaseModel::RandomSelect(G4double k,
560 const G4ParticleDefinition* particleDefinition)
561{
562 G4int particleTypeIndex = 0;
563 G4DNAGenericIonsManager *instance;
565
566 if (particleDefinition == G4Proton::ProtonDefinition())
567 particleTypeIndex = 0;
568
569 if (particleDefinition == instance->GetIon("alpha++"))
570 particleTypeIndex = 1;
571
572 if (particleDefinition == instance->GetIon("alpha+"))
573 particleTypeIndex = 2;
574
575 const G4int n = numberOfPartialCrossSections[particleTypeIndex];
576 G4double* values(new G4double[n]);
577 G4double value(0);
578 G4int i = n;
579
580 while (i > 0)
581 {
582 i--;
583 values[i] = PartialCrossSection(k, i, particleDefinition);
584 value += values[i];
585 }
586
587 value *= G4UniformRand();
588
589 i = n;
590 while (i > 0)
591 {
592 i--;
593
594 if (values[i] > value)
595 break;
596
597 value -= values[i];
598 }
599
600 delete[] values;
601
602 return i;
603}
604
605//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
606
607G4double G4DNADingfelderChargeDecreaseModel::Sum(G4double k,
608 const G4ParticleDefinition* particleDefinition)
609{
610 G4int particleTypeIndex = 0;
611 G4DNAGenericIonsManager* instance;
613
614 if (particleDefinition == G4Proton::ProtonDefinition())
615 particleTypeIndex = 0;
616
617 if (particleDefinition == instance->GetIon("alpha++"))
618 particleTypeIndex = 1;
619
620 if (particleDefinition == instance->GetIon("alpha+"))
621 particleTypeIndex = 2;
622
623 G4double totalCrossSection = 0.;
624
625 for (G4int i = 0; i < numberOfPartialCrossSections[particleTypeIndex]; i++)
626 {
627 totalCrossSection += PartialCrossSection(k, i, particleDefinition);
628 }
629
630 return totalCrossSection;
631}
632
@ eIonizedMolecule
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
@ fStopAndKill
double G4double
Definition: G4Types.hh:83
int G4int
Definition: G4Types.hh:85
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:52
static G4DNAChemistryManager * Instance()
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4DNADingfelderChargeDecreaseModel(const G4ParticleDefinition *p=0, const G4String &nam="DNADingfelderChargeDecreaseModel")
static G4DNAGenericIonsManager * Instance(void)
G4ParticleDefinition * GetIon(const G4String &name)
const std::vector< G4double > * GetNumMolPerVolTableFor(const G4Material *) const
Retrieve a table of molecular densities (number of molecules per unit volume) in the G4 unit system f...
static G4DNAMolecularMaterial * Instance()
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
size_t GetIndex() const
Definition: G4Material.hh:258
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:651
const G4Track * GetCurrentTrack() const
const G4String & GetParticleName() const
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:87
static G4Proton * Proton()
Definition: G4Proton.cc:92
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:757
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:133
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:652
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:645
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:764
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)